Investigating the topological structure of quenched lattice QCD with overlap fermions using a multi-probing approximation

The topological charge density and topological susceptibility are determined by a multi-probing approximation using overlap fermions in quenched SU(3) gauge theory. Then we investigate the topological structure of the quenched QCD vacuum, and compare it with results from the all-scale topological density. The results are consistent. Random permuted topological charge density is used to check whether these structures represent underlying ordered properties. The pseudoscalar glueball mass is extracted from the two-point correlation function of the topological charge density. We study 3 ensembles of different lattice spacing a with the same lattice volume 163×32. The results are compatible with the results of all-scale topological charge density, and the topological structures revealed by multi-probing are much closer to all-scale topological charge density than those from eigenmode expansion.


Introduction
The topological susceptibility χ of Yang-Mills theory is an important physical quantity in explaining the U (1) A anomaly and the mass of η , e.g. the study of this problem in the N c → ∞ limit brings up the famous Witten-Veneziano formula for the η mass [1][2][3]. In the continuum limit, the topological susceptibility as a function of four-momentum k, denoted by χ(k 2 ) in the Euclidean version, is defined as where the quantity C qq (x)= 0|T (q(x)q(0))|0 is the twopoint correlator of the topological charge density q(x) and |0 denotes the QCD vacuum. The contribution of the U (1) A anomaly to the η mass is given by [1,2,4] in which Y M indicates pure Yang-Mills theory, N f is the number of flavors, and F π is the corresponding π decay constant. This leads to χ Y M (0)≈ (180 MeV) 4 [2,4].
To calculate q(x) and χ(0) of SU (3) pure Yang-Mills theory on the lattice, we have to define the topological density q(x) first. There are several ways to define it. One usual way is to define q(x) with the field strength tensor F µν on the lattice as [5] q L (x)=− 1 32π 2 µνρσ µνρσ T rF µν (x)F ρσ (x), (3) in which subscript L denotes the lattice. This definition has the classical continuum limit q L (x) a→0 −−→ q(x). On the lattice, however, the direct calculated topological charge Q L = x a 4 q L (x) usually deviates from integral values, and Q L = Z(β)Q+O(a 2 ) [6][7][8], with the multiplicative renormalization coefficient Z(β) dependent on the inverse coupling constant β. This dependence can be computed perturbatively, or some kind of smoothing steps can be used to make the effective inverse coupling constant β eff larger and such that Z(β eff ) → 1 [8] is achieved. for χ(0), it has been found that multiplicative and additive renormalization terms both exist as χ L (0) = a 4 Z 2 (β)χ(0) + a 4 A(β) < S(x) > n.p. +P (β) [6], where S(x) is the action density, n.p. indicates nonperturbative, and A(β) and P (β) are ultra-violet effects that can be approximated in perturbation theory. We could also smooth configurations to approach integer values of Q L for calculating χ(0), but the smoothing procedure would influence q(x).
Another way to define q(x) is using a massless Dirac operator D(0,x,y) that obeys the Ginsparg-Wilson relation Dγ 5 + γ 5 D = aDγ 5 D. Here the first argument in D(m,x,y) denotes the mass of the quark while the second and third index are the space-time index on the lattice. For simplicity, we will also use D(0) to denote the massless Dirac operator with the space-time indices suppressed. With this notation, the topological charge density can be defined as q L (x)= 1 2 tr cd γ 5 D(0,x,x), where c and d are the index of color and Dirac spinor respectively. This is a proper approach that guarantees the topological charge Q L an integral value since Ginsparg-Wilson fermions satisfy the Atiyah-Singer index theorem [9] in which n L and n R are the number of left-and righthanded zero modes of D(0). A massless overlap fermion is exactly a Ginsparg-Wilson fermion [10,11], so in this paper we will use the massless overlap Dirac operator D ov (0) to calculate the topological charge density. The definition of D ov (0) reads, where D w (s) is the Wilson-Dirac operator, with parameter |s| < 1 for the optimization of the locality of the overlap operator [12]. D w (s) = D w −(1+s), where D w is the massless Wilson-Dirac operator with mass parameter m = −(1+s). In this work we choose parameter κ= 1 2(am+4) =0.21, which means s=0.619.
However, it will cost too much computer resources to calculate topological charge density using the definition in Eq. (4) directly, i.e., using point sources to calculate the trace. The resulting topological charge density q(x) is named all-scale topological charge density in Ref. [13]. An approximation method called eigenmode expansion has been proposed in the literature [13][14][15]: where ψ λ is the eigenmode of D(0) with eigenvalue λ, and λ cut is the supremum of the absolute eigenvalue of the corresponding eigenmodes included in the eigenmode expansion.
Only if all zero modes are included, which obey the relationships γ 5 ψ λ = ±ψ λ , D ov (0)ψ λ = 0, would q λcut (x) sum the individual local chiralities, and the topological charge of eigenmode expansion Q λcut = Q L . However, some features such as the positive core of two-point correlator of topological charge density C qq (x) would increase, and the height of peaks of C qq (x) would decrease. Even the negative region of C qq (x) would disappear when λ cut is too small [13]. Therefore, more eigenmodes should be considered, which would also cost too much computer time, and it would be unclear what would be discarded in the eigenmodes of larger λ cut .
In this paper we will use a multi-probing method [16] to approximate q(x), owing to the finite range of D(0) [13,17]. ∆Q, the deviation of Q L from an integer, is regarded as a parameter that signals the validity of this approximation. The lattice version of Eq. (1) is defined as with the lattice version of two-point correlation function . Noting that the true value of expectation of topological charge is 0, the lattice version of topological susceptibility at k 2 = 0 can be expressed as which obviously has the same accuracy as topological charge. It is believed that the U (1) A problem [18], topological susceptibility, η mass, θ dependence and spontaneous chiral symmetry breaking are all related to the vacuum fluctuation of topological charge density [1,19]. So, in Section 3 we attempt to investigate the local topological structures of the quenched QCD vacuum, hoping that this will reflect some ordered features of the quenched QCD vacuum. Since they can be explained by the instanton model, which was first introduced by Polyakov [20,21], it is reasonable to try to describe the structure of topological charge density in terms of the instanton liquid model, as the topological charge density shows sign-coherent unit-quantized lumps locally. In the past, several papers [22][23][24][25][26][27][28][29][30] have attempted to justify this idea, but these papers used methods such as cooling, smoothing and smearing, which will lead to configurations towards classical solutions, therefore the results should be taken with some reservations. On the other hand, it has been predicted that QCD vacuum does not resemble the instanton liquid model [19], and this idea was further supported by lattice evidence in Ref. [14]. In fact, there is no four-dimensional global structure of topological charge density but only lower dimensional structures survive [13][14][15]. In Section 3 we will discuss more details about the structure of topological charge density, and make comparisons with the results of allscale q(x) in Ref. [13], which were built up from three ensembles only with 53, 5 and 2 configurations, due to the large amount of computer time for calculating allscale topological charge density.

