Variance Reduction and Cluster Decomposition

TL;DR

This paper introduces a method based on the cluster decomposition principle to reduce statistical noise in lattice QCD calculations, significantly improving the precision of hadron mass and CP violation measurements.

Contribution

It demonstrates that the cluster decomposition principle can effectively reduce noise in lattice QCD calculations, leading to more accurate results for various physical quantities.

Findings

Statistical errors reduced by a factor of 3 to 4.

Systematic errors comparable to statistical errors for strangeness content.

Overall error decreased by 2 to 3 times.

Abstract

It is a common problem in lattice QCD calculation of the mass of the hadron with an annihilation channel that the signal falls off in time while the noise remains constant. In addition, the disconnected insertion calculation of the three-point function and the calculation of the neutron electric dipole moment with the $\theta$ term suffer from a noise problem due to the $\sqrt{V}$ fluctuation. We identify these problems to have the same origin and the $\sqrt{V}$ problem can be overcome by utilizing the cluster decomposition principle. We demonstrate this by considering the calculations of the glueball mass, the strangeness content in the nucleon, and the CP violation angle in the nucleon due to the $\theta$ term. It is found that for lattices with physical sizes of 4.5 - 5.5 fm, the statistical errors of these quantities can be reduced by a factor of 3 to 4. The systematic errors can be estimated from the Akaike information criterion. For the strangeness content, we find that the systematic error is of the same size as that of the statistical one when the cluster decomposition principle is utilized. This results in a 2 to 3 times reduction in the overall error.

Tables

Table 1. Table 1: The mean values, and the systematic and statistical errors are given for various fits. N b subscript 𝑁 𝑏 N_{b} is the number of bootstrap samples, N d subscript 𝑁 𝑑 N_{d} the number is data points in each fit, and N c subscript 𝑁 𝑐 N_{c} is the number of different combinations of the N d subscript 𝑁 𝑑 N_{d} data points.

$N_b$	$N_d$	$N_c$	mean	$E_{\mathrm{sys}}$	$E_{\mathrm{sta}}$
100	6	80	0.161	0.013	0017
100	8	80	0.160	0.015	0.015
100	10	80	0.163	0.019	0.012
100	8	60	0.161	0.015	0.015
100	8	80	0.160	0.015	0.015
100	8	100	0.161	0.015	0.015
50	8	80	0.160	0.015	0.015
100	8	80	0.160	0.015	0.015
200	8	80	0.160	0.015	0.015

Equations

$$\begin{array}{rl}& |⟨0|{\mathcal{B}}_1(x_1){\mathcal{B}}_2(x_2)|0{⟩}_s|\le A{r}^{-\frac{2}{3}}{e}^{-Mr}\end{array}$$

$${\int }_0^R{d}^{3}r{r}^{-\frac{3}{2}}{e}^{-Mr}=4\pi \left(\frac{\sqrt{\pi }\mathrm{erf}(\sqrt{MR})}{2M^{\frac{3}{2}}}-\frac{\sqrt{R}{e}^{-MR}}{M}\right),$$

$$C(R,t)=\frac{1}{V}⟨\sum _{\vec{x}}\sum _{r<R}{\mathcal{O}}_1(\vec{x}+\overrightarrow{r^{\prime }},t){\mathcal{O}}_2(\vec{x},0)⟩,$$

$$Var(R,t)$$

$$\frac{S/N(R_s)}{S/N(L)}\sim \sqrt{\frac{V}{V_{R_s}}}.$$

$$C(R,t)=⟨{\int }_{r<R}d\overrightarrow{r^{\prime }}K(\overrightarrow{r^{\prime }},t)⟩.$$

$$C_3(R,\tau ,t)=⟨\sum _{\vec{x}}\sum _{r<R}{\mathcal{O}}_N(\vec{x},t)S(\vec{x}+\overrightarrow{r^{\prime }},\tau ){\bar{\mathcal{O}}}_N(\mathcal{G},0)⟩,$$

$${\alpha }^1=\frac{\mathrm{Tr}[C_{3\mathcal{Q}}(t){\gamma }_5]}{\mathrm{Tr}[C_2(t){\Gamma }_e]}$$

