Chiral symmetry and taste symmetry from the eigenvalue spectrum of staggered Dirac operators

We investigate general properties of the eigenvalue spectrum for improved staggered quarks. We introduce a new chirality operator $[\gamma_5 \otimes 1]$ and a new shift operator $[1 \otimes \xi_5]$, which respect the same recursion relation as the $\gamma_5$ operator in the continuum. Then we show that matrix elements of the chirality operator sandwiched between two eigenstates of the staggered Dirac operator are related to those of the shift operator by the Ward identity of the conserved $U(1)_A$ symmetry of staggered fermion actions. We perform a numerical study in quenched QCD using HYP staggered quarks to demonstrate the Ward identity. We introduce a new concept of leakage patterns which collectively represent the matrix elements of the chirality operator and the shift operator sandwiched between two eigenstates of the staggered Dirac operator. The leakage pattern provides a new method to identify zero modes and non-zero modes in the Dirac eigenvalue spectrum. This method is as robust as the spectral flow method but requires much less computing power. Analysis using a machine learning technique confirms that the leakage pattern is universal, since the staggered Dirac eigenmodes on normal gauge configurations respect it. In addition, the leakage pattern can be used to determine a ratio of renormalization factors as a by-product. We conclude that it might be possible and realistic to measure the topological charge $Q$ using the Atiya-Singer index theorem and the leakage pattern of the chirality operator in the staggered fermion formalism.

We investigate general properties of the eigenvalue spectrum for improved staggered quarks. We introduce a new chirality operator [γ5 ⊗ 1] and a new shift operator [1 ⊗ ξ5], which respect the same recursion relation as the γ5 operator in the continuum. Then we show that matrix elements of the chirality operator sandwiched between two eigenstates of the staggered Dirac operator are related to those of the shift operator by the Ward identity of the conserved U (1)A symmetry of staggered fermion actions. We perform a numerical study in quenched QCD using HYP staggered quarks to demonstrate the Ward identity. We introduce a new concept of leakage patterns which collectively represent the matrix elements of the chirality operator and the shift operator sandwiched between two eigenstates of the staggered Dirac operator. The leakage pattern provides a new method to identify zero modes and non-zero modes in the Dirac eigenvalue spectrum. This method is as robust as the spectral flow method but requires much less computing power. Analysis using a machine learning technique confirms that the leakage pattern is universal, since the staggered Dirac eigenmodes on normal gauge configurations respect it. In addition, the leakage pattern can be used to determine a ratio of renormalization factors as a by-product. We conclude that it might be possible and realistic to measure the topological charge Q using the Atiya-Singer index theorem and the leakage pattern of the chirality operator in the staggered fermion formalism.

I. INTRODUCTION
It is important to understand the low-lying eigenvalue spectrum of the Dirac operator, which exhibits the topological Ward identity of the Atiya-Singer index theorem [1], the Banks-Casher relationship [2], and the universality of the distribution of the near-zero modes for fixed topological charge sectors [3,4]. Study on the eigenvalue spectrum of the Dirac operator is, by nature, highly nonperturbative. Hence, numerical tools available in lattice gauge theory provide a perfect playground to study diverse properties of the Dirac eigenvalue spectrum.
In lattice QCD, there are a number of popular methods to implement a discrete version of the continuum Dirac operator. We are interested in one particular class of lattice fermions that are widely used in the lattice QCD community: improved staggered quarks [5][6][7]. Here we study the eigenvalue spectrum of staggered Dirac operators in quenched QCD to show that the small eigenvalues near zero modes of the staggered Dirac operators * E-mail: wlee@snu.ac.kr reproduce the continuum properties very closely, which was originally noticed in Refs. [8][9][10]. To reach this conclusion, the authors of Refs. [8,9] performed a number of tests, verifying consistency of lattice data with (1) the Atiya-Singer index theorem that describes the chiral Ward identity relating the zero modes to the topological charge; (2) the Banks-Casher relationship that relates the chiral condensate to the density of eigenvalues at the zero modes; and (3) the universality of the small eigenvalue spectrum in the ε-regime predicted by random matrix theory. In addition, the authors of Refs. [11,12] used the spectral flow method of Adams [13] to identify the zero modes from a mixture with non-zero modes. The spectral flow method is robust but highly expensive in a computational sense.
Here, we introduce a new, advanced chirality operator [γ 5 ⊗ 1], which respects the continuum algebra of γ 5 . We show that matrix elements of this chirality operator between eigenstates are related to those of the shift operator [1 ⊗ ξ 5 ] through the Ward identity of the conserved U (1) A symmetry of staggered fermions. In addition, we introduce a new concept of leakage patterns to distinguish zero modes from non-zero modes. Using the leakage pattern of the chirality and shift operators, we show that one can measure the zero modes as reliably as when using the spectral flow method. Hence, one could determine the topological charge Q using the leakage pattern with much smaller computational cost than by using the spectral flow. We also show that it is possible to determine the ratio of renormalization constants Z P ×S /Z P ×P using the leakage pattern.
In Section II, we briefly review the continuum theory of the eigenvalue spectrum and its relation to the quark condensate ψ ψ . We also review the Atiya-Singer index theorem in brief. In Section III, we briefly review the eigenvalue spectrum of staggered Dirac operators that is obtained using the Lanczos algorithm. In Section IV, we briefly review the conserved U (1) A symmetry in the staggered fermion formalism and explain its role in the eigenvalue spectrum of staggered Dirac operators. We also present numerical examples to help readers to understand basic concepts and notation. In Section V, we define the chirality operator [γ 5 ⊗ 1] and the shift operator [1 ⊗ ξ 5 ]. We show that they respect the continuum recursion relation of γ 5 . Then we derive the chiral Ward identity of the U (1) A symmetry to show that the matrix elements of the chirality operator are related to those of the shift operator through the Ward identity. We discuss the eigenvalue spectrum in the continuum limit and introduce a new notation of quartet indices. Then we introduce the concept of leakage patterns for the chirality operator and the shift operator. We also present numerical examples to demonstrate that the leakage patterns of zero modes are completely different from those of non-zero modes. In Section VI, we review a machine learning technique and describe how to apply it to extract efficiently the quartet structure of non-zero modes using leakage patterns. In Section VII, we explain how the leakage pattern of the zero modes can be used to determine the ratio of the renormalization factors nonperturbatively. In Section VIII, we conclude. The appendices contain technical details on Lanczos algorithms and mathematical proofs, and more plots of leakage patterns for diverse topological charge values.