Multi-probing approximation
The massless overlap operator D ov (0) is not ultralocal but local. It is exponentially decaying with respect to lattice distance with a decay rate that is almost independent of inverse coupling constant β [13,17] D ov (0,x,x+r)∝exp(−µr/a) , µ>0, (9) in which r is the lattice Euclidean distance, and a is the lattice spacing, so we have D ov (0,x,x+r) D ov (0,x,x) when r/a is large enough. Therefore we can approximate topological charge density q(x) with massless overlap fermions as in which {r i } are some distances that satisfy exp(−µr i /a) 1, and c,d are indices for Dirac and color respectively.
In actual computing, we choose a set of sites {x i } on the lattice as a group of probes called a multi-probe. These sites are distributed uniformly at intervals equal to ∆r. The interval in lattice units ∆r/a should be an integer number that is a common divisor of space-time lengths of the lattice, or the interval is √ 2∆r when the sites in the same group are separated by odd and even into two new groups [16]. In detail, a group {x i } with intervals equal to ∆r or √ 2∆r is chosen by steps as follows: 1) First, for groups with intervals equal to ∆r, we pick the source site x 0 for a group {x i }. The spacetime coordinates are x 0 ={x 1 0 ,x 2 0 ,x 3 0 ,x 4 0 }, in which the superscripts 1,2,3 denote spatial notations, and 4 denotes the temporal notation. Every space-time coordinate x α 0 should satisfy x α 0 ∈[0,∆r−a] due to the interval we have chosen, where the superscript α stands for space-time index, which means that we can find (∆r/a) 4 different source sites, since every space-time coordinate has ∆r/a choices, and every source site corresponds to a group {x i }.
2) Second, the coordinates 0 +n α ×∆r, in which n α are nonnegative integers. Then any site of this group can be labeled by When we divide a group of sites {x i } into two groups according to whether 4 α=1 n α is odd or even, then we have groups with intervals equal to √ 2∆r. The whole set of lattice sites is then divided into N mp groups. Every group includes N x N y N z N t /N mp sites, and N mp =(∆r/a) 4 for interval ∆r and only has a non-zero value on sites {x i } belonging to the corresponding group. Based on Eq. (9), the diagonal element (γ 5 D(0,x,x))cd cd , x∈{x i } can now be probed by in which δ(c,d,x) is point source. Thus only 12N mp overlap operator multiplications with vector sources have to be performed with the multi-probing method, compared with 12N x N y N z N t if calculating all-scale q(x) directly, saving a substantial amount of computation, especially for large lattices.
In Table 1 we present the four ensembles of SU (3) pure gauge configurations in our study, which were generated with periodic boundaries using Iwasaki action.
Three of them have the same lattice size of N s = 16, N t =32, and the other has a smaller lattice size of N s =12, N t = 12. We use Wilson flow with w 0 = 0.1670(10)fm to set the scale of lattice spacing a [31][32][33]. In the calculations of topological charge density, we choose the scheme of intervals equal to 4 √ 2a to approximate the all-scale topological density q(x) for the ensembles with lattice size 16 3 ×32, which ensures that almost all configurations have ∆Q < 0.1, in which ∆Q is the deviation of Q = x q mp (x) from an integer, since the massless overlap operator D ov (0) should satisfy the Atiyah-Singer index theorem. Therefore only 4 4 ×2×3×4 = 6144 times of matrix multiplication to vector are needed. For calculating all-scale q(x) directly, 16 3 ×32×3×4 = 1572864 times will be needed, which means only 6144/1572864= 1/256 the time is needed for multi-probing method compared with calculating all-scale q(x) directly. In columns 5 and 6 we list the percentage of the configurations that satisfy the condition ∆Q below values 0.08 and 0.05 respectively, and the obtained values for χ 1 4 (0) using the multi-probing method are tabulated in last column. These values are compatible with phenomenological estimates χ Y M ≈ (180MeV) 4 [2,4] and the result χ Y M = (191±5MeV) 4 [34] in the continuum limit by using Eqs. (5)- (8). The ensemble including 3 configurations with lattice size of 12 3 ×12 is calculated with intervals equal to 4a, and the results q mp (x) are only used to compare with the all-scale q(x) by Eq. (12) in Fig. 2(a).
In Fig. 1, we present topological charge density computed both by all-scale q(x) directly and by the multiprobing method with intervals equal to 4 √ 2a on a 2D surface. Figure 1(b) only differs slightly from Fig. 1(a). The features of q(x) such as peaks, valleys and ridges are consistent. More detailed analyses of the structure of q(x) are discussed in Section 3.2.  In Fig. 2, we use ∆(q cut ) to investigate the deviation of q mp (x) from q(x), in which ∆(q cut ) is defined as: The all-scale topological densities q(x) of the three configurations with lattice size of 12 3 ×12 are computed over all lattice sites. We only use 2000 random sites for each of the ten configurations with lattice size of 16 3 ×32 due to the huge amounts of calculations needed. These configurations are from the two ensembles of β =9.45 shown in Table 1. The q max in Fig. 2 is the maximum abso-lute topological density of q(x) in the lattice. We use the density that has absolute value larger than q cut since the vacuum is thought to be mostly driven by the intense fields [15], and we will find that the topological structures from intense field are more important in Sections 3.2 and 3.3. Both subgraphs in Fig. 2 show that the deviation drops sharply when q cut increases just above 0, and only a few percent deviation remains when q cut reaches 0.05q max , which means that the large deviation comes from the lowest densities, which are influenced by short range fluctuations and which we are not interested in. The variance σ 2 of ∆(q cut ) also decreases as q cut increases. In the region q cut >0.05q max it is small enough to conjecture with confidence that all high q mp (x) only have a few percent deviation from q(x) when we use intervals of 4 √ 2a, even though only 2000 random sites from every configuration are used in Fig. 2(b). We find that the behaviors mentioned above are consistent among config-urations of different topological charge Q. Note that the magnitude in the left-hand panel is 4 times larger than that on the right.
The scheme of multi-probing with intervals of 4 √ 2a is presented for a 2D surface in Fig. 3(a). Sites with same number belong to the same multi-probing group, and the 3D sketch in Fig. 3(b) shows that the diagonal points, which are on the surfaces of a cube with side length 4a in the lattice, belong to the same multiprobing group. As for the scheme of multi-probing with intervals of 4a, vertices of a cube with side length equals 4a should belong to the same multi-probing group. In other words, sites with number 1 and 2 in Fig. 3(b) belong to the same group now, so the number of groups will be reduced by half. Therefore the schemes of multiprobing with interval equal to 4 √ 2a and 4a are correlated strongly. The same relationship exists universally between schemes with intervals of n √ 2a and na.   3 Analysing the topological charge density q(x) Topological susceptibility can be computed from topological charge individually, and the topological charge Q should be an integer value, which can be used to check the degree of approximation. Although the multi-probing method seems to approximate Q rather well, as indicated in Table 1, the performance of topological charge density, and especially the slope of the susceptibility, is still a concern. Thus, we need to find some other parameters to test the validity of the approximation, e.g. we may compare the two-point correlator C qq (r) from multi-probing and all-scale q(x). Other tools proposed in Ref. [13] will be utilized to analyze the vacuum structure of topological charge density. The results will also be compared with those obtained from all-scale q(x).