$$C_{3\mathcal{Q}}(t)=⟨\sum _{\vec{x}}{\mathcal{O}}_N(\vec{x},t){\bar{\mathcal{O}}}_N(\mathcal{G},0)\mathcal{Q}⟩.$$

$$\frac{S/N(R_s)}{S/N(L)}\sim \sqrt{\frac{A(L/2,3m_{\pi })}{A(R_s,3m_{\pi })}}.$$

$$L_{\max }=P(x|\vec{\theta },M),$$

$$L_{\max }=e^{-\frac{{\chi }_{\min }^2}{2}}.$$

$$\mathrm{AIC}=2k-2\ln (L_{\max })=2k+{\chi }_{\min }^2,$$

$$w_{i=}\frac{e^{-\frac{{\mathrm{AIC}}_i}{2}}}{{\sum }_i{e}^{-\frac{{\mathrm{AIC}}_i}{2}}},$$

$$\bar{x}=\sum _{i=1}^P{w}_i{x}_i,$$

$${\sigma }_{\mathrm{sys}}=\sqrt{\sum _{i=1}^P{w}_i{\sigma }_i^2},{\sigma }_i=(x_i-\bar{x}).$$

$$f_1(R)=C_0$$

$$f_2(R)=C_0+C_1\frac{\sqrt{R}{e}^{-mR}}{m}.$$

Full text

I Introduction

As the physical pion mass is accessible in lattice QCD simulations nowadays with larger physical volumes and several lattice spacings with different lattice actions to estimate the associated systematic errors, lattice QCD calculation is getting mature, particularly for flavor physics where the quark masses, heavy-light decay constants, CKM matrices, and strong coupling constant are reviewed and averaged by FLAG Aoki et al. (2017). On the other hand, the baryon physics is not as settled as that of mesons. Part of the reasons is illustrated in the Parisi-Lepage consideration of the signal-to-noise ratio of the nucleon two-point function. Since the variance of the nucleon propagator has three-pions as the lowest state in the correlator, the signal-to-noise (S/N) ratio is proportional to $e^{-(m_{N}-3/2m_{\pi})t}$ Parisi (1984); Lepage (1989) and noticeably grows exponentially with $t$ when the pion mass is close to the physical one in lattice calculations. This is why baryon physics is more noisy than that of mesons.

One special problem associated with the correlators of mesons involving annihilation channels or glueballs is that the signal falls off exponential with time, but the noise remains constant. Thus, after certain time separation, the signal falls below the noise and succumbs to the sign problem. Another aspect of the DI is observed in the DI three-point function involving a quark loop or the topological charge in the neutron electric dipole moment (nEDM) calculation from the $\theta$ term, the fluctuations of the quark loop and the topological charge are proportion to $\sqrt{V}$ which pose a challenge for calculations large volumes lattices. In this work, we shall show that the constant error and $\sqrt{V}$ fluctuation in the DI have the same origin and they can be ameliorated with the help of the property of the cluster decomposition principle so that the S/N ratio can be improved by a factor of $\sqrt{V/V_{R_{s}}}$ where $V_{R_{s}}$ is the volume with radius $R_{s}$ which is the effective correlation length between the operators.

II Cluster decomposition principle and variance reduction

One often invokes the locality argument to justify that experiments conducted on Earth is not affected by events on the Moon. This is a consequence of the cluster decomposition principle (CDP) in that if color-singlet operators in a correlator are separated by a large enough space-like distance, the correlator will be zero. In other words, the operators are not correlated in this circumstance. To be specific, it is shown Araki et al. (1962) that under the assumptions of translation invariance, stability of the vacuum, existence of a lowest non-zero mass and local commutativity, one has

[TABLE]