II. QUARK CONDENSATE IN THE CONTINUUM
In the continuum the quark condensate is given by where D is the Dirac operator, m is the quark mass, x is the space-time coordinate, V is the volume, and N f is the number of flavors with the same mass m. The trace is a sum over spin and color. Let us think of the eigenvalues of the Dirac operator. D is anti-Hermitian, so its eigenvalues are purely imaginary or zero.
where λ is a real eigenvalue, and u λ (x) is the corresponding eigenvector. By spectral decomposition [4], where we adopt the normalization convention Thanks to the chiral symmetry, Hence, if there exists u λ for λ = 0, then the parity partner eigenstate u −λ with negative eigenvalue −iλ must also exist. Now let us separate the zero mode contribution from the spectral decomposition.
Here, n + (n − ) is the number of right-handed (lefthanded) zero modes per flavor. Let us define the subtracted quark condensate ψ ψ sub : where the spectral density ρ s (λ) is Here, ρ s is a spectral density on a single gauge configuration with volume V . Now let us average over a full ensemble of gauge field configurations and take the limit of infinite volume (V → ∞). Then, in that limit, the spectral density ρ(λ) = ρ s (λ) has a well defined (smooth and continuous) value as λ → 0. We can define the chiral condensate as which is the Banks-Casher relation. The subtracted quark condensate ψ ψ sub is expected to behave well in the chiral limit, even though the contribution from the zero modes is divergent as a simple pole in the chiral limit. Hence, in the numerical study on the lattice, it is important to identify the would-be zero modes which correspond to the zero modes in the continuum limit, and to remove them in the calculation of the quark condensate. Before proceeding, let us briefly go through the index theorem. In the continuum theory in Euclidean space, the axial Ward identity [17] is Here A µ ≡ψγ µ γ 5 ψ is the axial vector current in the flavor singlet representation, P ≡ψγ 5 ψ is the corresponding pseudo-scalar operator, and q ≡ 1 32π 2 Tr[F µν F µν ] is the topological charge density (= winding number density). Now the topological charge Q is Using the spectral decomposition, we can rewrite Q as follows, Noting that Hence, only zero modes contribute to Q. For the zero modes, it is convenient to choose the helicity eigenstates as the basis vectors so that u L 0 |γ 5 |u L 0 = −1 and u R 0 |γ 5 |u R 0 = +1, where the superscripts L, R represent left-handed and right-handed helicity, respectively. Then deriving the index theorem is straightforward [1]: where n + (n − ) is the number of the right-handed (lefthanded) zero modes.

III. SPECTRAL DECOMPOSITION WITH STAGGERED FERMIONS
A number of improved versions of staggered fermions exist, such as HYP-smeared staggered fermions [5], asqtad improved staggered fermions [18], and highly improved staggered quarks (HISQ) [7]. Here we refer to all of them collectively as "staggered fermions." Staggered fermions have four tastes per flavor by construction [19]. Hence, the quark condensate for staggered fermions is defined as where χ represents a staggered quark field, D s is the staggered Dirac operator for a single valence flavor, V is the lattice volume, and N t is the number of tastes. We measure the quark condensate using a stochastic method.
where x, y are representative indices which represent the space-time coordinates, taste, and color indices collectively. Here ξ(x) represents either Gaussian random numbers or U (1) noise random numbers which satisfy a simple identity: where N ξ is the number of random vector samples. Staggered fermions have a taste symmetry SU (4) L ⊗ SU (4) R ⊗ U (1) V in the continuum limit at a = 0 [20]. However, this symmetry breaks down to a subgroup U (1) V ⊗ U (1) A on the lattice with a = 0 [19,20]. The remaining axial symmetry U (1) A plays an important role in protecting the quark mass from receiving an additive renormalization. In addition, it does not have any axial anomaly.
The Dirac operator (D s ) of staggered fermions is anti-Hermitian: D † s = −D s . Hence, its eigenvalues are purely imaginary: where λ is real. Here, the subscript s and superscript s represent staggered quarks.
In practice, when we obtain eigenvalues of D s numerically, we use the following relationship instead of Eq. (28): where the |g s λ 2 state is a mixture of the two eigenvectors |f s +λ and |f s −λ . In other words, |g s λ 2 = c 1 |f s +λ + c 2 |f s −λ (30) where the c i are complex numbers that satisfy the normalization condition The numerical algorithm is a variation of a Lanczos algorithm adapted for lattice QCD [21]. Details on the numerical algorithms as well as comprehensive references are given in Appendix A. Why do we obtain λ 2 instead of iλ? The first reason is that doing so allows us to use even-odd preconditioning [22], which makes Lanczos run on only even or odd sites on the lattice. This leads to two benefits: One is that there is a substantial gain in the speed of the code, and the other is that the code uses only half the memory otherwise required. Details on even-odd preconditioning are described in Appendix B. The second reason is that obtaining λ 2 instead of iλ allows us to implement polynomial acceleration algorithms [23] into Lanczos more easily, since the eigenvalues of D † s D s are positive definite and have a lower bound λ 2 > 0. Note that staggered fermions have would-be zero modes whose eigenvalues are small and positive (λ 2 > 0) on rough gauge configurations. There are no exact zero modes (λ = 0) with staggered fermions on rough gauge configurations [24]. Details of our implementation of polynomial acceleration are described in Appendix A.
Hence, we use the Lanczos algorithm to solve Eq. (29) for the eigenvector |g s λ 2 as well as the corresponding eigenvalue λ 2 . We obtain |f s +λ and |f s −λ by using projection operators defined as where P + is the projection operator that selects only the |f s +λ component and removes the |f s −λ component. Then and the orthonormal eigenvectors are

IV. CHIRAL SYMMETRY OF STAGGERED FERMIONS
The two vectors |f s ±λ are related to each other through a chiral Ward identity of staggered fermions. Here we address this issue of the chiral symmetry of staggered fermions and its consequences.