C qq (r) from multi-probing method
The reflection positivity [35] requires that the twopoint correlator of topological charge density should be negative [4,36,37]: Together with the massless overlap operator D ov (0) having a finite range, C qq (r) with overlap fermions will have a positive core, which has been studied extensively already.
In Figs. 4(a) and 4(b) the two-point correlator from our work displays a positive core and negative tails, just as all-scale q(x) does in Figs. 4(c), 4(d) (the results of Fig. 14 in Ref. [13]). The details shown in Figs. 4(b) and 4(d) demonstrate that the ranges of positive cores are consistent with β =9.45 (a≈0.1 fm) in Fig. 4(b) and with β = 8.45,8.60(a ≈ 0.1 fm) in Fig. 4(d) with the same lattice size 16 3 ×32. From Fig. 4(b), ensembles with different lattice spacing a have different positive cores in physical units, but are almost the same in lattice units, which means it is independent of β in lattice units. This effect is because the finite range of overlap operator D ov (0) has little dependence on β. In other words, as lattice spacing a decreases, the positive core of C qq (x) will disappear in physical units. When a → 0, only the contact term C qq (0) > 0, and the reflection positivity condition Eq. (13) will be restored in the continuum limit. (c) Fig. 14(a) in Ref. [13] (d) Fig. 14(b) in Ref. [13]  Inspired by the instanton model, a functional form of the two-point correlator C qq (x) in the negative region can be used to extract the pseudoscalar particle mass [38]: where K 1 (z) is the modified Bessel function, which has the asymptotic form We can use Eqs. (14) (15) to fit the data of C qq (r). In the fitting procedure the amplitude and mass are treated as free parameters and then we can extract the mass of the pseudoscalar glueball in quenched QCD, since only glueballs exist in pure gauge theory. The fitting range is determined by the effective mass m eff (r) which is defined by C qq (r+∆r) C qq (r) = r r+∆r where we set ∆r = 0.5a and C qq (r) is averaged over [r− ∆r 2 ,r+ ∆r 2 ]. In Fig. 5 we present the effective mass plateaus of 2 ensembles from Table 1 with errors computed by the jackknife method, except for ensemble a = 0.1161 fm, whose lattice spacing may be too coarse. The starting and ending points of the fitting range vary in [3.25a,5.25a] for the a = 0.0769 fm ensemble and in [3.5a,4.5a] for the a = 0.0963 fm ensemble. They are the plateaus shown in Fig. 5. The fitting range is set to be greater than or equal to ∆r = 0.5a, then the final fitting range is decided by the χ 2 /d.o.f. that is nearest to 1.0 [39].   Table 2 presents the fitted results for the masses for ensembles a=0.0963, 0.0769 fm, both of which are compatible with the value 2560(35) MeV from Ref. [40], though the fitting range is determined somewhat subjectively. The mass has little dependence on the lattice spacing a.