for a large enough space-like distance $r=|x_{1}-x_{2}|$, where $\langle 0|\mathcal{B}_{1}(x_{1})\mathcal{B}_{2}(x_{2})|0\rangle_{s}\equiv\langle 0|\mathcal{B}_{1}(x_{1})\mathcal{B}_{2}(x_{2})|0\rangle-\langle 0|\mathcal{B}_{1}(x_{1})|0\rangle\langle 0|\mathcal{B}_{2}(x_{2})|0\rangle$ is the vacuum-subtracted correlation function. $\mathcal{B}_{1}(x_{1})$ and $\mathcal{B}_{2}(x_{2})$ are two color-singlet operator clusters whose centers are at $x_{1}$ and $x_{2}$ respectively, $M$ is the smallest non-zero mass for the correlator, and $A$ is a constant. This is the asymptotic behavior of a boson propagator $K_{1}(r)/r$. This means the correlation between two operator clusters far apart with large enough space-like distance $r$ tends to be zero at least as fast as $r^{-\frac{2}{3}}e^{-Mr}$. Given that the longest correlation length in QCD is $1/m_{\pi}$, one has $M\geq m_{\pi}$. Since the Euclidean separation is always ‘space-like’, the cluster decomposition principle (CDP) is applicable to the Euclidean correlators. Some of the recent attempts to reduce variances in the calculations of strangeness in the nucleon Bali et al. (2009), the $\rho$ meson mass Chen et al. (2015), the light-by-light contribution in the muonic $g-2$ Blum et al. (2016a), the factorization of fermion determinant Cè et al. (2017), and reweighting of nEDM calculation with topological charge density Shintani et al. (2016) have applied concepts similar or related to that of CDP. In this work, we prove that, applying the CDP explicitly, the error of an DI correlator can be improved by a factor of $\sqrt{V/V_{R_{s}}}$.

In evaluating the correlators, one often takes a volume sum over the three-dimensional coordinates. To estimate at what distance the large distance behavior saturates, we integrate the fall-off to a cut off distance $R$,

[TABLE]

where ${\rm erf}$ is the error function. Since the kernel of the integral decays very quickly, the integral has already gained more than $99.5\%$ of its total value for $R=8/M$. Assuming the fall-off behavior dominates the volume-integrated correlator, we consider $R_{s}\sim\frac{8}{M}$ as an effective cutoff and the correlation with separation $r>R_{s}$ has negligible contribution.

To test the principle of cluster decomposition with lattice data, we consider the two-point correlator for a fixed $t$ with a cutoff of $R$ in the relative coordinate between the two color-singlet operators $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$

[TABLE]

where $r=\sqrt{|\vec{r^{\prime}}|^{2}+t^{2}}$. The correlation functions in the present work are calculated using valence overlap fermions on the RBC-UKQCD $2+1$ flavor domain-wall configurations. More detailed definitions and numerical implementations can be found in previous works  Li et al. (2010); Gong et al. (2013); Sufian et al. (2017); Yang et al. (2016a).

We examine the nucleon two-point function first on the $48^{3}\times 96$ lattice (48I) with the physical sea quark mass  Blum et al. (2016b). We use 3 valence quark masses corresponding to pion masses $70$ MeV, $149$ MeV and $260$ MeV respectively.

The results for the nucleon correlators at $t=9$ for three different valence quark masses are plotted in Fig. 1 as a function of $R$ which is the cutoff of the Euclidean distance $r$ between the point source and the sink. We see that the nucleon correlator basically saturates after $R\sim 15=1.71$ fm with $a=0.114$ fm for the three cases. This agrees well with our earlier estimate of a saturation radius $R_{s}=8/M$ which corresponds to $\sim 1.66$ fm. This shows that the CDP works and Eq. (2) gives a good estimate of $R_{s}$.

Since the signal of the correlators falls off exponentially with $r$, summing over $r$ beyond the saturating radius $R_{s}$ does not change the signal and will only gather noise. Let’s consider the disconnected insertion next and see how the S/N ratio can be improved with this observation. In the case of the DI, the variance of the correlator in Eq. (3) $\frac{1}{V^{2}}\langle|\sum_{\vec{x}}\sum_{r<R}\mathcal{O}_{1}(\vec{x}+\vec{r^{\prime}},t)\,\mathcal{O}_{2}(\vec{x},0)|^{2}\rangle$ can have a vacuum insertion in addition to the exponential fall off in $t$ due to the $\mathcal{O}^{\dagger}\mathcal{O}$ operator.

[TABLE]

