-
Lattice QCD is a gauge invariant and nonperturbative regularization scheme for QCD that was introduced by K. Wilson in 1974 [1]. Many quantities can be obtained from the first principles by using a finite space-time lattice to simulate the interactions between quarks and gluons. The calculation of diagrams with disconnected quark loops is one of the most challenging problems in lattice QCD. Disconnected loops are needed in calculations of many lattice QCD observables, such as the nucleon electromagnetic form factors [2, 3], pi-NN coupling and pion polarizability [4], the mass of pseudoscalar flavor-singlet mesons [5, 6], strangeness and charm content of the nucleon [7, 8], hadronic scattering lengths and structure functions [9], and electron or muon hadronic
g−2 loop contributions [10]. In the lattice QCD, the expectation value of disconnected quark loops can be written as [11]⟨¯ψΓψ⟩=−Tr(ΓM−1),
(1) where
Γ∈{1,γμ,γ5,γ5γμ,σμν,μ,ν=1,2,3,4} , and M is the Dirac fermion operator. In this work, we use Wilson's Dirac operator, which can be written asM=1−κD,κ=12(am+4),
(2) D(n|m)αβab=4∑μ=1(1−γμ)αβUμ(n)abδn+ˆμ,m+4∑μ=1(1+γμ)αβU−μ(n)abδn−ˆμ,m,
(3) where a is the lattice spacing, and D is referred to as the hopping matrix. The real number
κ is the hopping parameter [12]. To calculate the disconnected diagrams, we need to solve the equationMx=b,
(4) where M is the Dirac matrix with dimension
K×K , and b is the source vector of dimensionK×1 . The solution of this equation isx=M−1b.
(5) The evaluation of disconnected quark loops requires
M−1 connecting arbitrary pairs of lattice points. Generally, Wilson's Dirac matrix is a large sparse matrix and its typical dimension K is from106 up to109 . Its direct evaluation is prohibitively expensive, both in terms of computer time and memory.In order to calculate disconnected quark loops, specific techniques must be introduced. Unbiased noise methods are traditionally used to estimate the inverse matrix [13, 14]. The truncated solver method (TSM) [15] and the spin explicit method (SEM) [16, 17] were introduced to reduce the stochastic error. A probing technique is also a way to deal with this problem [18]. However, we found that all stochastic methods result in large errors in the calculations of disconnected quark loops
Tr(ΓM−1) . Hence, we introduce the symmetric multi-probing source (SMP) method to calculate all disconnected quark loops, and compare the results with theZ(2) noise and spin-color explicit (SCE) methods. All comparisons are based on the point source results which are taken as exact. -
In general, the inverse of a large sparse matrix can be calculated by using the unbiased stochastic method. Here, we briefly review the
Z(2) noise method. We assume L column noise vectorsb1 ,b2 ,b3 ,… ,bL , which have the following two properties [13],⟨bi⟩=1LL∑l=1bli=O(1/√L),
(6) ⟨bibj⟩=1LL∑l=1bliblj=δij+O(1/√L),
(7) where
bli is thei -th entry in the noise vector l. The stochastic average<⋅⋅⋅> is taken over the ensemble of noise vectors L.IfZ(2) noise vectors are used in Eq. (4), we obtainxli=∑kM−1ikblk.
(8) The inverse matrix element,
M−1ij , is given as<bjxi>=∑kM−1ik<bjbk>=M−1ij.
(9) The variance of the
Z(2) noise method is [19]σ2≡1LK∑i≠j|M−1ij|2.
(10) It can be seen that the stochastic error of the
Z(2) noise estimate results only from the off-diagonal elements of the inverse matrix.Therefore, we obtain
Tr(ΓM−1) Tr(ΓM−1)=∑jΓM−1jj=∑i,jΓM−1ijδij=∑i,j1L∑lΓM−1ijbliblj=∑i1L∑lbli(Γxli).
(11) In order to reduce the error of the
Z(2) noise method, we applied an independent stochastic inversion of the spin and color components (spin-color explicit method, SCE method) for eachZ(2) noise vector, similar to SEM [16]. -
The SMP source vector
ϕP is introduced as follows [20]ϕP(S(x,P),α,a)=∑y∈S(x,P)ψ(y,α,a),
(12) where
α is the Dirac index and a the color index.S(x,P) represents the sites with the same color of x obtained by the symmetric coloring schemeP(nsd,nsd,nsd,ntd,mode) , wherens andnt are the spatial and temporal sizes of the lattice. d is the distance parameter, andmode=0,1,2 corresponds to the Normal, Split and Combined mode. x is the seed site at(x1,x2,x3,x4) and y are the other lattice sites belonging to the setS(x,P) .ψ is the normalized point source vector. The number of SMP sourcesNSMP that cover all lattice sites isNSMP={12d4mode=0,24d4mode=1,6d4mode=2.
(13) Applying the SMP source in Eq. (4) to calculate the trace of
ΓM−1 , we obtainTr(ΓM−1)SMP=∑x,α,aψ(x,α,a)(ΓM−1)ϕP(S(x,P),α,a)=∑x,α,aψ(x,α,a)(ΓM−1)ψ(x,α,a)+y≠x∑y∈S(x,P)∑α,aψ(x,α,a)(ΓM−1)ψ(y,α,a)=Tr(ΓM−1)+y≠x∑y∈S(x,P)∑α,aψ(x,α,a)(ΓM−1)ψ(y,α,a),
(14) where the second term in the last line is the sum of off-diagonal elements of
ΓM−1 . Considering the space-time locality of M, this term can be regarded as the error in the calculation ofTr(ΓM−1) . If we choose the proper scheme P of the SMP source, the error is quite small. In this case, we can neglect the error term and getTr(ΓM−1)≈Tr(ΓM−1)SMP.
(15) -
It is very time consuming to solve Eq. (4) if b is the point source (p-s) vector which runs over all lattice sites of a large lattice. In this work, we use the point source method to evaluate
Tr(ΓM−1) , and take the result as exact, and use it for comparison. We only work with ensembles of Iwasaki pure gauge configurations, where the volume isL3σ×Lτ= 123×24 . The lattice spacing of the ensembles isa≈0.1 fm. The analysis is performed on25 configurations withκ=0.151 , corresponding tomπ≈488 MeV. We show the results forTr(M−1) andTr(γ5M−1) in the main body of this paper. The results forTr(γ3M−1) ,Tr(γ5γ1M−1) andTr(σ34M−1) , selected as representative, are presented in appendix A.Using the identity
M=γ5M†γ5 [21], we obtainTr(γ5M−1)=Tr(γ5γ5(M−1)†γ5)=Tr((M−1)†γ5)=Tr((γ5M−1)†).
(16) Therefore,
Tr(γ5M−1) is a real number, andTr(M−1) is also real. The absolute errorΔr discussed in this work is defined asΔr=|rp−rm|,
(17) where
rp is the result for the point source andrm is the result for any other method.We also define the method error
σ asσ2=1Nconf(Nconf−1)Nconf∑i=1Δr2i.
(18) Δri is the absolute error of the i-th configuration (conf), andNconf is the number of configurations. In the definition ofσ , we userp and not the average value⟨r⟩c of all configurations, to eliminate the fluctuations due to different configurations. Therefore, this definition gives only the error of the method. -
The results and absolute errors of the first ten configurations are shown as representative in this paper. However, we show the method error for all 25 configurations. For the SMP method, the results for
Tr(M−1) are shown in Table 1, and in Table 2 forTr(γ5M−1) . The number of source vectors for all methods is the same.conf. p-s N=15552(6,0) N=3072(4,0) N=1536(4,2) N=384(2,1) N=192(2,0) N=96(2,2) 1 453821.05 453821.93 453834.19 453831.06 453821.88 453785.60 453890.51 2 454042.56 454049.54 454036.58 454040.99 454020.87 453919.91 453945.70 3 454216.18 454218.24 454219.28 454218.60 454223.99 454389.67 454395.52 4 454002.41 453999.82 453998.76 453994.74 453964.91 454033.19 454069.26 5 454354.48 454357.04 454337.75 454319.32 454291.41 454229.23 454220.83 6 454301.64 454299.51 454308.33 454298.46 454274.90 454210.82 454231.96 7 454085.53 454087.16 454089.11 454085.17 454083.09 454065.35 454038.91 8 454059.63 454060.60 454054.26 454053.96 454016.87 454104.69 454054.37 9 453944.56 453941.69 453936.50 453930.76 453919.66 453778.34 453765.85 10 453914.30 453915.53 453924.14 453923.83 453936.33 453826.90 453825.25 Table 1. Results for
Tr(M−1) for the SMP method. N is the number of SMP source vectors. The numbers in the parentheses represent the parameters d andmode of the SMP source; for example,(6,0) denotesd=6 andmode=0 . The SMP method is a good approximation when the number of source vectors is large.Figure 1. (color online) Absolute error of
Tr(M−1) for the SMP method. The numbers in the parentheses represent the parameters d andmode of the SMP source vectors; for example,(6,0) denotesd=6 andmode=0 .conf. p-s N=15552(6,0) N=3072(4,0) N=1536(4,2) N=384(2,1) N=192(2,0) N=96(2,2) 1 −18.94 −18.37 −22.18 −15.61 −10.08 −21.66 −22.73 2 202.86 203.21 202.90 212.26 192.83 184.64 188.23 3 183.32 184.33 181.43 180.09 178.41 178.20 148.52 4 −170.65 −170.89 −163.42 −160.16 −186.75 −172.39 −206.26 5 −586.77 −583.80 −563.81 −555.51 −559.26 −569.70 −519.76 6 74.85 73.18 73.35 75.99 81.00 96.71 124.26 7 −288.10 −290.81 −288.23 −284.46 −291.66 −280.63 −271.28 8 141.12 141.86 137.16 138.44 153.09 167.57 164.22 9 −338.58 −335.66 −342.16 −329.44 −338.10 −332.02 −351.82 10 −223.22 −220.96 −224.15 −221.93 −237.32 −218.63 −222.10 Table 2. Results for
Tr(γ5M−1) for the SMP method. N is the number of SMP source vectors. The numbers in the parentheses represent the parameters d andmode of the SMP source.Figure 2. (color online) Absolute error of
Tr(M−1) for theZ(2) noise method.NZ is the number ofZ(2) source vectors.Figure 3. (color online) Absolute error of
Tr(M−1) for the SCE method.NS is the number of SCE source vectors.Figure 4.
(color online) σ ofTr(M−1) as a function of the number of source vectors N for the three different methods.In Figs. 1, 2 and 3, the absolute errors of
Tr(M−1) for the SMP,Z(2) noise and SCE methods are shown. As expected, as the number of source vectors increases, the absolute errors become smaller for all three methods.In Fig. 4,
⟨Δr⟩ ofTr(M−1) for the SMP and SCE methods are compared with theZ(2) noise method for the same number of source vectors, indicating that when the number of source vectors is small, the absolute errors of SMP and SCE are larger than ofZ(2) . The figure also shows that the SCE method does not bring any improvement in the evaluation of the scalar disconnected quark loops compared with theZ(2) noise method. Hence, theZ(2) noise method is an efficient method for the calculation ofTr(M−1) .As the number of source vectors increases, the results of the SMP, SCE and
Z(2) noise methods become more and more accurate. The results forTr(M−1) show that theZ(2) noise method is a good choice for evaluating the scalar disconnected loops.The results for
Tr(γ5M−1) for the SMP method are shown in Table 2. The absolute errors ofTr(γ5M−1) are presented in Figs. 5, 6 and 7. As the number of source vectors increases, the absolute errors of all methods decrease, as expected. Fig. 6 shows that theZ(2) noise method gives a too large error ofTr(γ5M−1) , especially for small N.In Fig. 8, the method error of
Tr(γ5M−1) is shown, indicating thatσ decreases with increasing number of source vectors for all three methods. The SCE method results in a smallerσ than theZ(2) noise method for the same number of source vectors. Hence, the SCE method is an obvious improvement compared with theZ(2) noise method. Also, the SMP method gives a much smaller method error than theZ(2) noise method, especially for a small number of source vectors.All these results indicate that the
Z(2) noise method gives a larger error than the SMP and SCE methods for the same number of source vectors in the calculation ofTr(γ5M−1) . Compared with theZ(2) method, the SCE method indeed improves the calculation ofTr(γ5M−1) . On the other hand, the absolute errors of the SMP method are smaller than the SCE method. The results of the SMP method are quite precise when the parameter d is large enough. Hence, the SMP method has a considerable advantage over the SCE andZ(2) noise methods in the estimation ofTr(γ5M−1) .In Table 3,
σ/σZ(2) of the SMP and SCE methods for different N are presented. In the caseΓ=1 , theZ(2) method is the best of the three methods, while forΓ=γ5 , the SMP method gives better results than the others.σ for the SMP method is 22.15% ~ 38.31% of the Z(2) noise method, and about half of the SCE method. Therefore, the SMP method is a considerable improvement in the estimation of pseudoscalar disconnected loops.Γ method N=96 N=192 N=384 N=1536 N=3072 N=15552 1 SMP 199.22% 371.77% 179.12% 141.37% 138.09% 105.91% SCE 99.56% 148.38% 176.85% 161.88% 164.25% 143.99% γ5 SMP 35.77% 38.31% 25.13% 32.53% 23.36% 22.15% SCE 65.87% 81.97% 65.56% 56.26% 51.62% 52.64% Table 3.
σ/σZ(2) for the SMP and SCE methods for different N. -
We have calculated in this work all disconnected quark loops using the SMP, SCE and
Z(2) noise methods, and presented in the paper a detailed analysis of the results for the scalar and pseudoscalar disconnected quark loops. As expected, the absolute errors of the SMP, SCE andZ(2) noise methods decrease as the number of source vectors in the evaluation of the disconnected quark loops increases. In the case of scalar disconnected loops, it was shown that the SCE method does not have an advantage over theZ(2) noise method. The absolute errors of theZ(2) noise method are smaller than of the SMP and SCE methods even if the number of source vectors in the calculation of scalar disconnected diagrams is small.We also found that the SMP method can improve the precision of the calculation of
Tr(γ5M−1) by about a factor of 2.5-4.6 compared to theZ(2) noise method, and that the SCE method does not have an advantage over the SMP method. The results show that of the three methods, SMP is the best for the calculation ofTr(γ5M−1) . We believe that the SMP method is suitable due to the operator hermiticity, but the reason why it gives a significant improvement in the calculation of pseudoscalar disconnected quark loops requires further study.Numerical simulations have been performed on the Tianhe-2 supercomputer at the National Supercomputer Centre in Guangzhou (NSCC-GZ), China.
-
In this appendix, we show
<Δr> in Figs. A1, A2, A3 andσ/σZ(2) in Table A1 ofTr(γ3M−1) ,Tr(γ5γ1M−1) andTr(σ34M−1) , which indicate that the SMP method does not have a significant advantage compared with theZ(2) noise and SCE methods. As the number of source vectors increases, the absolute errors of all methods decrease. When the number of source vectors is large enough, all methods give quite similar results.Figure A2. (color online) Method error of
Tr(γ5γ1M−1) as a function of the number of source vectors N.Figure A3. (color online) Method error of
Tr(σ34M−1) as a function of the number of source vectors N.Γ method N=96 N=192 N=384 N=1536 N=3072 N=15552 γ3 SMP 205.18% 269.00% 174.39% 129.72% 123.95% 101.65% SCE 134.94% 176.29% 153.68% 174.92% 142.95% 154.08% γ5γ1 SMP 66.81% 63.66% 55.92% 37.60% 52.31% 47.30% SCE 43.32% 52.49% 46.56% 41.89% 57.33% 48.15% σ34 SMP 72.78% 54.20% 47.71% 64.61% 70.62% 73.00% SCE 84.09% 60.01% 79.27% 87.96% 90.24% 133.93% Table A1.
σ/σZ(2) for the SMP and SCE methods for different N.
conf. | p-s | ![]() | ![]() | ![]() | ![]() | ![]() | ![]() |
1 | 453821.05 | 453821.93 | 453834.19 | 453831.06 | 453821.88 | 453785.60 | 453890.51 |
2 | 454042.56 | 454049.54 | 454036.58 | 454040.99 | 454020.87 | 453919.91 | 453945.70 |
3 | 454216.18 | 454218.24 | 454219.28 | 454218.60 | 454223.99 | 454389.67 | 454395.52 |
4 | 454002.41 | 453999.82 | 453998.76 | 453994.74 | 453964.91 | 454033.19 | 454069.26 |
5 | 454354.48 | 454357.04 | 454337.75 | 454319.32 | 454291.41 | 454229.23 | 454220.83 |
6 | 454301.64 | 454299.51 | 454308.33 | 454298.46 | 454274.90 | 454210.82 | 454231.96 |
7 | 454085.53 | 454087.16 | 454089.11 | 454085.17 | 454083.09 | 454065.35 | 454038.91 |
8 | 454059.63 | 454060.60 | 454054.26 | 454053.96 | 454016.87 | 454104.69 | 454054.37 |
9 | 453944.56 | 453941.69 | 453936.50 | 453930.76 | 453919.66 | 453778.34 | 453765.85 |
10 | 453914.30 | 453915.53 | 453924.14 | 453923.83 | 453936.33 | 453826.90 | 453825.25 |