Clusters of topological charge density from multi-probing method
We use the cluster analysis proposed in Ref. [13] to investigate the sign-coherent structure of topological charge density (the word "sign-coherent" will be omitted below for simplicity): 1. The cluster on the lattice is composed of sites that are link-connected of same sign of the topological charge density q(x). q cut means that only sites with |q| > q cut are taken into account, and q max is the maximum value of |q(x)|.

The connectivity is defined by
where r max = L x 2 +L y 2 +L z 2 +L t 2 /2 is the largest distance available on the periodic lattice, and Θ c (x) is a characteristic function of the cluster c, namely Θ c (x) = 1 , if x ∈ c ,else Θ c (x) = 0. The value q cut =q perc for which f (r max )=0 is called the onset of percolation. When q cut <q perc ,f (r max )>0, indicating that global clusters exist.
3. The Euclidean distance between the two largest clusters c and c is defined by : 4. A random walk method is utilized to define the fractal dimension d * of a cluster. Starting from the site that has the maximum q max inside a cluster, and after τ steps of independent random walks, the probability of returning to the starting point P (τ )=(2πτ ) −d * /2 . The walk jumps to a linked site in the same cluster with equal probability at every step.
Since the physics of the vacuum is thought to be mostly affected by intense fields, and the super-longdistance might be relevant to the physics of the QCD vacuum such as Goldstone boson propagation [15], we are concerned about the large topological charge density (q cut is large). At the same time global structures still survive, i.e. f (r max )>0.
In Fig. 6, the number of clusters, which we call cluster multiplicity, increases as q cut increases till it reaches the maximum, then decreases. The fractional size of the largest cluster, which denote the size of the largest cluster relative to total volume occupied by all clusters, reaches a minimum near the maximum point of cluster multiplicity, then rises. The trend of fractional size of the largest cluster is exactly opposite that of cluster multiplicity. At q cut =0, the two largest clusters with opposite sign both have fractional size of the largest cluster ≈ 0.5. Except the two largest clusters, the dozens of clusters left are only fragments with volume of 1 or 2, and their total volume compared with the lattice volume 16 3 ×32 is so small that it can be ignored.
In Fig. 7, the distance d(c,c )≈ √ 3a between the two largest clusters at q cut = 0 is independent of β, like the results shown in Fig. 7(b), and just as for the positive core of two-point correlator of topological charge density. This means that at q cut = 0, the two largest clusters of opposite sign tangle with each other closely. When a→0, they are just like two tangled 3D papers dominating the whole 4D space. After the onset of percolation, the distance between the two largest clusters reaches a plateau, which is also independent of β. This has been proved in Fig. 7(b) from Ref. [13]. The turning point of cluster multiplicity, fractional size of the largest cluster and the plateau of the distance between two largest clusters all happen near the onset of percolation at q cut ≈0.2q max . Fig. 6. (color online) The left-hand plots are from the multi-probing method, and the right-hand plots from all-scale q(x) in Ref. [13]. The first row shows the total number of clusters, and the second row the size of the largest cluster relative to all clusters. Fig. 7. (color online) The left-hand plots are from the multi-probing method, and the right-hand plots from allscale q(x) in Ref. [13]. The first row shows the distance between the two largest clusters, and the second row the connectivity f (rmax) of the largest cluster.    Fig. 9. (color online) The average fractional size of largest clusters relative to occupied volume by total clusters, separated by sections of size of clusters: the first 2 largest clusters, the 3rd to 10% largest clusters, 10% to 30% largest clusters, 30% to 50% largest clusters, and the remaining 50%.
In Fig. 8 we present the effective fractal dimension of the largest cluster using a random walk algorithm. Lower than 4 dimensional structure of topological charge density can be found. As q cut increases, the effective fractal dimension drops, but even with q cut at 0.20q max −0.25q max , when percolation is just starting, global clusters with lower effective dimension d * = 1−1.4 have already appeared. This is consistent with the result from all-scale q(x) shown in Fig. 8(d). This structure was first found as a one-dimensional skeleton in Ref. [15].
The trends of the properties of clusters with respect to q cut /q max from our work shown on the left-hand side of Figs. 6-8 are consistent with the trends shown on the right-hand sides in the corresponding centers, which we quote from Ref. [13]. Most of the quantitative values are also compatible, except for the cluster multiplicity. Especially when q cut =0, dozens of clusters still survive, while in all-scale q(x) there are only two global clusters with opposite sign. Also, the magnitude of cluster multiplicity in Fig. 6(a) is about twice that in Fig. 6(b), which shows a scene much more chaotic than all-scale q(x) does. This extra chaos or clusters might come from the multiprobing approximation method itself. After doing more study on the data and considering that Fig. 6(c) is consistent with Fig. 6(d), we find that these extra small clusters add to the total number of clusters, but they do not substantially affect the fraction of the largest cluster to the total occupied volume. Especially at q cut = 0, there are few small fragments composed of only one or two sites, and the total volume of these fragments can be ignored with respect to the whole lattice space or the two largest clusters with opposite sign, which are the clusters that all-scale q(x) only have at q cut =0.
Note that in Figs. 6(b), 6(d), 7(b), 7(d), lattice ensembles at β = 8.10 with 53 configurations, at β = 8.45 with 5 configurations, and 2 configurations at β = 8.60, the numbers of configurations are so small because the direct all-scale q(x) calculation consumes too much computer resources. All the 3 ensembles used in this work consist of 200 configurations.
In Fig. 9 we present the fractional size, which is separated by different sections of size of clusters, relative to the total volume occupied by clusters. The clusters are ranked according to the volume of clusters. From Fig. 9, most volume is occupied by the first 50% of largest clusters in the three different ensembles below the onset of the percolation q cut = 0.2q max , then the fractional size of extra clusters should not be a concern. According to Figs. 6-8, however, above the onset of percolation the extra clusters from the multi-probing approximation still do not influence the role of the largest clusters substantially.
It is easy to understand why the multi-probing method will bring extra clusters. Because we use the clusters of topological charge density as tools to investigate the structures of the quenched QCD vacuum, with the hope that the structures would reflect ordered properties of q(x), and the multi-probing method will bring superpositions from different probing sites to the topological charge density, which can be regarded as disordered disturbance, then noise will be introduced to the ordered structure of clusters.
Note that the properties of clusters of topological charge density from the multi-probing method shown above are much closer to the all-scale topological charge density than those from the eigenmode expansion shown in Figs. 20 and 21 from Ref. [13].