where $r_{1}=\sqrt{|\vec{r_{1}^{\prime}}|^{2}+t^{2}}$ and $r_{2}=\sqrt{|\vec{r_{2}^{\prime}}|^{2}+t^{2}}$ respectively. For the case where $\vec{r_{1}^{\prime}}$ and $\vec{r_{2}^{\prime}}$ are integrated over the whole lattice volume, the sum over the positions $\vec{x},\vec{y},\vec{x}+\vec{r^{\prime}_{1}}$ and $\vec{y}+\vec{r^{\prime}_{2}}$ can be carried out independently. Consequently, $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$ in the DI fluctuate independently which leads to a variance which is the product of their respective variances. In this case, the leading vacuum insertion is a constant, independent of $t$. This is the reason why the noise remains constant over $t$ in DI. On the other hand, the constant variance is reduced to $V_{R_{s}}/V$ when $r$ it is integrated to $R_{s}$, while the signal is not compromised. The sub-leading contribution (denoted by …) in Eq. (4) has an exponential decay in $t$ with a mass in the scalar channel. It is clear that to leading order, the ratio of the cutoff S/N at $R_{s}$ to that without cutoff is

[TABLE]

We shall consider several DI examples involving volume summations over two or more coordinates. Since the convoluted sum with a relative coordinate in Eq. (3) can be expensive, we shall invoke the standard convolution theorem by calculating the product of two functions $\tilde{K}(\vec{p},t)=\tilde{\mathcal{O}}_{1}(-\vec{p})\tilde{\mathcal{O}}_{2}(\vec{p}),$ where $\tilde{\mathcal{O}}_{1}(-\vec{p})/\tilde{\mathcal{O}}_{2}(\vec{p})$ is the Fourier transforms of $\mathcal{O}_{1}(\vec{x})/\mathcal{O}_{2}(\vec{x})$ in each configuration on their respective time slices. Then

[TABLE]

where $K(\vec{r^{\prime}},t)$ is the Fourier transform of $\tilde{K}(\vec{p},t)$. In this way, the cost of the double-summation, which is order $V^{2}$, is reduced to that of the fast Fourier transform (FFT) which is in the order of $V\,{\rm log}V$.

II.1 Scalar matrix element of the strange quark

The first example is the disconnected insertion for the nucleon matrix element with a scalar loop which involves a three-point function and can be expressed as

[TABLE]

where $S$ is the vacuum-subtracted scalar loop, $\mathcal{G}$ denotes the source grid for increasing statistics Li et al. (2010). Note here $r=\sqrt{|\vec{r^{\prime}}|^{2}+(t-\tau)^{2}}$ is the 4-D distance and $r_{x}$ is the spatial separation between the loop and the sink. Since the low-modes dominate the strangeness in the nucleon Yang et al. (2016b), we calculate the strange quark loop with low-modes only to illustrate the CDP effect. The sum over the spatial relative coordinate between the scalar quark loop $S$ and the sink interpolation operator $O_{N}(\vec{x},t)$ is carried out through the convolution in Eq. (6). This calculation is done on the domain-wall $32^{3}\times 64$ (32ID) lattice Blum et al. (2016b) with pion mass $\sim 170$ MeV and the lattice size is $4.6$ fm.

The upper panel of Fig. 2 gives the value and the error of the ratio of three-to-two point functions $C_{3}/C_{2}(R,\tau,t)$ in Eq. (7) as a function of $R$ at $\tau=4$ and $t=9$. The nucleon source-sink separation is $1.29$ fm in this case. We see the error grows after $R$ is greater than $\sim 12\,(\sim 1.7{\rm fm})$ while the central value remains constant within errors. This behavior reflects the fact that the signal without summing over $|\vec{r^{\prime}}|$ falls off exponentially with $r$, while the error remains constant as shown in the lower panel.

Fig. 3 shows $DI(scalar)$– the disconnected three-point to two-point function ratio to obtain the scalar matrix element for the strange quark in the nucleon as a function of $\tau-t/2$ for two source-sink separations at 1.00 fm (upper panel) and 1.57 fm (lower panel). Two results with cutoffs of $R_{s}=27$ and $R_{s}=12$ are plotted. $R_{s}$ is the cutoff radius for the relative coordinate between the sink and the quark loop in the spatial sum. It can be seen that the central values of the two cutoffs are all consistent within errors, while the errors with cutoff $R_{s}=12$ are smaller than the ones with cutoff $27$, which includes the whole spatial volume, by a factor of 4 or so. Thus, cutting off the spatial sum at the saturation distance is equivalent to gaining $\sim 16$ times more statistics.

