{'fraction_non_alphanumeric': 0.027895059723575477, 'fraction_numerical': 0.011728194579769766, 'mean_word_length': 4.420970266040689, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Light Quark Mass Reweighting

We present a systematic study of the effectiveness of light quark mass reweighting. This method allows a single lattice QCD ensemble, generated with a specific value of the dynamical light quark mass, to be used to determine results for other, nearby light dynamical quark masses. We study two gauge field ensembles generated with 2+1 flavors of dynamical domain wall fermions with light quark masses m_l=0.02 (m_\pi=620 MeV) and m_l=0.01 (m_\pi=420 MeV). We reweight each ensemble to determine results which could be computed directly from the other and check the consistency of the reweighted results with the direct results. The large difference between the 0.02 and 0.01 light quark masses suggests that this is an aggressive application of reweighting as can be seen from fluctuations in the magnitude of the reweighting factor by four orders of magnitude. Never-the-less, a comparison of the reweighed topological charge, average plaquette, residual mass, pion mass, pion decay constant, and scalar correlator between these two ensembles shows agreement well described by the statistical errors. The issues of the effective number of configurations and finite sample size bias are discussed. An examination of the topological charge distribution implies that it is more favorable to reweight from heavier mass to lighter quark mass.


I. INTRODUCTION
Generating ensembles of gauge field configurations with dynamical quarks is usually the most expensive part of a lattice QCD calculation. Because of the renormalization of the quark mass, we are usually not able to identify in advance the input bare quark mass that will correspond to a target renormalized mass. Quark mass reweighting is a powerful technique that allows us to fine tune the sea quark masses to their physical or other desired values after an ensemble has been generated, avoiding the computationally expensive generation of new ensembles. As is explained below, a reweighted lattice quantity is computed by averaging the product of that quantity and a reweighting factor over the given ensemble.
The reweighting factor is chosen to reproduce the effects of the change in the action that would result from changing from the simulated to the reweighted quark mass.
Reweighting of the strange quark mass has been widely applied to accomplish such a fine tuning of the physical strange quark mass [1][2][3] or to study dependence on the dynamical strange quark mass [4]. Reweighting of the light quark mass has been used less frequently, presumably because the light quark masses used in most current lattice calculations are sufficiently heavy that reweighting methods may not be able to bridge the gap between the simulated and physical masses. However, Aoki et al. [5] have applied light quark mass reweighting to obtain results at the physical light quark mass. In the future, we expect there will be an increasing use of mass reweighting to adjust the light quark mass to its physical value. Another important application is the use of light quark mass reweighting to avoid difficulties caused by a small quark mass such as long autocorrelation times or a loss of numerical stability at small mass. For example in the case Wilson fermions, one intentionally performs a simulation at larger light quark mass and then reweights to the physical mass [6,7]. This can also potentially reduce the computational cost.
As an important and widely used technique in lattice QCD, quark mass reweighting has been carefully studied and examined, for example, using Wilson [6] or overlap [8] fermions.
However, to better understand how mass reweighting performs, we believe that addition study of the resulting statistical fluctuations and a systematic comparison of a reweighted ensemble with one directly generated without reweighting are needed. As has been emphasized in earlier work such as Refs. [6][7][8], a small overlap in configuration space between the simulated ensemble and the target ensemble can lead to large fluctuations in the reweighting factors and unreliable results because of the lack of sufficient statistics. There are two issues.
First, for a fixed sample size N, the reweighting factor is determined from what may be a biased estimator. (We will discuss this in detail in section IV). If the reweighting factors fluctuate too much, this bias may lead to a large systematic error. Second, without a sufficiently large number of configurations with non-negligible weight, the error estimate may not be reliable, or even when reliable, the estimated may be too large for the result to be useful.
In this paper, we provide a systematic check of the reweighting technique by using two ensembles generated with the Iwasaki gauge action and the Domain Wall Fermion (DWF) action for the quarks. In order to see the limitations of reweighting, we consider a case were the method is used to achieve a large change in fermion mass. Specifically we reweight the light input quark mass from 0.02 (m π = 620 MeV) to 0.01(m π = 420 MeV) and then check the reweighted results for many important physical quantities against those computed directly on the dynamical m l = 0.01 ensemble. We also do this in the reverse direction, reweighting the m l = 0.01 ensemble to 0.02. With the ensemble size up to 1000 configurations (a total of 10,000 molecular dynamics time units), we find that even for such a large range of quark mass reweighting, the reweighted results and those obtained directly on the dynamically generated ensemble are consistent, demonstrating the success of the reweighting technique.
We study the behavior of reweighting on both short and long distance quantities. The paper is organized as follows: We begin in Sec. II with a description of the reweighting method, a detailed study of the reweighting factors and an analysis of the effective number of configuration that results after reweighting. We then discuss the reweighting of the topological charge and susceptibility in Sec. III. In Sec. IV, we discuss the reweighting of a short distance physical quantity, the average plaquette, and the problem of reweighting bias.
A second short distance quantity which shows significant dynamical quark mass dependence, the residual mass, is examined in Sec. V as well as two important physical observables: the pion mass and decay constant,. Finally the effect of reweighting on the scalar correlation function is presented in Sec. VI and in Sec. VII we further discuss and summarize our results.

II. LIGHT QUARK MASS REWEIGHTING METHOD
We follow closely the strategy presented in Ref [6]. Reweighting an ensemble of configurations from the light sea quark mass m 1 to m 2 requires the evaluation of a reweighting factor w(U; m 1 , m 2 ) for each configuration on which measurements are to be performed. We can see that the expectation value for any operator O on an ensemble with light sea quark mass m 2 can be calculated from an ensemble generated with sea quark mass m 1 as follows.
One begins with the usual expression for the expectation value of the operator O computed for the case of a sea quark mass m 2 : Here D(m) represents the DWF Dirac operator for a fermion of mass m, m s is the strange quark mass and the determinant factors with mass m = 1.0 are the standard Pauli-Villars contributions which appear in the DWF formulation. The quantity Z 2 appearing in the denominator is the partition function for the sea quark mass m 2 : simply the integral in Eq. (1) with the factor in the numerator representing the operator O removed.
Next one multiplies and divides by factors of Z 1 and det{D † (U, m 1 )D(U, m 1 )} obtaining Here we have introduced the reweighting factor w(U; m 1 , m 2 ) for each configuration U: and performed a similar manipulation to write the ratio Z 1 /Z 2 as w(m 1 , m 2 ) 1 . In the above formulae, we have used the notation D(U, m) to make explicit the dependence of the DWF Dirac operator on the gauge field variable U. The reweighting factor w(U; m 1 , m 2 ) can be determined from the stochastic average: where the average appearing in Eq. (6) is to be performed over the stochastic variable ξ drawn from the random Gaussian distribution exp{−ξ † ξ}.
We have purposely arranged the product of the four D(U, m) operators in the exponent of Eq. (6) in a form that will minimize the cost of the needed Conjugate Gradient (CG) inversion. Let us call the exponent of that equation f (U; m 1 , m 2 ; ξ) so that w(U; m1, m2) = If m 2 is close to m 1 , then the solution for η will usually converge faster than will the solution for a combination such as [D(m 2 ) † D(m 2 )] −1 D(m 1 )ξ, which would have resulted from a different arrangement than that used above.
It was pointed out in Ref [6] that the average of stochastic sources and the average of gauge fields can be performed in any order, so we can perform the average needed in Eq. (6) using any number of random hits N hit . By doing so, we introduce extra noise into our system, but no systematic error. To reduce the noise following the idea of determinant factorization proposed by Hasenbush [6,9], we calculate the reweighting factor by carrying out N step steps. For each step we change the mass by ∆m = (m 2 − m 1 )/N step . Following this approach, the reweight factor for a given configuration is then calculated from the product where for each step, N hit hits are performed by summing over the N hit Gaussian variables Our goal is to calculate the reweighting factor with sufficient accuracy that the uncertainty introduced by the finite number of the random hits is smaller than the fluctuation of w(U; m1, m2) among the underlying gauge configurations. The total number of DWF Dirac inversions required by Eq. (9) is N t = N hit N step . By increasing N step , we reduce the fluctuation among the independent steps. Averaging over an additional N hit hits gives a better estimate of the average. From our experiments with a fixed number of total inversions, we find that smaller errors can be achieved from more steps than more hits. These results confirm similar conclusions reached by the RBC and UKQCD collaborations when reweighting in the strange quark mass.
This observation can be verified by the following explicit calculation. Consider a single gauge configuration and evaluate the reweighting factor w as a product of the reweighting factors w s , 1 ≤ s ≤ N step corresponding to N step steps. In our stochastic evaluation the factor w s corresponding to the step s is given by a function w s (ξ) of a set of random variables ξ where a different stochastic vector ξ will be used for each step. Define the quantity σ s from the ratio of averages: where the average is taken over the stochastic vector ξ in the fixed gauge background. The factor 1/N 2 step has been extracted from the exponent on the right hand side because the exponent in Eq. (6) will be reduced in size by 1/N step when a fixed mass interval is divided into N step segments. Since the reweighting factor w is linear in each of the independent stochastic factors w s (ξ s ), it is a straight forward exercise to express the fluctuations in w in terms of the quantities σ s . If w is obtained from a single hit: where in the final line we have made the simplifying assumption that σ s is independent of s.
If this process is repeated N hit times and the results averaged (N hit hits) then the resulting error from the stochastic sampling becomes A smaller error results if the N hit hits are averaged at each step, as is indicated in Eq. (9).
For this case a series of steps similar to those in Eqs. (11)- (15) gives: where we write the simpler expression that results if we continue to use uniform steps and assume that the corresponding quantities σ s are independent of s.
While the fluctuations given in Eq. (17) are smaller than those in Eq. (16), both expressions are decreasing functions of N step for a fixed total number N t of Dirac inversions, making it favorable to use larger N step instead of larger N hit . When N step is sufficient large, both of the above formulae become approximately wσ/ √ N t . In this limit it is unimportant how N t is divided between steps and hits. A second advantage of introducing a number of Hasenbusch steps is that we can then reweight to many different intermediate mass values in a single reweighting calculation. A third advantage can be seen from Eq. Because of the large shift in mass that we are attempting to achieve with reweighting, the reweighting factors presented in Fig. 1 vary in size by a few orders of magnitude. Such large fluctuations in the reweighting factors may cause concern that such a substantial reweighting is doomed to fail. However, the large fluctuations of these factors toward lower values are not important. No mater how small they are, their only consequence is that the corresponding configurations make a negligible contribution in the reweighted ensemble.
By simply counting the number of configurations that have reweighting factor larger than 1, we see a variation in size of 100 and 200 for the 0.02 and 0.01 ensembles respectively.
Thus, the reweighting may be more successful than Fig. 1 naively suggests. The efficiency of reweighting can be estimated in a more systematic way as follows.
Suppose that we have a specific ensemble of N configurations. The reweighted value for any physical quantity O based on Eq. (3) is where the quantities appearing on the right-hand side can be approximately calculated from a single ensemble of N configurations according to the equations: The quantities C(l) and W (l) appearing in Eq. (23) are the autocorrelation functions for the measured quantity, O i and the the reweighting factors w i and can be estimated from a single ensemble using the formulae: for the case l ≥ 0. When evaluating the autocorrelation time with a finite sample size N, we must impose a maximum length L max ≪ N. Here, we use the smallest value of l at which C(l + 1) becomes negative.
If we interpret the error given in Eq. (19) as the size of the fluctuations among our samples divided by the square root of an effective number of configurations, δO 2 /N eff , we can determine N eff to be: In the case of no reweighting(w i ≡ 1), this expression reverts to our normal expression for the effective number of configurations N/τ corr . Even though different physical quantities have different autocorrelation times, the reweighted auto correlation times τ corr in Eq. (23) are usually very close to 1 for our data, because C(l) and W (l) are both small numbers except when l=0. Therefore, the estimated effective number of configurations can be calculated from the reweighting factors alone. We find N eff = 48 for the m l = 0.02 ensemble (reweighting to 0.01) and N eff = 63 for the m l = 0.01 ensemble (reweighting to 0.02). Clearly, the effective sample size decreases dramatically because of our large range of mass reweighting. However, with a sample size of roughly 50 configurations, we can expect that the reweighting still works well for many interesting physical quantities.
To justify our assumption that τ corr ∼ 1, we calculated the autocorrelation function C(l) for a few of quantities that we will study later. Figure 2 shows the autocorrelation function  Table II. We can see that even for the topological charge which has the longest autocorrelation time, the reweighted integrated autocorrelation time is less than 2.
Since most quantities have significantly shorter autocorrelation time compared to that of the topological charge, we can expect that τ → 1, thus supporting our estimate of the number of effective configurations. As can be seen in Tab. II, the reweighted ensembles typically show a decreased autocorrelation time. This should be expected because reweighting effectively thins the initial ensemble. Note this decrease in τ corr partially compensates for the decrease in N eff arising from the ω 2 /ω 2 factor inEq. (26).

III. TOPOLOGICAL CHARGE AND SUSCEPTIBILITY
The topological charge typically has the longest autocorrelation time among the usual observables and represents long distance physics. We will examine the reweighting effect on it first. We will later study short distance quantities such as the average plaquette and residual mass.
As the first step in determining the topological charge we use Wilson flow [11,12] developed by Lüscher to smooth or cool the gauge configurations. As a function of the Wilson flow time, the gauge fields obey the following differential equation and initial conditions: where S w (U) = 1 g 2 0 p ReTr{1 − U(p)} is the Wilson gauge action and {T a } 1≤a≤8 the eight hermitian generators of SU(3). When acting on a particular link matrix U(y, ν) the matrix X a = T a if y = x and ν = µ and is zero otherwise. We follow the Runge-Kutta integration scheme (see Appendix C of Ref [12] ). This realizes a continuous version of the stout smearing technique developed by C. Morningstar and M. Peardon [13]. We have used an integration time step δt = 0.04, which gives an error smaller than 10 −5 for t ≤ 16 as seen by comparing the results to those obtained with δt = 0.005 for dozens of configurations. All results given here use δt = 0.04 and the integration error should be negligible for the topological charge calculation of interest.
We then determine the topological charge using the "5Li" (5-loop-improved) algorithm [14] on these "cooled" configurations, finding values of topological charge very close to integers. Figure 3 shows a few typical evolution curves for the topological charge Q as a function of the flow time t. Notice that there are ambiguities for some of the configurations: the topological charge converges to one integer, stays there for a while and then moves away converging to a neighboring integer. The problem of defining the topological charge on the lattice and its ambiguities are not the focus of this paper. We will simply adopt the Wilson flow prescription and use the value of Q at t = 16 as our result for the topological charge.  given by the probability function P w (Q) computed from the relation: where the Kronecker delta factor δ Q i ′ ,Q ensures that sum over i ′ only includes those configurations having topological charge Q.
The most important property of the topological charge distribution (given its central value of zero) is the width Q 2 of the distribution, which is related to the topological susceptibility. The calculated values of Q 2 with the average Q and the reweighted values Q 2 w and Q w are summarized below the distributions shown in Fig. 4. The quoted errors are calculated from the fluctuations seen among blocks of 100 trajectories to remove any effect of autocorrelation.
We can see very clearly how we have successfully reweighted one distribution to another.  Fig. 4 (bottom row, right-hand column) may look troublesome. We find that its contribution comes to a large extent from a single configuration that has the largest reweighting factor in Fig. 1. Therefore, this large spike has very large inherent uncertainly.
To evaluate the quality of the data set, we (incorrectly) dropped this configuration and performed a complete analysis, repeating all the calculations in this paper. This experiment revealed no surprises, giving results consistent with those presented here which included all configurations. Nonetheless, the unphysical topological charge distribution for the ensemble reweighted to m l = 0.02 ensemble suggests that our calculation is not far from being sensitive  both long distance and short distance physics.
So far we have ignored the bias which arises from the finite sample size used in our stochastic estimate of the reweighted quantity O w N given by the ratio in Eq. (18). Let us take the average plaquette as an example and see how large it is. Suppose the sample size is N. The reweighted result P w N is biased in the sense that the expectation value P w N will not agree with the true value P w = lim N →∞ P w N . It is straightforward to work out the difference, keeping only the lowest order terms: We can see that the bias is proportional to 1/N, and the coefficient is related to the correlation of w i and P i . We might attempt to estimate the bias given in this equation for our single sample by omitting the average over samples, . . . and by examining the l = i term which might be expected to be dominant. This term can be written as This expression hints that the bias is negative for a positive correlation between w i and P i , and vice versa. As should be expected, if there is no correlation between P i and w i , the coefficient vanishes. Taking all of our data and dividing them into small samples, we calculate the expectation value of P w N . The results are shown in Fig. 6 and suggest that the bias for the m l = 0.02 ensemble is negative and the bias for the m l = 0.01 ensemble is positive although these effects are at or below the one σ level. These results also suggest that the bias is negligible for reasonably large ensemble size N > 20.

V. PION MASS AND DECAY CONSTANT
Let us next study the effect of reweighting on some additional, interesting physical quantities. In this section, we will examine the sea quark dependence of the residual mass, the pion mass, and the pion decay constant. In the following, for the m l = 0.02(0.01) ensemble, we fix the valence quark mass at 0.01(0.02) and gradually reweight the sea quark mass To obtain m π and f π , we calculated three different correlation functions. Using a gaugefixed pseudo-scalar wall source, we choose the sink to be either a gauge-fixed pseudo-scalar wall sink, a pseudo-scalar point sink or a point sink formed from the zeroth component of the axial current. These correlation functions are identified as C sink 1 source 2 (t) which for our three cases becomes C W W P P (t), C LW P P (t) and C LW AP (t). Here the label 'W' means a gauge-fixed wall source or sink, 'L' means a point sink, 'P' means pseudo-scalar operator, 'A' means the zeroth component of the axial current operator. We perform a simultaneous fit to all three correlators at large times, using the functional form This fit give us m π directly. From the amplitudes N s 1 s 2 O 1 O 2 , we construct f π following Ref [2], where Z A is the renormalization constant of the local axial current A µ , calculated in Ref. [15].
Since Z A is to be evaluated in the chiral limit and thus does not depend on sea quark mass, we simply set Z A = 1 for the present calculation. The reweighted results show a clear trend, although with increasing errors, as we reweight the sea quark mass. In both cases the result obtained in the final step of the reweighting procedure agrees with the directly calculated result at the 1 σ level. Figure. 9 shows similar results for f π , presented in the same way as the data for m π . in Fig. 8. Similar conclusions can be drawn.

VI. SCALAR CORRELATOR
Finally, let us look at something less successful but equally interesting. The properties of the iso-vector, scalar particle are known to be difficult to calculate because its mass is much larger than the mass of the pion, so large are statistics are needed. For the reweighted data, the effective number of configurations drops to less than one hundred, resulting in such a small ensemble size that it is not possible to determine the mass of this scalar state.
However, by simply looking at the scalar correlator itself, we can see the dramatic effect of reweighting: the partially quenched violation of positivity disappears as the sea quark mass is reweighted to the unitary value.
We calculate the scalar correlator using an identical, gauge-fixed wall source and sink so that the resulting correlation function should be positive. A well known artifact of partial quenching appears in the case where m val < m sea : the correlator becomes negative for large source sink time separations when the lighter, negative norm πη ′ intermediate state dominates [16]. This behavior can be clearly seen in the left panel of Fig. 10 where the valence quark mass is 0.01, and the sea quark mass is 0.02. The correlator becomes negative for t > 4. After reweighting the sea quark mass from 0.02 down to the unitary value of 0.01, the correlator data points increase and becomes positive as seen previously in Ref. [6]. In

VII. SUMMARY AND DISCUSSION
We have shown how reweighting can be used to vary the light sea quark mass while working with a single ensemble generated with a single sea quark mass. We first successfully reweighted the topological charge distribution from one of our light sea quark masses to the another. We found that it is better to reweight from a heavier mass to a lighter mass, because this corresponds to reweighting from a wider distribution to one that is more narrow.
In this case, the poor statistics present in the tail of the distribution is suppressed instead of magnified by the reweighting factor. Next we studied the reweighting of the average plaquette, examining results obtained for a sequence of intermediate masses. In addition, we examined the bias introduced by our stochastic finite-sample reweighted measurement of the average plaquette. We found that the bias is negligible for fairly large ensemble with N > 20. After this study of both long distance and short distance quantities, we examined the behavior under reweighting of the residual mass, pion mass and pion decay constant. In each case, reweighting revealed how the given physical quantity varied with the sea quark mass when computed on an ensemble with a fixed quark ass. In each case, the final step of reweighting was chosen to duplicate the sea quark mass for which had been used to generate a second ensemble. This allowed a direct comparison between the reweighted and a directly generated result and in each case the two results agreed within errors. Finally, reweighting was shown to dramatically to correct the partial quenching artifacts seen in the calculation of the iso-scalar, scalar correlator initially computed with m val = 0.01 on an m l = 0.02 ensemble.
The error on the quantities we have studied was generally increased 3-4 times by reweighting, in rough agreement with our estimate of the number N eff of effective configurations. It should be emphasized that the change in the quark mass from m l = 0.02 to 0.01 and the reverse is a large change to be accomplished by reweighting. This is certainly suggested by the large fluctuations of more than 4 orders of magnitude seen in the reweight factors.
The resulting efficiency in the reweighting procedure, N eff /N, is less than 10%, causing the error after reweighting to increase by a factor of 4. However, such an aggressive use of reweighting may be appropriate for a test of the method, more clearly revealing its power and potential limitations.
In real applications of reweighting, we may choose to be more conservative and reweight over a smaller range of mass than considered here to keep the efficiency high. Had we only reweighted by half of the mass difference, reweighting both ensembles to m l = 0.015, the effective number of configurations would increase dramatically from 48 to 308 for the m l = 002 ensemble, and from 63 to 386 for the m l = 0.01 ensemble. Depending on the physics being studied, we must ensure that N eff is large enough to give statistically useful results. For example, N eff < 100, is apparently insufficient for a study of the iso-vector, scalar particle mass. In conclusion, reweighting has been shown to work remarkably well, even for a large change in the light quark mass, with the resulting errors well described by an effective number of configurations given in Eq. 26.