Random permutation of topological charge density
Though we have multi-probing observed structures of topological charge density, there is still the question of where these stuctures come from. For example, do these structures arise from the ordered properties of the quenched QCD vacuum or can they arise from disordered vacuum and eventually reveal nothing? To get some information about this, we permute the topological charge density randomly through space-time to gen-erate stochastically distributed q(x), which should no longer contain information of the ordered properties of the quenched QCD vacuum, and then investigate the consequences of this permutation.   Figure 10 shows the results of the permuted topological charge density. The cluster multiplicity becomes much larger than in Fig. 6(a). Even at q cut = 0 there are still hundreds of clusters. The onset of percolation is found near 0.15q max now, and the fractional size of the largest cluster and the connectivity drop quickly. The distance is √ 2a, which is less than the results above in Fig. 7(a). These results suggest that the clusters tangle with each other more closely. Above the onset of percolation, the distance between the two largest clusters drops, when it reached a plateau before. The effective dimension shown in Figs. 10(e), 10(f) and 10(g) drops faster than those shown in Figs. 8(a), 8(b) and 8(c) and becomes negligible at q cut = 0.4q max . On the other hand, the changes shown in Fig. 10 also demon-strate that the differences between the results of multiprobing and all-scale q(x) seen above really come from disordered disturbance on the lattice. Therefore we can be sure that the properties of clusters from multi-probing method do reflect the underlying ordered structure of QCD vacuum, especially the structure recognized between q cut = 0.15q max −0.20q max where the structure of clusters is similar to all-scale q(x) but different from the permuted case.

Summary
In this paper, using the multi-probing method with overlap fermions to approach the topological charge density, the topological aspects of quenched QCD were inves-tigated. We studied the topological charge, topological susceptibility and the structure of topological charge density in space-time. It is found that the cluster properties are consistent with the results obtained from allscale topological charge density. After comparing them with clusters of randomly permuted topological charge, we conclude that the global structure that appears near the onset of percolation with effective dimension near 1 may well reflect the underlying ordered properties of the quenched QCD vacuum, since there are still global clusters surviving when the sites that have topological charge density of absolute value below q cut have been discarded. The q mp (x) approximated by the multi-probing method has been proved to deviate little from the exact value q(x) even when q cut is much lower than the onset of percolation, therefore the structures found near the onset of percolation are reliable.
This work was mainly run on the Tianhe-2 supercomputer at NSCC in Guangzhou.