II.2 Glueball mass

Next, we consider the glueball correlators in Eq. (3) on the 48I lattice with $La=5.5$ fm. The correlators from the scalar $E^{2}$ and $B^{2}$ operators and the pseudoscalar $E\cdot B$ operator are presented in Fig. 4, where they are plotted as a function of $R$ in Eq. (3) at $t=4$. We note the scalar correlators saturate after $R\sim 9\,(\sim 1.0{\rm fm})$ and the pseudo scalar one saturates after $R\sim 12\,(\sim 1.2{\rm fm})$, which can be understood in terms of the different ground state masses in these two channels. Again, comparing the error at $R=9$ to that at $R=24$, the latter includes the whole spatial volume, for the scalar case, the error is reduced by a factor of $\sim 4$ which is in reasonable agreement with the prediction of $\sim(\frac{24}{9})^{\frac{3}{2}}=4.4$ from Eq. (5). For the pseudoscalar case, the improvement is around 3 times and is consistent with the estimate of $\sim(\frac{24}{12})^{\frac{3}{2}}=2.8$.

II.3 Neutron electric dipole moment:

Finally, we examine the CP-violation phase $\alpha^{1}$ on the same 48I lattice which is needed for calculating the neutron electric dipole moment (nEDM). The phase is defined as

[TABLE]

for large enough $t$, where $C_{2}(t)$ is the common nucleon two-point function, $\Gamma_{e}=\frac{1+\gamma_{4}}{2}$ is the parity projector, $C_{3\mathcal{Q}}(t)$ is the nucleon propagator weighted with the total topological charge $\mathcal{Q}$

[TABLE]

We can turn the total topological charge into the summation of its density, i.e. $\mathcal{Q}=\sum_{x}q(x)$ where we use the plaquette definition for $q(x)$. Then the expression of $C_{3\mathcal{Q}}(t)$ with a cutoff $R$ can be cast in the same form as in Eq. (7), except the scalar quark quark loop $S$ is replaced with the local topological charge $q(x)$ and the sum of the topological charge density is over the four sphere with a radius $R$.

The result of $\alpha^{1}$ as a function of $R$ in Fig. 5 shows that the signal saturates after $R\sim 16$. Cutting off the sum of $q(x)$ at this R leads to a factor of $~{}\sim 3.6$ times reduction in error compared to the case of reweighting with the total topological charge as in Eq. (9). This example indicates that for four-dimensional sums, our new method employing the CDP can also improve the S/N. As we illustrated in the introduction, the nEDM from the $\theta$ term suffers from a $\sqrt{V}$ problem. It is shown here it is related to the vacuum insertion in the variance. This problem is resolved by turning the topological charge into a 4-D sum of the local charge density and applying the CDP by cutting off the relative 4-D distance in the sum.

III Systematic and statistical errors

So far, we have taken a simple cutoff $R_{s}=8/M$ to illustrate the efficacy of the variance reduction. This ad hoc choice inevitably incurs a systematic error. Since the asymptotic behavior of the integral of the correlator as a function of the separation $R$ is similar to that of the effective mass, we fit it as such and apply the Akaike information criterion (AIC) Akaike (1974) to obtain the statistical error and the estimated systematic error by using analysis with different fitting windows and models. The details of the application of AIC to the strange matrix element is provided in the Appendix as an example. It turns out that both the systematic and statistical errors are stable against multiple choices of windows and two fitting models. The statistical error is close to that at the cutoff distance when the plateau emerges (i.e. $R_{s}=8/M$). Using a representative fit with 2 fitting formulas and 80 combinations of 8 data points each for a total of 160 fits and 100 bootstrap samplings, we obtain the value of the strange matrix element that we considered earlier to be 0.160 (15) (15). The statistical error (first one) and the systematic error (second one) are practically the same. This is to be compared with the original value of 0.143(45) without taking the CDP into account. The systematic error is to be added to the total systematic error of the calculation.