A. Notation and definitions
Let us begin with notation and definitions. For staggered fermions, there are two independent methods to transcribe operators to the lattice: One is the Golterman method [19,25,26], and the other is the Kluberg-Stern method [27][28][29][30]. In Appendix D, we explain how to construct chirality operators using both the Golterman method and the Kluberg-Stern method, and we compare the two methods. The comparison is summarized in Table XI of Appendix D. Since the Kluberg-Stern method respects the recursion relationship, uniqueness of chirality, and the Ward identity while the Golterman method does not, we adopt the former method to construct bilinear operators. Accordingly we define staggered bilinear operators as where χ b are staggered quark fields, and a, b are color indices. Here the coordinate x A = 2x+A, x is a coordinate of the hypercube, and A, B are hypercubic vectors with A µ , B µ ∈ {0, 1}. The spin-taste matrices are where γ S represents the Dirac spin matrix, and ξ T represents the 4 × 4 taste matrix. In addition, where P SU (3) represents the SU (3) projection, and C represents the complete set of shortest paths from x A to x B . V (x, y) represents the HYP-smeared fat link [5,6] for HYP staggered fermions, the Fat7 fat link [6,[31][32][33] for asqtad staggered fermions or HISQ, and the thin gauge link for unimproved staggered fermions. The conserved U (1) A axial symmetry transformation is where Γ is often called "distance parity," and Under the U (1) A transformation, the staggered Dirac operator transforms as follows, Therefore, and f s −λ can be obtained from f s +λ through Γ transformation as follows, In general, there is no constraint for the real phase θ, so we expect its probability distribution to be random. In practice, however, we make use of even-odd preconditioning, and we obtain the odd site fermion fields (|g o ) from the even site fermion fields (|g e ) with the relation |g o = η D oe |g e , where D oe is a sector of D s that connects even site fields to odd site fields, and η is a random complex number. Hence the distribution of θ depends on our choice of η. In our numerical study, we set η = 1. Then θ is given by Details on the even-odd preconditioning and the derivation of Eq. (48) are explained in Appendix B. We expect that if there exists an eigenvector |f s +λ , there must be a corresponding parity partner |f s −λ due to the exact chiral symmetry Γ . In other words, the Ward identity of Eq. (47) comes directly from the conserved U (1) A axial symmetry.

B. Numerical Examples
We now use numerical examples to demonstrate how the above theory works in quenched QCD. In Table I, details of the gauge configurations are presented.
We measure the topological charge Q using gauge links. We use the Q(5Li) operator defined in Ref. [40,41] after 10 ∼ 30 iterations of APE smearing with α = 0.45 [42][43][44]. We show an example of the eigenvalue spectrum for Q = 0 in Fig. 1. Since Q = 0, we do not expect to find any zero modes for this gauge configuration. In Fig. 1 (a), we show eigenvalues λ 2 for the eigenvectors |g s λ 2 defined in Eq. (29). Here we observe eight-fold degeneracy for non-zero eigenmodes due to the conserved U (1) A axial symmetry. Here λ 2 = −λ 1 and, in general, λ 2n = −λ 2n−1 for integer n > 0. In other words, λ 2n is the parity partner of λ 2n−1 . For each λ i , there exists four-fold degeneracy due to approximate SU (4) taste symmetry. For each of these four-fold degenerate eigenvalues (for example λ 1 , λ 3 , λ 5 , λ 7 in Fig. 1 (a)), there  exists a parity partner eigenvalue due to the U (1) A symmetry: λ 2 = −λ 1 , λ 4 = −λ 3 , λ 6 = −λ 5 , and λ 8 = −λ 7 (refer to Fig. 1 (b)). Let us turn to an example with Q = −1. Since Q = −1, we expect to observe four-fold would-be zero modes. The gauge configurations are so rough that we expect to observe not exact zero modes but would-be zero modes. In Fig. 2, we demonstrate how the would-be zero modes behave on a gauge configuration with Q = −1. As one can see in Figs. 2 (a) and 2 (b), we find four-fold degenerate would-be zero modes: λ 1 , λ 2 , λ 3 , λ 4 . Thanks to the U (1) A chiral Ward identity in Eq. (47), we find that λ 2 = −λ 1 and λ 4 = −λ 3 . As in the case with Q = 0, we find that the non-zero eigenmodes are eight-fold degenerate. This pattern of four-fold degeneracy for would-be zero modes and eight-fold degeneracy for non-zero modes is also observed in the cases with Q = −2 and Q = −3, which are presented in Appendix C. At this point, the reader may have already concluded that we can distinguish would-be zero modes of staggered quarks from non-zero modes by counting the degeneracy of the eigenvalues [8,9,45]. This conclusion is true but can lead to wrong answers in practice. The reason is that, on large lattices, the eigenvalues are so dense that visually distinguishing between 4-fold and 8-fold degeneracies is typically impossible. Hence, we need a practical and robust method to identify would-be zero modes and non-zero modes of staggered fermions. The introduction of such a method is the main subject of the next section, Sec. V.
Using the chiral Ward identity of Eq. (47), we can measure the phase θ numerically. In Fig. 3, we show numerical results (red circles) for θ. Here, the blue line represents the theoretical prediction given in Eq. (48). We find the results are consistent with the prediction within numerical precision.

V. CHIRALITY MEASUREMENT
To simplify the notation, we introduce the following convention for eigenvalue indices, where |f j = |f s λj is defined in Eq. (28). Using the Kluberg-Stern method explained in Appendix D, we define the chirality operator and abbreviations as follows.
The phase θ as a function of λ. Red circles represent numerical results for θ. The blue line represents the prediction from the theory. Here we use a gauge configuration with Q = −1 for the measurement.
where x A , x B and [γ 5 ⊗ 1] are defined in Eqs. (38)- (40) in Subsection IV A, and λ i and λ j represent eigenvalues of D s . The chirality operator [γ 5 ⊗ 1] satisfies the same relationships as the continuum chirality operator γ 5 . [ where n is a non-negative integer. A rigorous proof of Eqs. (53)-(56) is given in Appendix E. Our definition of the chirality operator [γ 5 ⊗1] uses the Kluberg-Stern method, and is different from that used in Refs. [8,13,24], which adopt the Golterman method. The old chirality operator (the Golterman method) of Refs. [8,13,24] does not satisfy the recursion relations of Eqs. (53)- (56). In addition, it does not satisfy the chiral Ward identities of Eqs. (62)-(64). This difference is addressed in Appendices D and E. The bottom line is that the conventional chirality operator (the Golterman method) does not satisfy the recursion relationships in Eqs. (53)-(56), even though it is classified according to the true irreducible representations of the lattice rotational symmetry group [25,26,28].
For further discussion we define another operator [1 ⊗ ξ 5 ], which we call the "(maximal) shift operator," where x A , x B and [1 ⊗ ξ 5 ] are defined in Eqs. (38)- (40) in Subsection IV A, and λ i and λ j represent eigenvalues of D s . This shift operator satisfies the following recursion relations: where n is a non-negative integer. The conserved U (1) A symmetry transformation can be expressed in terms of the chirality operator and the shift operator as follows, A rigorous proof of Eq. (62) is given in Appendix E.
In addition, the chirality and shift operators satisfy the following relations: A rigorous proof of Eqs. (63)-(64) is also given in Appendix E. Therefore, we obtain the Ward identities: where Hence, we define the spectral decomposition by where we use Eqs. (50) and (51). Similarly, where we use Eqs. (57) and (58). Thanks to the Ward identities of Eq. (65), we obtain Similarly, Applying Γ to both sides of Eq. (67), we obtain Hence, we obtain another Ward identity: Similarly, we can obtain the Ward identities: We can summarize all the results of Eqs. (69)-(75) in the following form: In addition, Hermiticity ensures that we can interchange λ i and λ j , which gives the final form of the chiral Ward identities.
The quantity (|Γ 5 | i j ) 2 for i = j represents the leakage probability of the chirality operator. We call |Γ 5 | i j the leakage parameter for the chirality operator. Similarly, the quantity (|Ξ 5 | i j ) 2 for i = j represents the leakage probability of the shift operator, and we call |Ξ 5 | i j the leakage parameter for the shift operator. By examining the leakage pattern, we can distinguish zero modes and non-zero modes, which is the main subject of the next subsections.

A. Eigenvalue spectrum of Ds in the continuum
Here we consider staggered quark actions at a = 0. We define a general form of the shift operator corresponding to a generator of the SU (4) taste symmetry: where ξ µ respects the Clifford algebra {ξ µ , ξ ν } = 2δ µν in Euclidean space. Let us consider the following quantity W 1 in the continuum: Since the SU (4) taste symmetry is exactly conserved in the continuum, we know that Hence, we find the following Ward identity: and therefore Hence, in the continuum, we find the following properties of the eigenvalue spectrum: In other words, if the eigenvalues are different, there is no leakage ((Ξ F ) n = 0) between the two eigenmodes.
• If λ j ≡ λ = λ n , (Ξ F ) n = 0 is possible. In this case, the eigenvalues are degenerate and belong to a quartet such that where m is a taste index which represents the fourfold degeneracy for the eigenvalue λ j , |f and |f n are linear combinations of the quartet {|f j,m } , and the eigenvectors for different m are orthogonal to each other by construction due to the Lanczos algorithm.
• We know that the staggered fermion field χ c (x A ) is mapped into the continuum fermion field ψ c α;t (x), where α represents a Dirac spinor index, c represents a color index, and t = 1, 2, 3, 4 represents a taste index. Hence, for a given eigenvalue λ j , there remain four degrees of freedom which come from the taste index. Accordingly, for a given eigenvalue λ j , there are four degenerate eigenstates |f j,m with m = 1, 2, 3, 4.
• If we know all four eigenstates {|f j,m } for a certain eigenvalue λ j , we find that This is because the SU (4) group generators are traceless in the fundamental representation.
However, on the lattice at a = 0, the taste symmetry is broken by terms of order a 2 α n s with n ≥ 1 , which is explained in Ref. [20]. In addition, for a = 0, and λ j,m = λ j,m in general for m = m , which reflects the taste symmetry breaking effect at a = 0. We know that λ j,m = λ j,m for all m, m in the continuum (a = 0) due to the exact taste symmetry. Hence, on the finite lattice, we expect a small deviation from the above continuum properties. A good barometer to measure this effect is to monitor T 5 and measure how much it deviates from zero (the continuum value). Another direct barometer to measure effects of taste symmetry breaking is the leakage S 5 from one quartet (λ ) to another quartet (λ j ) with λ = λ j .
The size of the leakage S 5 indicates directly how much the taste symmetry is broken at a = 0, since S 5 = 0 in the continuum. We present numerical results for T 5 and S 5 in the next subsection.

B. Numerical study on chirality and leakage
Here we use dual notations for the eigenmodes: One is the serial index i for λ i and the other is the quartet index j with taste index m for λ j,m . The serial index is convenient for the plots, tables, and leakage patterns such as |Γ 5 | a b , while the quartet index is convenient to explain the eigenstates classified by the taste symmetry group. The one-to-one mapping from the serial index i to the quartet indices j, m is given in Table II for the quartet index j = 0, ±1 when Q = ±1. The mapping for the quartet index j = ±2 (non-zero modes) is similar.
TABLE II. One-to-one mapping of serial index i of the λi eigenstate into a quartet index j and taste index m for λj,m = λi. Here λ2n = −λ2n−1 and λ−j,m = −λ+j,m. Here we assume Q = ±1.
In Fig. 4, we present the leakage pattern of the zero mode of λ 1 and its parity partner λ 2 = −λ 1 . Since Q = −1 in Fig. 4, we expect to observe four-fold degenerate would-be zero modes within a single quartet (quartet index j = 0).
In the continuum limit (a = 0), the SU (4) taste symmetry becomes exactly conserved and so would-be zero modes become exact zero modes. However, at finite lattice spacing a = 0, the gauge configuration is sufficiently rough that would-be zero modes have non-zero eigenvalues: λ 2 = −λ 1 , λ 4 = −λ 3 , and λ 1 = λ 3 for λ 1 , λ 3 > 0. In Fig. 4 (a), we show the leakage pattern of We find that there is, in practice, no leakage, and so the only non-zero component is The other components are practically zero. In Figs. 4 (b), 4 (c), and 4 (d), we find that the Ward identity of Eqs. (78) and (79) is well-respected by the numerical results. In other words, the Ward identity is satisfied within the numerical precision of the computer. Please refer to Table  III for numerical details. This confirms that the theoretical prediction from the Ward identity in Eqs. (78) and (79) is correct.
In Fig. 4 (a), we find that there is small leakage into other quartets (j = ±1, ±2). The typical size of leakage between off-diagonal elements of the would-be zero modes, the j = 0 quartet, (e.g. |Γ 5 | 3 1 ) is of order 10 −3 . We also observe small leakage patterns of order Leakage pattern for would-be zero modes at Q = −1.
Here, the red bar represents leakage to λi=2n−1 > 0 with i odd, and the blue bar represents leakage to its parity partner λi=2n = −λ2n−1 with i even.
10 −2 ∼ 10 −3 from the would-be zero modes, the j = 0 quartet, to the non-zero modes, the j = ±1, ±2 quartets (e.g. |Γ 5 | 5 1 ). Now let us consider non-zero modes in the j = +1 quartet. In Fig. 5, we present the leakage pattern for the non-zero modes of λ 5 and its parity partner λ 6 = −λ 5 . Even in the continuum limit (a = 0), λ 5 = 0, and so it is a non-zero mode. Thanks to the approximate SU (4) taste symmetry and the exact U (1) A axial symmetry, there will be eight-fold degeneracy in the family of eight eigenstates composed of the j = +1 quartet, to which λ 5 belongs, and the j = −1 quartet (the parity partners). These eight-fold degenerate modes are designated together as j = ±1 quartets in Fig. 5, where they are a set of {λ i } with 5 ≤ i ≤ 12.  Fig. 4. Here, j represents a quartet index for the λi eigenstate. The leakage represents Let us scrutinize the leakage pattern of the non-zero mode λ 5 = λ j=+1,m=1 . In Fig. 5 (a), first note that there is practically no leakage in the Γ 5 chirality measurement from λ 5 into λ 2n−1 with n > 0. In other words, |Γ 5 | 2n−1 5 = |Γ 5 (λ 2n−1 , λ 5 )| ∼ = 0. This implies that the chirality operator on the non-zero mode with λ > 0 leaks into only the parity partner modes with λ < 0. Second, note that the nontrivial leakage goes to those eigenstates in the j = −1 quartet such as  Table IV, we present numerical values of the |Γ 5 | i 5 shown in Fig. 5 (a) . Let us examine the Γ 5 = [γ 5 ⊗ 1] leakage pattern of the j = +1 quartet of the non-zero modes {λ 5 , λ 7 , λ 9 , λ 11 }. In Fig. 6, we find that the chirality measurement vanishes; (Γ 5 ) i i = Γ 5 (λ i , λ i ) = 0 for λ i in the j = +1 quartet of non-zero modes. We also find that the Γ 5 leakage of λ +1,m > 0 of the j = +1 quartet goes to the parity partners with λ −1,m < 0 of the j = −1 quartet, and the leakage to other quartets such as j = ±2 is negligibly small compared to the leakage to the j = −1 quartet. The numerical values of |Γ 5 | −1,m +1,m are summarized in Table V.
Let us examine the Ξ 5 = [1 ⊗ ξ 5 ] leakage pattern of the j = +1 quartet of the non-zero modes {λ 5 , λ 7 , λ 9 , λ 11 }. In Fig. 7, we find that the Ξ 5 leakage from the j = +1 quartet to the j = −1 quartet (parity partners) vanishes in practice. Since the leakage pattern of Ξ 5 is related to the leakage pattern of Γ 5 by the Ward identity Fig. 7 can be obtained from Fig. 6 using the Ward identity. We find that the Ξ 5 leakage from the j = +1 quartet to other quartets such as j = ±2 quartets is negligibly small compared to leakage to itself (the j = +1 quartet). Leakage patterns of the Γ 5 chirality and Ξ 5 shift operators for diverse topological charges are shown in Appendix G. Let us summarize the leakage pattern for would-be zero modes and that for non-zero modes. We first begin with the leakage pattern for the zero modes.   into its parity partner mode with eigenvalue −λ, and no leakage into any other eigenmodes.
The leakage pattern for nonzero modes is 1. A non-zero mode of staggered fermions appears as an eight-fold degeneracy composed of a quartet (+j quartet) and its parity partner quartet (−j quartet). In other words, non-zero eigenmodes can be grouped into sets with eight elements in each set. This is due to the approximate SU (4) taste symmetry and the conserved U (1) A axial symmetry.
2. In the chirality Γ 5 = [γ 5 ⊗1] measurement, the nonzero mode with eigenvalue λ j,m has no leakage to its own quartet (j quartet), but has leakage only to the parity partner (−j quartet) with λ −j,m . It has no leakage to any eigenmode which belongs to other quartets with = ±j.
3. In the shift Ξ 5 = [1⊗ξ 5 ] measurement, the non-zero mode with λ j,m has no leakage to its parity partner (−j quartet) at all. But it has leakage only to the eigenstates in its own (+j) quartet. This pattern comes directly from the Ward identity. In other words, the Ξ 5 leakage pattern is a mirror image reflecting Γ 5 through the mirror of the Ward identity. Ξ 5 has no leakage to any eigenmode which belongs to other quartets with = ±j. is identical to that of |Ξ 5 | +j,m ,m . In Appendix F, we provide more examples to demonstrate our claim that the leakage pattern for zero modes holds in general. In Appendix G, we give more examples to demonstrate our claim that the leakage pattern for non-zero modes is valid in general. We have repeated numerical tests over hundreds of zero modes and tens of thousands of nonzero modes. We performed the numerical study on hundreds of gauge configurations and find that the above leakage pattern is valid for all of them except for those gauge configurations with unstable topological charge.
1. We find a number of gauge configurations which do not have a stable topological charge.
2. We find about 10 gauge configurations with unstable topological charge among 100 gauge configurations with 12 4 lattice geometry at β = 4.6.
3. We find about 8 gauge configurations with unstable topological charge among 300 gauge configurations with 20 4 lattice geometry at β = 5.0.  Table I. Nq represents the number of quartets used to obtain the statistical error. Here j = 0 represents would-be zero modes, and j > 0 represents non-zero modes. In Table VI, we present results for T 5 defined in Eq. (92), which is a direct barometer to estimate the effect of taste symmetry breaking. If the taste symmetry is exactly conserved, then T 5 must vanish. Hence, a nontrivial value of T 5 indicates the size of taste symmetry breaking. In Table VI, we find that |Re(T 5 )| is of order 10 −3 , while |Im(T 5 )| is essentially zero. This indicates that the effect of taste symmetry breaking is very small (in the sub-percent level within each quartet). In Fig. 8, we present S 5 defined in Eq. (93) as a function of | − j| with , j ≥ 0. Here | − j| = 1 represents a pair of nearest neighbor quartets, | − j| = 2 represents a pair of next-to-nearest neighbor quartets, and so on. The values of S 5 are the same size as the statistical errors. This indicates that the taste symmetry breaking results in simply random noise added to the physical signal (S 5 = 0). For | − j| = 1, the noise is ≈ 7%, and for | − j| = 2, the noise is ≈ 3%. We find that the noise decreases as | − j| increases. The numerical values of S 5 in Fig. 8 are presented in Table VII.  Here, we measure S5 between two different quartets ( = j and , j ≥ 0). Np represents the number of ( , j) pairs with = j.

VI. MACHINE LEARNING
In previous sections, we have shown that the U (1) A symmetry of staggered fermions induces the chiral Ward identities in Eq. (78), and we have also noted that the approximate SU (4) taste symmetry brings in the quar-  Table VII. tet behavior of the eigenvalue spectrum. Furthermore, a combined effect of those symmetries gives us distinctive leakage patterns for the chirality operator Γ 5 and the shift operator Ξ 5 . In this section, we apply a machine learning technique to the following tasks.
1. We want to know how much the non-zero modes respect the quartet classification rules, which come from the SU (4) taste symmetry.
2. We want to know how efficiently we can measure the topological charge Q using the index theorem from the quartet structure of the non-zero modes.
3. We want to detect any anomalous behavior of the eigenvalue spectrum, which does not follow the standard leakage pattern of the non-zero modes.

We want to figure out what causes the anomalous behavior of the eigenvalue spectrum.
Let us explain our sampling method for the machine learning. In Fig. 9, we show matrix elements |Γ 5 | i j on a gauge configuration with Q = 2. Fig. 9(a) is for the 200 lowest eigenmodes, and Fig. 9(b) is a zoomed-in version of Fig. 9(a) for the 32 lowest eigenmodes. Here the depth of the blue color represents the size of the matrix element |Γ 5 | i j , and i, j run from zero to 199. We identify two zero mode quartets (red boxes) by looking at the magnitude of the diagonal components. These two quartets have the same chirality (n − = 2), which is consistent with the topological charge Q = 2. Excluding the would-be zero modes, we randomly choose a 15 × 15 sub-matrix of |Γ 5 | i j along the diagonal line of |Γ 5 | i j matrix elements. This 15 × 15 sub-matrix is the smallest square sub-matrix of |Γ 5 | which contains all elements of only one quartet of non-zero modes and its parity partner quartet.
In Fig. 10, we present 8 different classes for arbitrary samples. Our purpose for the machine learning is to find borders (black lines) of the non-zero mode quartets (or octets when the parity partners are included) in each sample. We classify arbitrary samples into eight different classes according to the location of the border lines. Each class is labeled as in Fig. 10. We use a deep learning model which combines the multi-layer perceptron (MLP) [46] and the convolutional neural network (CNN) [46]. In Table VIII, we present our basic setup for the machine learning. We use the gauge configuration ensemble described in Table I  around ten 15 × 15 matrix samples from the 200 lowest eigenmodes without overlapping. We make popular and suitable choices for the loss function [49], optimization method [50], and activation functions [51] relevant to our purpose, which are summarized in Table VIII. The best hyper-parameters such as the number of layers and the number of units for each layer are determined by using the Keras Tuner [47]. The accuracy of classification per gauge configuration is obtained by averaging the accuracies of the machine learning (ML) prediction for all the samples on a single gauge configuration. Our best model achieves an average accuracy of 96.5(156)% for 142 test gauge configura-  [48] activation function for hidden layers ReLU [46] activation function for output layer softmax [46]  tions. The hyper-parameters which represent the structure of the neural network model are given in Table IX. Among the test set, we find five gauge configurations on which the average accuracy per gauge configuration is lower than 50%. Data show that some ghost (unphysical) eigenvectors are present in the eigenvalue spectrum on these gauge configurations, so that the ML prediction gives a wrong answer not due to failure of the ML algorithm but due to human mistakes in labeling quartet samples based on the eigenvalue index. Excluding these five gauge configurations, we achieve the average accuracy of 99.4(23)%. Considering that all samples generated on the same gauge configuration are connected by the eigenvalue index (or quartet index), this average accuracy of 99.4% implies that one can in the end find completely correct quartet groups for all the normal gauge configurations of the test set. It also demonstrates our claim that the leakage pattern is universal over all normal gauge configuration ensembles. Details of the results of this ML research will be reported separately in Ref. [52].

VII. ZERO MODES AND RENORMALIZATION
As explained in Sec. V, we know that there is practically no leakage for the zero modes in the chirality measurement. Hence, it is possible to determine the renormalization factor κ P by imposing the index theorem as follows. For Q = 0, where S 0 is the entire set of zero modes, and where and the subscript [· · · ] R ([· · · ] B ) represents a renormalized (bare) operator. The Z P ×S and Z P ×P are the renormalization factors for the bilinear operators O S and O P , respectively. One advantage of this scheme is that κ P is independent of valence quark masses, even though we perform the measurement with arbitrary masses for valence quarks. Numerical results for κ P are summarized in Table X.
There are a few key issues in the physical interpretation of κ P .
• Since the topological charge Q and sum C 0 are independent of renormalization scale, κ P must be independent of the renormalization scale µ.
• This means that the scale dependence of Z P ×S (µ) must cancel off that of Z P ×P (µ).
• It would be nice to cross-check this property of κ P in the RI-MOM scheme [53] and in the RI-SMOM scheme [54].

VIII. CONCLUSION
We study general properties of the eigenvalue spectrum of Dirac operators in the staggered fermion formalism.
As an example, we use the Dirac operator for HYP staggered quarks. In Section V, we introduce a new chirality operator Γ 5 and a new shift operator Ξ 5 and prove that they respect the continuum recursion relationships, as given in Eqs. (53)- (56) and Eqs. (60)- (61). Using these operators with nice chiral properties, we find that the leakage pattern of |Γ 5 | −j,m ,m is related to that of |Ξ 5 | j,m ,m through the Ward identity of the conserved U (1) A symmetry. We find that the leakage pattern of Γ 5 and Ξ 5 for the zero modes is quite different from that for the non-zero modes. This difference in leakage pattern allows us to distinguish the zero modes from the non-zero modes even though we do not know a priori about the topological charge. We find that using the leakage pattern of Γ 5 and Ξ 5 , one can determine the topological charge as reliably as when using standard field theoretical methods such as the cooling method.
We use a machine learning (ML) technique to check the universality of this leakage pattern over the entire ensemble of gauge configurations (refer to Table I). Our best-trained deep learning model identifies the quartet of non-zero modes with 98.7(34)% accuracy using a single normal gauge configuration. Choosing the highest probability prediction of the ML and comparing the prediction with the known answer, we find that the ML can identify all quartet groups on an eigenvalue spectrum correctly. In addition, the ML technique detects wrong answers resulting from human input mistakes since the ML prediction disagrees with a wrong answer by giving the prediction with low accuracy (< 50%). This reassures us that the ML technique is highly reliable at identifying anomalous gauge configurations with defects such as violation of the index theorem and ghost eigenmodes.
Once we identify the zero modes, it is also possible to determine the ratio of renormalization factors κ P = Z P ×S (µ)/Z P ×P (µ) from the chirality measurement of Γ 5 .
The leakage pattern is a new concept introduced in this paper. It can be used to study the low lying eigenvalue spectrum of staggered Dirac operators systematically. It helps us understand how to extract the taste symmetry and chiral symmetry from the staggered eigenvalue spectrum. The leakage pattern will help us to dig out related physics more efficiently, such as topological charge, index theorem, Banks-Casher relation, and nonperturbative renormalization.

ACKNOWLEDGMENTS
We would like to express sincere gratitude to Eduardo Follana for helpful discussion and providing his code to cross-check results of our code. We also thank Jon A. Bailey for helpful comments on the manuscript. We also thank Stephen R. Sharpe for helpful discussion. Lanczos is a numerical algorithm for calculating eigenvalues and eigenvectors of a Hermitian matrix [21]. It transforms an n × n Hermitian matrix H to a tridiagonal matrix T through a unitary transformation Q, which is represented by Here columns of Q are composed of basis vectors of the nth Krylov subspace K n (H, b) generated by H and a starting vector b of our choice. Each iteration of Lanczos computes a column of Q and T in sequence. At the end, diagonalizing the tridiagonal matrix T yields eigenvalues and eigenvectors of H. In principle, Lanczos is a direct method that takes n iterations to construct the n × n tridiagonal matrix T . However, since these columns of T are computed in order, a sequence of m < n iterations also constructs an m×m tridiagonal matrix T which is a submatrix of T . In practice, the real benefit of Lanczos is that eigenvalues of T approximate eigenvalues of T . As iteration continues, and the size of the submatrix T increases, eigenvalues of T converge to eigenvalues of T . The convergence behavior is somewhat complicated. The eigenvalues converge to the largest, the smallest, or the most sparse eigenvalue first. The speed of convergence depends on the density of eigenvalues. The less dense, the faster the convergence.
In this paper, we make use of two popular improvement techniques of Lanczos: (1) implicit restart [55], and (2) polynomial acceleration with Chebyshev polynomials [56]. The implicit restart method gets rid of converged eigenvalues in the middle of the Lanczos iteration. It takes effect as if we restarted the Lanczos with a shifted matrix H given by where λ i are eigenvalues we want to remove. Then H is still Hermitian but does not have such eigenvalues λ i . Hence, Lanczos with H converges to remaining eigenvalues faster. The implicit restarting procedure gives us a new submatrix, which has a dimension ((m−r)×(m−r)) reduced by the number of eigenvalues we have removed (r). Then we iterate Lanczos r times to refill the submatrix and restore the structure of the m × m matrix.
We repeat the implicit restart to obtain a new submatrix of dimension (m − r) × (m − r), and so on. This procedure allows us to control the size of the submatrix, the computational cost, and the memory usage, while the submatrix T contains (m − r) eigenmodes that are more precise (much closer to the true eigenmodes of the full matrix H) for each iteration. A polynomial operation on a matrix changes the eigenvalue spectrum accordingly while retaining the eigenvectors. Since the polynomial of a Hermitian matrix is also Hermitian, Lanczos is still available to calculate its eigenvalues and eigenvectors. By choosing a proper polynomial, one can manipulate the density of the eigenvalue spectrum so that the convergence to the desired eigenvalues is accelerated. A Chebyshev polynomial is a popular choice for this purpose. Using the Chebyshev polynomial, we want to map the first region of eigenmodes of no interest to [−1, 1] and map the second region of eigenmodes of our interest to [−∞, −1]. In the interval [−1, 1], the eigenvalues are dense enough that Lanczos does not converge. In addition, the Chebyshev polynomial rapidly changes in the second region so that the density of eigenmodes is low enough to more quickly accelerate the convergence of Lanczos. Here we apply the Chebyshev polynomial for D † s D s , whose eigenvalues are λ 2 ≥ 0. We set the lower bound of the first region to a value somewhat greater than the largest eigenvalue of interest. This strategy will not only suppress high unwanted eigenmodes, but also accelerate the speed of Lanczos for the low eigenmodes of interest.
Numerical stability is essential for the Lanczos algorithm. Each Lanczos iteration generates Lanczos vectors, which are column vectors of the unitary matrix Q in Eq. (A1). After several iterations, however, Lanczos vectors lose their orthogonality due to gradual loss of numerical precision. If not addressed, this loss would induce spurious ghost eigenvalues [57]. A straightforward prescription to solve the problem is performing a reorthogonalization for every calculation of Lanczos vectors. There are also alternative approaches to eliminate the ghost eigenvalues without reorthogonalization, such as the Cullum-Willoughby method [58,59]. Here we choose the first solution and perform the full reorthogonalization for each Lanczos iteration.
For a large scale simulation using Lanczos, Multi-Grid Lanczos [60] and Block Lanczos [61] are available. Multi-Grid Lanczos is also based on the implicit restart and Chebyshev acceleration. In addition, Multi-Grid Lanczos reduces the memory requirement significantly by compressing the eigenvectors using their local coherence [62]. A spatially-blocked deflation subspace is constructed from some of the lowest eigenvectors of the Dirac operator. Then the coherence of eigenvectors allows us to represent other eigenvectors on this subspace and to run Lanczos with much less memory. Meanwhile, Block Lanczos utilizes the Split Grid method [61]. This algorithm deals with multiple starting vectors for Lanczos, where the Split Grid method divides the domain of the Dirac operator application into multiple smaller domains so that each partial domain runs in parallel on a partial grid (lattice) with a lower surface to volume ratio than that of the full grid. Hence, one can optimize the offnode communication by adjusting the block (grid) size. This approach would give a significant speed-up compared with our method. We plan to implement Multi-Grid Lanczos and Block Lanczos in the near future.
Appendix B: Even-odd preconditioning and phase ambiguity Even-odd preconditioning reorders a fermion field χ(x) so that even site fermion fields are obtained first, and odd site fermion fields are obtained from them: where χ e (χ o ) is the fermion field collection on even (odd) sites. On this basis, the massless staggered Dirac operator D s can be represented as a block matrix: where D oe (D eo ) relates even (odd) site fermion fields to odd (even) site fermion fields. Since D † s = −D s , we also find that D † oe = −D eo and D † eo = −D oe . On this basis, D † s D s is expressed as Hence, the eigenvalue equation of D † s D s (Eq. (29)) can be divided into two eigenvalue equations as follows, where |g e(o) is the collection of even (odd) site components of |g s λ 2 . Here we omit the superscript s and the subscript λ 2 for notational simplicity. Now let us multiply D oe from the left on both sides of Eq. (B5). Then we find that which is identical to Eq. (B6). Hence, we find that |g o = η D oe |g e where η = re iα is an arbitrary complex number with r > 0 and 0 ≤ α < 2π. Here r represents the scaling behavior and α represents a random phase. Since −D eo D oe (= D † oe D oe ) is Hermitian and positive semi-definite, one can solve Eq. (B5) using the Lanczos algorithm introduced in Appendix A. From the result for |g e , it is straightforward to obtain the eigenvector |g s λ 2 of Eq. (29) since Now we apply the projection operator P + , defined in Eq. (32), to |g s λ 2 . Using Eq. (B5), we find that Similarly, for the projection operator P − , defined in Eq. (33), we find that Since η only appears in the overall factor for both cases, it gives only the relative phase difference between the normalized eigenvectors |f s ±λ defined in Eqs. (36) and (37).
We can proceed further to obtain the eigenvectors |f s ±λ . The norm of |χ + is given by Hence, |f s +λ is where Similarly, These results for |f s ±λ indicate that the phase difference θ for the Γ transformation defined in Eq. (47) depends on the value of η.
In our numerical study, we set η to η = re iα = 1: r = 1 and α = 0. Hence, the relative random phase between |f s ±λ states is removed by hand. Therefore, our value of θ defined in Eq. (47) includes a bias from our choice of η = 1.
where the coordinate x A = 2x+A, x is a coordinate of the hypercube, and A is a hypercubic vector with A µ ∈ {0, 1} with µ = 1, 2, 3, 4. M γ5⊗1 is defined as where p 1 = p 2 = p 3 = p 4 , and P is the set of all permutations of {1,2,3,4}. The symmetric shift operatorD µ is defined as whereμ is a unit vector in the µ direction of Euclidean space. V µ (y) is a (smeared) gauge link used in Eq. (40). ρ γ5⊗1 is defined as whereĀ µ = (A µ + 1) mod 2. In the Golterman method, the chirality operator [γ 5 ⊗ 1] Gol connects a site A=(0,0,0,0) ofχ to the 16 sites B=(1,1,1,1), (-1,1,1,1), ..., (-1,-1,-1,-1) of χ. As a consequence, unlike the continuum chirality operator which satisfies Using the definition of bilinear operators in Eq. (38), obtained with the Kluberg-Stern method, we define the chirality operator as Using the definition in Eq. (39) with γ S = γ 5 and ξ T = 1, we find Hence, the phase of the Kluberg-Stern operator is identical to that of the Golterman operator, which is in general true for the whole set of bilinear operators. We have freedom to choose U (x A , xĀ) ab to make the chirality operator gauge-invariant. Here, we set U (x A , xĀ) ab to where P SU (3) represents the SU (3) projection, and C represents the complete set of shortest paths from x A to xĀ.
Here, the SU (3) projection is crucial to make the chirality operator satisfy the continuum recursion relation: The Kluberg-Stern operator without the SU (3) projection P SU (3) contains the Golterman operator as a leading term: where O irrel represents irrelevant operators of higher dimension. For example, O irrel includes a 4-dimensional operator: where p 1 = p 2 = p 3 = µ, and P µ is the set of all permutations of {p a |p a = µ}. Technical details of the derivation of a complete set of irrelevant operators are explained in Ref. [63]. All the irrelevant operators have tastes different from 1 (ξ T = 1), and they contain at least one derivative D µ , which leads to higher dimension operators. As a consequence, their contribution to the chirality vanishes in the continuum limit a → 0. The recursion relation in Eq. (D9) is essential to define the chirality value uniquely for the staggered fermion formulation.
for all positive n ∈ Z. Hence, in the case of the Kluberg-Stern operators with the SU (3) projection, we can define the chirality value uniquely without any ambiguity. However, in the case of the Golterman operators, it is not possible to define the chirality value uniquely due to the following ambiguity: In addition, the Golterman operator does not satisfy the Ward identity, while the Kluberg-Stern operator respects it, A rigorous proof is given in Theorem E.3. However, in the continuum limit, they converge to a unique value: since the contribution from all irrelevant operators vanishes in the continuum. We summarize the differences between the Golterman method and the Kluberg-Stern method in Table XI. We define the chirality operator First let us prove the following theorem.
Proof. Let us first rewrite [γ 5 ⊗ 1] 2 as follows, We know that whereĀ µ = (A µ + 1) mod 2, and Similarly, we find that Hence, we can rewrite Eq. (E5) as follows, where we use the helpful identity η 5 (Ā) = η 5 (A Using the results of Eq. (E4), we can prove the recursion relationship as follows, Using the results of Eq. (E4), we can prove another recursion relationship as follows, Finally, we can prove the following theorem.
Using Eq. (E20), we can prove that for integer n > 0, by induction. At this stage, it will be trivial to prove that The next two theorems concern the chiral Ward identities.
Proof. Using the results of Eq. (E6), we find that whereĀ µ = (A µ + 1) mod 2. Let us rewrite (1 ⊗ ξ 5 ) AB as follows, whereĀ µ = (A µ + 1) mod 2, and Hence, we find that This is a proof of the first part of the theorem. Similarly, This is a proof of the second part of the theorem. (Q.E.D.) We prove the Ward identities in Eqs. (63)-(64) as follows.
Proof. Using the results of Theorem E.3, we know the following Ward identity: Let us begin with the case Q = −2. In Fig. 13, we show leakage patterns of the chirality operator for the first set of zero modes at Q = −2. In Fig. 14, we present the leakage patterns of the shift operator for the first set of zero modes at Q = −2. By comparing Fig. 13 with Fig. 14, we find that the chiral Ward identities of Eqs. (78) and (79) are well-respected. In Fig. 15, we show leakage patterns of the chirality operator for the second set of zero modes at Q = −2. In Fig. 16, we present the leakage patterns of the shift operator for the second set of zero modes at Q = −2. By comparing Fig. 15 with Fig. 16, we find that the chiral Ward identities of Eqs. (78) and (79) are well-preserved. Now let us consider an example with Q = −3. The leakage patterns for the first and second sets of zero modes are similar to those at Q = −2. Hence, we choose the third set of zero modes as our example. In Fig. 17, we show leakage patterns of the chirality operator for the third set of zero modes at Q = −3. In Fig. 18, we present the leakage pattern of the shift operator for the third set of zero modes at Q = −3. By comparing Fig. 17 with Fig. 18, we find that the chiral Ward identities of Eqs. (78) and (79) are well-preserved. Appendix G: Examples for the leakage pattern of non-zero modes Let us begin with an example with Q = 0. Since the gauge configuration with Q = 0 usually has no zero mode (n − = n + = 0), it is relatively easy to study nonzero modes. In Fig. 19, we present leakage patterns of the chirality operator Γ 5 = [γ 5 ⊗ 1] for non-zero modes {λ 1 , λ 3 , λ 5 , λ 7 } = {λ j,m | j = +1, m = 1, 2, 3, 4} in the j = +1 quartet when Q = 0. The results show that the Γ 5 leakages for non-zero modes λ +1,m mostly go into their parity partners {λ 2 , λ 4 , λ 6 , λ 8 } = {λ j,m | j = −1, m = 1, 2, 3, 4} in the j = −1 quartet. Meanwhile, the leakages to other quartets such as j = ±2, ±3 are negligibly small compared to those of the j = −1 quartet elements. This observation is consistent with that for Q = −1 in Fig. 6. In Fig. 20, we present leakage patterns of the shift operator Ξ 5 = [1 ⊗ ξ 5 ] for the non-zero modes {λ 1 , λ 3 , λ 5 , λ 7 } of λ +1,m in the j = +1 quartet when Q = 0. For the Ξ 5 operator, we find the great part of leakages are from nonzero modes λ +1,m to other elements within the j = +1 quartet. Meanwhile, there are only negligible leakages to parity partner quartet elements (j = −1) and other quartets with j = ±2, ±3, and so on. This observation corresponds to the case Q = −1 in Fig. 7. We also find that the leakages of Γ 5 in Fig. 19 and Ξ 5 in Fig. 20 are related to each other by the Ward identity of Eq. (95).