IV Discussion and summary

Regarding the nucleon correlator in Fig. 1, we notice that there is no conspicuous increase of the error as a function of $R$. This is because, unlike the DI, the variance does not have an vacuum insertion for the CI. The leading contribution to the variance is expected to be $e^{-3m_{\pi}r}$. This has a longer range than that of the signal which falls off with the nucleon mass. Therefore, in principle, one would expect some gain in the S/N when $R$ is cut off at $R_{s}$. Therefore, the corresponding ratio of S/N in Eq. (5) is

[TABLE]

For the 48I lattice in Fig. 1, this ratio is 1.13 for the physical pion mass with the cutoff $R_{s}=8/M$. This is not nearly as much a gain as in the DI where the variance is dominated by the vacuum insertion. In the CI case, the noise saturates at $\sim 8/(3m_{\pi})=3.5$ fm. There is no gain for a lattice with a size larger than this.

In summary, we have shown that the exponential fall off of the Cluster Decomposition Principle (CDP) seems to hold numerically for the several correlators that we examined. For the disconnected insertions (DI), we find that the vacuum insertion dominates the variance so that the relevant operators fluctuate independently and is independent of the time separation. This explains why the signal fall off exponentially, while the error remains constant in the DI. To demonstrate the efficacy of employing the CDP to reduce the variance, we have restricted the volume sum of the relative coordinate between the operators to the saturation radius $R_{s}$ to show that there is an effective gain of $V/V_{s}$ in statistics without compromising the signal. This applies to all DI cases. For the cases we have considered, namely the glueball mass, the strangeness content in the nucleon, and the CP violation angle in the nucleon due to the $\theta$ term, we found that for lattices with a physical sizes of 4.5 - 5.5 fm, the errors of these quantities can be reduced a factor of 3 to 4. We have applied the Akaike information criterion (AIC) Akaike (1974) in the Appendix to estimate the systematic and statistical errors incurred by applying the CDP. For the strangeness content, we find that the systematic error is practically of the same size as that of the statistical one when the cluster decomposition is taken into account. This results in a 2 to 3 times reduction in the overall error. For the connected insertions, there is no vacuum insertion in the variance, the gain in statistics is limited.

V Acknowledgement

This work is supported in part by the U.S. DOE Grant No. $\text{DE-SC}0013065$. This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. This work used the Extreme Science and Engineering Discovery Environment (XSEDE) Stampede at TACC through allocation MCA01S007, which is supported by National Science Foundation grant number ACI-1053575 Towns et al. (2014). We also thank National Energy Research Scientific Computing Center (NERSC) for providing HPC resources that have contributed to the research results reported within this paper. We acknowledge the facilities of the USQCD Collaboration used for this research in part, which are funded by the Office of Science of the U.S. Department of Energy.

VI Appendix

VI.1 Error estimate using Akaike information criterion (AIC)

AIC is founded in information theory. It is an estimator of the relative relevance of the statistical models that are used to describe a given set of data Akaike (1974); AIC .

The usual task is to fit the data generated by some unknown process (function) $f$ with different trial models. However, we cannot tell which model is a better representation of $f$ with certainty in practice, because we do not know $f$. Akaike (1974) showed that we can estimate, via AIC, how much more (or less) information is lost when comparing one model to another. The estimate, though, is only valid asymptotically. If the number of data points is small, then some correction is often necessary (e.g., AICc) AIC .

VI.2 Definition of the AIC value

Considering a model $M$ for some data $x$ with $k$ parameters, the maximum value of the likelihood function for the model $L_{\max}$ is represented by the following probability

[TABLE]

where $\vec{\theta}$ is the parameter vector that maximizes the likelihood function. The maximum likelihood function is related to the minimum $\chi^{2}$ values in the standard fitting via

[TABLE]

The AIC value of the model is defined as (see ref. AIC and the reference therein for details)

[TABLE]

This favors models with the minimum AIC value and penalizes those with many fitting parameters.

VI.3 Practical application

In practice, we take a weighted average of all the models as is carried out in Ref. Borsanyi et al. (2015); Durr et al. (2016). The normalized weight for each model is

[TABLE]

where $i$ is the index of the models.

When handling the systematic errors of lattice calculations, we usually need to take consideration of various fitting formulas with different combinations of data points used for the fit. The AIC method can be helpful in these cases. Let’s consider a case where we have $N$ configurations and, for each configuration, there are $M$ data points (e.g. different time separation $t$ of a hadron two-point correlator). We plan to use $P$ models to fit the data (the number $P$ includes different fit ranges and different combinations of data points) in order to obtain the mean value, the systematic error and the statistical error of some model parameters (e.g. the mass of the ground state). The detailed procedure are given as follows:

The mean value is the weighted average of all the $P$ fit models (formulas and combinations of data points).

[TABLE]

$x_{i}$ is the fitting result from each model and $w_{i}$ is the normalized weight in Eq. (14). 2. 2.

The systematic error is taken to be the standard error of the weighted mean,

[TABLE] 3. 3.

The statistical error can be obtained from bootstrap resampling. We first do $N_{b}$ times bootstrap operation, in each bootstrap sample, we fit these $P$ models and save the weighted mean value. After that, we will have $N_{b}$ AIC weighted mean values. The bootstrap error of these weighted mean values gives the final statistical error.

In summary, we need to do $(1+N_{b})\times P$ times of correlated fits to obtain all the relevant results.

VI.4 The cluster decomposition case

In the cluster decomposition case, we shall consider the strange matrix element as a function of the cutoff radius $R$ as illustrated in the following figure. We can use the AIC method to estimate the mean value, the systematic error and the statistical error of the ratio.

To apply the AIC method to this particular case, we need first to determine our fitting formulas. In view of the fact that the ratio between the three-point function and the two-point function falls off exponentially as a function of the relative separation between the quark loop and the sink of the nucleon propagator, the accumulated sum of the separation with a cutoff $R$ is expected to be a constant after certain $R$, such as $R_{s}=8/M$, as illustrated in the upper panel of Fig. 2. This is much like fitting the effective mass plot to isolate the ground state. The two formulas (models) we will use are

[TABLE]

and

[TABLE]

The first one is the asymptotic form for $R\rightarrow\infty$ and the second one is the form commensurate with that from the cluster decomposition principle to cover more range of $R$ in the fitting.

Then, we need to choose the combinations of data points. To make sure that every combination has the same weight, we set the number of data points (marked as $N_{d}$) contained in each combination to be equal. And to enlarge the number of different combinations (marked as $N_{c}$) we can have, we do not force the points in each combination to be contiguous. (e.g., if the total range of data points is $R\in[8,27]$, and we set $N_{d}=4$, $N_{c}=5$, the possible combinations can be $[8,9,12,15]$, $[15,18,21,22]$, $[15,16,17,18]$, $[10,12,20,22]$ and $[20,22,24,27]$.) If $N_{c}$ is large enough, combinations will include both contiguous data points and noncontiguous data point. So in this sense, this is a more general way to choose data points.

Having the formulas and combinations, we can then proceed to do the fittings. The final results might be affected by 3 factors: $N_{b}$ (number of bootstraps), $N_{d}$ (number of data points) and $N_{c}$ (number of combinations). We vary these 3 factors to check if the results are stable in these fits.

The data range used is $R\in[10,27]$. After several hundreds of thousands of correlated fits, the results are show in the above table. The values and errors are pretty stable no matter how we vary $N_{b}$, $N_{d}$ or $N_{c}$. We decide to take the representative results from $N_{b}=100,N_{d}=8$, and $N_{c}=80$ which gives the value of $0.160(15)(15)$ as our final estimation. The systematic error (the second parenthesis) is comparable to the statistical one (the first parenthesis) from 100 bootstrap samples. This is to be compared to the original value of 0.143(45) when the sum over the relative coordinate between the quark loop and the sink of the nucleon propagator is carried out to cover the whole spatial volume. The total number of analysis $P$ is the product of the number of models (2 in this case) and $N_{c}$, the number of different combinations of $N_{d}$ data points. In this case, $P=2\times 80=160$.