Many-body localization of zero modes

The celebrated Dyson singularity signals the relative delocalization of single-particle wave functions at the zero-energy symmetry point of disordered systems with a chiral symmetry. Here we show that analogous zero modes in interacting quantum systems can fully localize at sufficiently large disorder, but do so less strongly than nonzero modes, as signified by their real-space and Fock-space localization characteristics. We demonstrate this effect in a spin-1 Ising chain, which naturally provides a chiral symmetry in an odd-dimensional Hilbert space, thereby guaranteeing the existence of a many-body zero mode at all disorder strengths. In the localized phase, the bipartite entanglement entropy of the zero mode follows an area law, but is enhanced by a system-size-independent factor of order unity when compared to the nonzero modes. Analytically, this feature can be attributed to a specific zero-mode hybridization pattern on neighboring spins. The zero mode also displays a symmetry-induced even-odd and spin-orientation fragmentation of excitations, characterized by real-space spin correlation functions, which generalizes the sublattice polarization of topological zero modes in noninteracting systems, and holds at any disorder strength.

The celebrated Dyson singularity signals the relative delocalization of single-particle wave functions at the zero-energy symmetry point of disordered systems with a chiral symmetry. Here we show that analogous zero modes in interacting quantum systems can fully localize at sufficiently large disorder, but do so less strongly than nonzero modes, as signified by their real-space and Fock-space localization characteristics. We demonstrate this effect in a spin-1 Ising chain, which naturally provides a chiral symmetry in an odd-dimensional Hilbert space, thereby guaranteeing the existence of a many-body zero mode at all disorder strengths. In the localized phase, the bipartite entanglement entropy of the zero mode follows an area law, but is enhanced by a system-size-independent factor of order unity when compared to the nonzero modes. Analytically, this feature can be attributed to a specific zero-mode hybridization pattern on neighboring spins. The zero mode also displays a symmetry-induced even-odd and spin-orientation fragmentation of excitations, characterized by real-space spin correlation functions, which generalizes the sublattice polarization of topological zero modes in noninteracting systems, and holds at any disorder strength.

I. INTRODUCTION
Complex quantum systems owe their rich physical properties to the intricate interplay of symmetries, disorder, and interactions. This interplay is mirrored by the intertwined concepts and frameworks which have been developed to capture these aspects. Symmetry-reduced representations of quantum systems were established at the beginning of quantum mechanics [1], whilst the absence of unitary symmetries allows for complex wave dynamics even for low numbers of degrees of freedom [2]. This ties both to semiclassical descriptions of classically chaotic systems as well as to statistical descriptions of structureless noninteracting disordered systems, whose properties are captured by random-matrix theory [3]. The latter provides a natural framework to classify complex quantum systems also in accordance with their invariance under antiunitary symmetries, as pioneered by Wigner and Dyson [4,5] and completed with the ten Altland-Zirnbauer universality classes [6]. This ten-fold way also underpins the topological classification of electronic band structures in periodic systems of different spatial dimensionality [7]. Besides time-reversal symmetry, this classification also accounts for antiunitary charge-conjugation symmetries and the combination of the two into a unitary chiral symmetry, as originally encountered in random-matrix descriptions of Dirac systems [8].
For structureless systems, random-matrix theory describes wave functions ergodically spreading across the whole system, but this can be amended to also provide statistical descriptions of Anderson-localized systems, in particular in one-dimensional or quasi-one-dimensional geometries [9]. Again, these descriptions can be organized according to the universality classes of the ten-fold way [10], and then account for topological phenomena in nonperiodic, disordered settings. The most striking effect amongst these is the possibility of such systems to be less localized near spectral symmetry points. This phe-nomenon was first realized by Dyson [11], who noted that a one-dimensional system with hopping disorder develops a logarithmically diverging density of states around the band center-the so-called Dyson singularity, which goes along with anomalously localized states that exhibit a stretched-exponential spatial profile. Within the classification framework described above, this relative delocalization phenomenon becomes tied to the existence of a topologically protected zero mode in a chirally symmetric system with an odd-dimensional Hilbert space [12], and also occurs in higher-dimensional systems, where the anomalously localized states can resemble those at a metal-insulator transition [13,14]. Weaker analogues of such anomalous localization characteristics also occur in absence of spectral symmetries [15,16]. Such robust features deserve attention as they significantly broaden the scope for topologically protected quantum phenomena to realistic, disordered systems, both conceptually as well as in practical terms.
For many-body systems, interactions provide a significant complication of all of these aspects, with much recent effort devoted to the question of ergodic versus many-body localized behavior [17][18][19]. The ergodic phase is again well captured by random-matrix theory [20], reflecting its original setting of nuclear physics [21]. Topological aspects have been extensively studied for gapped ground-state physics [22], but topological order can also emerge in excited states, where it competes with the many-body localized phase [23,24]. In particular, it has been found that such novel phases can be induced by a particle-hole symmetry that pairs up excited states around the centre of the spectrum [24][25][26]. However, the role of many-body zero modes pinned to such spectral symmetry points is much less understood, both conceptually as well as concerning the identification of concrete phenomena.
Here, we clarify this role by drawing motivation from the Dyson singularity. We first identify a natural disordered many-body system displaying an analogous chi-arXiv:2001.07653v2 [cond-mat.str-el] 3 Apr 2020 ral many-body zero mode, consisting of a simple spin-1 Ising chain with a random transverse field. We then address the question of its localization properties both in real space and in Fock space, where we identify two localization phenomena that characterize the zero mode. (i) In real space the spin correlations of the zero mode fragment into 5 independent sectors. This fragmentation occurs both with respect to the even and odd sublattice, which generalizes the sublattice polarization of chiral zero modes in noninteracting systems, as well as with respect to the orientations of the correlated spin, which is specific to the chosen many-body context. This phenomenon holds at all disorder strengths. (ii) In contrast to the noninteracting case, the zero mode still localizes at strong disorder, both in real space as well as in Fock space, as indicated, e.g., by an area law for the entanglement entropy, a large inverse participation ratio, and short-ranged spatial correlations. However, these measures also indicate that quantitatively, the zero mode is noticeably less localized than the nonzero modes-a phenomenon that for brevity we refer to as 'relative delocalization' in the remainder of this work. In particular, the bipartite entanglement entropy is significantly enhanced for the zero mode by a system-size-independent factor of order unity, whilst the inverse participation ratio is correspondingly reduced. Thereby, the zero mode attains characteristics that set it apart from all other states in the system, even if they may be very close in energy.
In Sec. II, we describe the spin-1 Ising chain, which provides a natural model for the described phenomena as it combines a chiral symmetry with an odd-dimensional Hilbert space, and always features a zero mode in one of the two spin-parity sectors. In Sec. III we first discuss the real-space fragmentation of the spin correlations, as this follows directly from the symmetry constraints and holds at all disorder strengths, which we show analytically and illustrate numerically. In this section we also identify the spin correlations that are most characteristic to quantify the localization properties of zero modes and nonzero modes, to which we then turn to in Secs. IV and V. In Sec. IV we demonstrate the relative delocalization of the zero mode numerically based on both Fockspace and real-space measures. The analytical explanation of this relative delocalization is provided in Sec. V, where we identify the dominant hybridization patterns of zero modes and nonzero modes. Enforced by the chiral-symmetry constraints, the dominant zero-mode hybridization configurations involve three spin states on neighboring spins, whilst those of nonzero modes only involve two spin states, so that the Fock-space localization characteristics of these modes fundamentally differ. In Sec. VI we summarize and discuss the results and put them into further context.

A. Model
To demonstrate the effects outlined in this work, we require a many-body system in which a chiral symmetry is manifest for a system with an odd-dimensional Hilbert space. This is naturally provided by a spin-1 Ising chain with a transverse magnetic field, given by the Hamiltonian Here J is the coupling strength between adjacent spins, described by the matrices

B. Symmetries
The primary reason we choose to study the transverse Ising model owes to the fact that it possesses a chiral symmetry with a unitary involution fulfilling X X † = X 2 = 1 1. To see that this chiral symmetry is manifest for all sizes of spin, consider the spin rotation operator with individual rotation matrices U a n (ϕ) = exp (iϕS a n ) .
The operator (5) is well-defined in infinite chains and finite chains with open boundary conditions, while periodic boundary conditions require that the chain consists of an even number of spins, as we will thus assume. The operator X rotates all spins by π about axes alternatingly aligned with x and y, which inverts the sign of all on-site terms ∼ S z n , as well as exactly one of the two spin operators in the interaction terms ∼ S x n S x n+1 , in accordance with the requirements of Eq. (4). The factor i 2SN in the definition (5) makes sure that X 2 = 1 1 holds for all system lengths and arbitrary spin.
The assignment of sites as even or odd provides a gauge freedom, allowing the definition of an alternative chiral symmetry operator It follows that the system also possesses an ordinary symmetry [P, H] = 0, arising from which inverts the sign of all operators S x n in the Hamiltonian and leaves the operators S z n invariant. In the basis (3), so that P represents the total spin parity of the system. Therefore, we can divide the Hilbert space into sectors of even and odd parity, which are of size N + = N − = N /2 for half-integer spins, while for integer spins. Finally, as S x n and S z n can always be represented by real matrices, the Hamiltonian displays a time-reversal symmetry T HT ≡ H * = H, where T is antiunitary and fulfills T 2 = 1. Therefore, the system also possesses a charge-conjugation (particle-hole) symmetry CHC = −H, where C = T X is an antiunitary operator with C 2 = 1.
For spins of size 1/2, it is well known that the rotations defined in Eq. (6) can be written as U a (π) = iσ a with the usual 2 × 2-dimensional Pauli matrices σ a , where one exploits the fact that σ 2 a = 1 1. In the case of spin 1, where the spin operators are given by Eq. (2), we obtain in contrast where In other words, the set of these chiral symmetry operators, with the identity operator, is isomorphic to the Klein four-group, which describes the symmetries of a nonsquare rectangle. We can also exploit the relations [U a , S a ] = 0 and {U a , S b } = 0, provided that a = b.

C. The zero mode
A direct consequence of the chiral symmetry is that eigenstates |ψ k with energies E k are paired with eigenstates |ψk = X |ψ k with energy Ek = −E k . The exception are zero modes with energy E k = 0, for which we formally identify the indicesk = k. Even in the case of degeneracy of these zero modes, they can always be chosen to fulfill X |ψ k = σ|ψ k , where σ = ±1 distinguishes between two types of zero modes. The number ν σ of modes of each type is then constrained by the signature of X , according to tr X = ν + − ν − , which serves as a topological index. In particular, in a system with an odd overall Hilbert space dimension, at least one zero mode is always guaranteed to exist, as it is impossible to pair up all states.
In the Ising chain, the Hilbert space dimension N is even in the case of half-integer spins, and so are the two parity sectors of dimensionality N ± = (2S + 1) N −1 as soon as N ≥ 2. In contrast, according to Eq. (10), for integer spins N − is even but N + is odd. Hence, at least one zero mode, denoted as |ψ 0 , is guaranteed to exist in the even-parity sector of chains with integer spins. In the basis (3), this can be attributed to the existence of the state |0 (the state where s z n = 0 for all spins), which is the only basis state that is left invariant under the operation with the chiral operator X , which connects all other basis states in pairs with index s ands = −s.
The zero mode possesses an energy that remains pinned to zero no matter the disorder or interaction strength. This invariance with respect to parameter variations does not occur for any other eigenstate, which leaves the question if this has any bearings on the localization characteristics, in analogy to what is known from single-particle systems. Therefore, the key question explored in this work is whether the zero mode displays different localization characteristics to the modes with finite energy.

D. Numerical techniques
We will address the localization properties of the zero mode both by analytical and numerical approaches. The numerical results are obtained by exact diagonalization from the positive parity sector in chains with an even number N of spins of size S = 1, where we apply periodic boundary conditions. As the effective Hilbert space dimension rises as (3 N + 1)/2, and only a single zero mode is present in each realization, we obtain disorder averages from chains of limited lengths up to N = 8, but also show results from individual realizations with N = 10. For nonzero modes we collect data from the middle 10% of the spectrum. Quantities assigned to the zero mode are denoted in the form Q 0 , while those of nonzero modes are denoted in the form Q =0 . Disorder averages of any quantity are denoted by an overline, and are obtained from 10000 realizations. Where focusing on individual disorder strengths, we use values W = 1 for weak disorder (ergodic regime), W = 8 for moderate disorder, and W = 20 for strong disorder (localized regime).

III. FRAGMENTATION OF THE ZERO-MODE CORRELATIONS
We start with a general key characteristic of the zero mode, which relates to its real-space structure and holds at all strengths of disorder. We first introduce the spin correlation matrix that captures this structure, and discuss its general properties. We then show that the spin correlations of the zero mode fragment into five independent elementary patterns, whilst for nonzero modes there are only two, and verify and illustrate these patterns numerically.

A. Spin correlation matrix
The real-space spin structure in a given energy eigenstate |ψ k is captured by the correlations which we consider as blocks of a Hermitian matrix C k of dimension 3N × 3N . We term the eigenvalues and eigenvectors of this correlation matrix the correlation eigenvalues and eigenvectors, in distinction to the energy eigenvalues and eigenvectors associated with the Hamiltonian.
The following features are useful to note. (i) According to the relation (S x n ) 2 + (S y n ) 2 + (S z n ) 2 = 21 1, the trace tr C k = 2N is fixed.
(ii) The matrix is well-behaved under local changes (S x n , S y n , S z n ) → (S x n , S y n , S z n )O T n of the spin basis by an orthogonal transformation O n (hence, basis changes that are compatible with the Lie algebra), which transform the correlation matrix as C k,nm → O n C k,nm O T m . This leaves the eigenvalues of C k invariant, while the corresponding eigenvectors automatically adapt to the chosen local spin orientations.
(iii) Since the eigenstates of the Hamiltonian have a fixed parity, the expectation values S z n S x m = S z n S y m = 0. Therefore, the spin correlation matrix decomposes into a direct sum C k = ∆ k ⊕ Z k , given by the block (iv) Utilizing the unitary matrix we can further introduce the transformed matrix with spin-ladder operators S ± n = 2 −1/2 (S x n ± iS y n ). Recalling the analogy between spin-ladder operators and fermionic creation and annihilation operators for systems with spin 1/2, this expression resembles a one-particle density matrix (OPDM), equipped with a Bogoliubov-Nambu structure that is appropriate for a system with a nonconserved particle number.
(v) In a canonical basis state |s , the correlation matrix ∆ s is block-diagonal, with each block having a correlation eigenvalue 1 and an eigenvalue 1−(s z n ) 2 , so that the correlation eigenvectors are localized on individual spins. The correlation matrix Z s then has elements Z s,nm = s z n s z m , and hence is of rank 1, with a single finite eigenvalue Z max s = n (s z n ) 2 (as indicated, we interpret this as the maximal eigenvalue). This counts the number of spins with a nonzero z component.
These features imply that fully localized states are characterized by an approximately quantized correlation spectrum, in close analogy to the OPDM occupation spectrum in a many-body localized system [27][28][29]. In contrast, in an ergodic state represented by a random superposition of basis states, the correlation matrix selfaverages to C erg k ∼ (2/3)1 1. This results in a correlation spectrum centered around the single value 2/3, smoothed out by the influence of the residual off-diagonal elements of C k , which is in close analogy to the smooth OPDM occupation spectrum in an ergodic many-body system.

B. Zero-mode correlations
As we show next, for the zero mode the correlation matrix ∆ 0 further decomposes into four sectors, each pertaining the S x or S y component and additionally confined to the sublattice of even or odd sites. This structure follows directly from the symmetry constraints, and hence holds at all strengths of disorder.
To arrive at these features, we first note that for all states time-reversal symmetry implies as this amounts to an expectation value of a Hermitian operator with imaginary matrix elements, evaluated with a real-valued eigenvector. This constraint does not apply for n = m as the matrix product S y n S x n is not Hermitian (it is furthermore not simply related to S z n , in contrast to the case of spin 1/2). However, for the zero mode the chiral symmetry further implies and analogously which are relations that hold for all n and m. In combination with the constraint (17) from time-reversal symmetry, these relations imply that the blocks ∆ 0,nm are all diagonal, and furthermore vanish if n − m is odd. Thus, for the zero mode the ∆ correlation matrix decomposes into four independent blocks, where the superscripts denote the supporting spin component and sublattice. Including the spin correlations from Z 0 , we can, therefore, identify five independent elementary spin correlation patterns for the zero mode.

C. Numerical illustration
This structure of the spin correlations is illustrated in Fig. 1, where we show correlation eigenstates with minimal and maximal correlation eigenvalues in a typical individual disorder realization at (a) weak, (b) moderate and (c) strong disorder (W = 1, 8, 20, respectively). Subpanels (i-iv) show the four types of correlation eigenvectors from ∆ 0 , while panel (v) shows correlation eigenvectors from Z 0 . The position of these eigenvectors in the occupation spectrum is depicted in the adjacent subpanels (α-γ).
Whilst the predicted fragmented structure holds at all disorder strengths, the correlation eigenvectors from ∆ 0 display a noticeable trend from being extended over the whole system for weak disorder, to becoming highly localized on individual spins at strong disorder. In contrast, we notice that the Z-correlation eigenvectors more sensitively quantify the hybridization of neighboring spins, a feature that will be important in the subsequent sections. In conjunction, the correlation spectra from ∆ 0 and Z 0 both move away from the ergodic value 2/3, approaching the quantized values 1 and 0, respectively, as expected for a many-body localized state.
For comparison, Fig. 2 shows the analogous spin correlation features in a representative nonzero mode. Note that subpanels (i) and (ii) now refer to the x and y components of the same ∆-correlation eigenvectors, as these correlations no longer separate. Furthermore, each of these eigenvectors now populates both the even and odd sublattices. Otherwise, we notice the same qualitative tendencies as for the zero mode-the ∆ spin correlations again become highly localized for strong disorder, whilst the Z correlations remain more extended, and the corresponding correlation spectra move away from their ergodic values 2/3 to quantized values, which now depend on the number of finite spins in the approached basis state |s . In the example, this state has four finite spins, so that there are four nearly-vanishing ∆correlation eigenvalues, and a dominant Z-correlation eigenvalue approaching the value 4.
By surveying different examples we can certify that these qualitative features are typical for individual states in fixed disorder realizations, with the variations at moderate and strong disorder pointing to different spin hybridization patterns. As indicated above, the Z 0eigenvector with maximal eigenvalue Z max 0 is particularly useful to characterize the excitation patterns of the zero mode above the reference state |0 for the zero mode, and analogously above the reference states |s for nonzero modes. These insights will inform our discussion of the quantitative differences in their localization, on which we focus in the following sections.

IV. ZERO-MODE DELOCALIZATION
We now turn to the second key feature of the zero mode, which pertains to the fact that it is less localized than the nonzero modes. In this section, we establish this feature based on numerical results, while the theoretical explanation is provided in the following section.

A. Measures of localization
To address this question, we consider a number of complementary indicators of localization, whose general properties we summarize first.
As a general measure of localization, we consider the bipartite von Neumann entanglement entropy [30]. This is defined for each normalized eigenstate |ψ k as where ρ (k) = tr B |ψ k ψ k | is the reduced density matrix of a subsystem A, obtained by tracing out the complement B. We take A to be a contiguous subchain of length N A = N/2, hence half the length of the total system. In delocalized states, the von Neumann entropy is large, and should be well approximated by Page's law for completely ergodic states [20], S k N A ln 3 − 1 2 . Therefore, the entropy grows linearly with the system size, which manifests a volume law. In contrast, in a localized state the von Neumann entropy is expected to be small, and on average independent of the system size, which manifests an area law. The value of the entropy can then be taken as a proxy for the effective localization length [31]. In a basis state |s , the entanglement entropy S s = 0 vanishes as these states are all separable.
To quantify the degree of Fock-space localization, we make use of the inverse participation ratio (IPR) in the basis (3),     The corresponding correlation eigenvalue spectra are shown in the adjacent panels (α-γ). For weak disorder, these eigenvalues lie around the ergodic value 2/3, whilst for strong disorder they approach quantized values 1 (for ∆) and 0 (for Z).
In the case of perfect Fock-space localization, the IPR goes to unity, whilst in the case of complete delocalization, the IPR goes to 1/N . We also consider the intensity of the states with the special state |0 , which we expect to become large for the zero mode at large disorder, whilst for ergodic states again I k 1/N .

B. Numerical results
At strong disorder, the zero mode is expected to have a large overlap I 0 = | 0|ψ 0 | 2 with the state |0 , in which the contribution from the field h vanishes. This is verified in Fig. 3(a), which shows that the disorder-averaged I 0 rises sharply at disorder strengths W 4. In contrast, the corresponding average I =0 ∼ O(N −1 ) for the nonzero modes is negligible for all disorder strengths. Nonetheless, overall the zero mode is noticeably less localized in Fock space than the nonzero modes, as evidenced in Fig.  3(b) by an inverse participation ratio IPR 0 that is reduced relative to IPR =0 , and in Fig. 3(c) by a bipartite entanglement entropy S 0 that is increased relative to S =0 . Therefore, for strong disorder the nonzero modes approach basis states |s with s = 0 more quickly than the zero mode approaches the basis state |0 . As we explain in Sec. V, the residual hybridization of basis states can be quantified by the maximal Z spin-correlation eigenvalue Z max k , whose average is shown in Fig. 3(d).
In Fig. 4, we show the disorder-averaged entanglement entropy S k as a function of the mode index k, obtained by ordering all states by their energy and centering the Spin correlations for a nonzero mode close to the band center, in analogy to Fig. 1. Note that the correlation eigenvectors displayed in subpanels (i) and (ii) now belong to the same correlation eigenvalues, hence, represent their x and y components, as these no longer separate, also not with respect to the sublattice. Therefore, only two types of elementary spin correlations exist for such nonzero modes. As shown in subpanels (α) and (β), the correlation spectra again become quantized for strong disorder, reflecting the number of spins with finite sz in the approached basis state |s .
resulting index at the zero mode. The entropy of the zero mode is clearly enhanced in the localized regime, by an amount that is independent of the accessible system sizes, hence remaining consistent with an area law. This well-confined relative delocalization peak also confirms that the enhancement is restricted to exact zero modes, and not shared, e.g., by nonzero modes very close to the band center. As shown by the statistical distribution functions of the entanglement entropy in Fig. 5, this delocalizing tendency can be attributed to an accumulation of zero modes with entropy S 0 slightly above 1, to be identified as S 0 ln 3 in the following section. This accumulation is already well pronounced at moderate values of disorder (panel b), and is well defined at very large values of disorder (panel c), suggesting that it arises from a specific delocalization mechanism. In contrast, nonzero modes display accumulations at smaller characteristic values of the entropy, to be identified as ln 2 and (ln 8)/2, which hints towards a competition of several distinct delocalization mechanisms. We will identify the underlying hybridization patterns in the following section.

V. DIMER HYBRIDIZATION
We explain the relative delocalization of the zero mode based on quasi-degenerate perturbation theory at relatively large disorder. This reveals a characteristic dimer hybridization pattern involving three collective ba-sis states localized on neighboring spins, whilst nonzero modes support a much wider range of hybridization patterns. These characteristics define the typical features of all modes in the strongly localized regime, where any hybridization is absent. The question is how the modes gradually delocalize due to resonant interactions at weaker disorder. We show that this involves distinct hybridization processes on adjacent spins, leading to characteristic features in the entropy, inverse participation ratio, and spin correlations.
We first identify the resonance conditions in general terms, and then derive the hybridization patterns and their characteristic signatures, which we further support with numerical results. , which quantifies the residual hybridization of these states in the strongly localized regime, as further discussed in Sec. V.

B. Resonance conditions
In first-order perturbation theory, the hybridization of the zero mode with other states |s is strongly suppressed by energy denominators E (0) s . In principle, hybridization can set in for states with individual |h n | ≤ J, for which individual spins can align freely. However, at least two sites need to be involved to retain positive parity, and furthermore these configurations have vanishing perturbation matrix elements unless sites neighbor each other. On the other hand, it should then suffice that |h n | − |h n+1 | J, instead of both |h n |, |h n+1 | J individually, implying that such disorder configurations should be dominant as they require fewer constraints.
We can verify the above reasoning by examining all excitations patterns above the background state |0 .
Amongst the excitations involving neighboring spins (hence relevant in the first order of the perturbation), only two patterns are allowed by parity, chiral, and timereversal symmetry, namely those obtained from state |0 by terms generated via application of the matrix combinations iS x n S y n+1 and iS y n S x n+1 . These can be conveniently combined into excitation operatorŝ leading to the perturbative ansatz where φ ± n are the amplitudes of the two excitation fields. Expanding the condition H|ψ 0 = 0 in orders of the relative interaction strength, we then obtain the perturbatively closed coupled equations whereupon Thus, one of the two fields becomes large when |h n | |h n+1 |, which agrees with the resonance conditions identified above.

C. Zero-mode hybridization patterns
To describe such resonant disorder configurations more accurately, we resort to quasi-degenerate perturbation theory in the subspace of the hybridizing spins, taken without loss of generality as (n, n + 1) = (1, 2). We start with dimer hybridizations of even parity, assuming initially that they are embedded into a chain where the remaining spins are unhybridized, We first consider the vicinity of the resonance condition (30), where we write h 1 =h + δ/2, h 2 =h − δ/2 whilst setting J = 1. Ordering the even-parity states as |1, 1 , |1, −1 , |0, 0 , |−1, 1 , |−1, −1 , the reduced Hamiltonian then separates into three sectors, with states |1, 1 and |−1, −1 gapped out by an energy ±2h, whilst the zero mode is contained in the quasi-degenerate sector spanned by the states |1, −1 , |0, 0 , |−1, 1 . Diagonalizing this sector, we find two states of finite energy ± δ 2 + 1/2, to which we will come back later, as well as a zero mode of vanishing energy, which we will call the dimer zero mode. Near the resonance condition (31), the same considerations apply upon writing h 1 =h + δ/2, h 2 = −h + δ/2 with the roles of the states (|1, 1 , |−1, −1 ) and (|1, −1 , |−1, 1 ) interchanged, leading to zero-mode hybridizations In both cases, the bipartite entanglement entropy of the dimer zero mode is given by and the inverse participation ratio is given by On the dimer, the Z correlation matrix has a single finite eigenvalue which we interpret as the maximal eigenvalue as for the remaining spins Z 0,nn = 0 vanishes [32].
The top row in Fig. 6 displays these characteristics of the dimer zero mode as a function of the detuning δ. The entropy has a stationary point at δ = 0 with the value S 0 = ln 2, where IPR 0 = 1/2 and Z max 0 = 2, and two stationary points at δ = ±1/2 with the value S = ln 3, where IPR 0 = 1/3 and Z max 0 = 4/3. We note that several of these hybridization patterns can be embedded along different positions of the zero mode. The entropy then arrives from the dimers spanning the bipartite partition point, and still adheres to Eq. (38). The resulting inverse participation ratio is the product of those of all hybridized dimers, so that Eq. (39) provides an upper bound for IPR 0 . Furthermore, the Z correlation matrix decomposes into independent blocks, so that Eq. (40) provides a lower bound for the maximal Z correlation eigenvalue Z max 0 .
Overall, we therefore arrive at seven hybridization patterns and two nonhybridized states, reflecting the full di-   Fig. 7(a,b). Note the different scale on the horizontal axis in panel (a).
mensionality of the dimer subspace. For the second resonance case h 1 =h + δ/2, h 2 = −h + δ/2, the same considerations apply upon interchanging dimer basis states |s, s ↔ |s, −s . The characteristic features of these finite-energy hybridization patterns are shown in the bottom row of Fig. 6. Hybridizations based on |ψ 0 dimer still produce the same entropies and inverse participation ratios as for the zero mode, while the largest eigenvalue of the Z correlation matrix now arises from the remainder of the chain, where it counts the number of finite spins, thus giving rise to the straight lines at even integers. For the hybridizations |ψ +,± dimer of even parity, the entropy is stationary around δ = 0, where S k = (ln 8)/2 whilst the inverse participation ratio takes the value IPR k = 3/8. For the hybridizations |ψ −,±,k dimer of odd parity, a similar behavior is observed with stationary entropies S k = ln 2 and inverse participation ratios IPR k = 1/2. In both these hybridization patterns, the eigenvalue Z max k depends both on the hybridization strength δ and the number of finite spins in the remainder of the chain. Considering that several hybridized dimers can occur along the chain, we again interpret the predicted inverse participation ratios and Z correlation eigenvalues as upper bounds and lower bounds, respectively.

E. Summary and numerical verification
Summarizing the results from this section, we arrive at the following detailed predictions.
For the zero mode, delocalization occurs via dimer hybridization patterns with typical entropies S 0 ∼ ln 3 or S 0 ∼ ln 2, as already observed numerically in the upper panels of Fig. 5. Entropies S 0 ∼ ln 3 are further expected to correlate with inverse participation ratios bounded as IPR 0 1/3 and leading Z correlation eigenvalues bounded by Z max 0 4/3, while for entropies S 0 ∼ ln 2 we expect IPR 0 2 and Z max 0 ∼ 2. More generally, these quantities should be correlated as shown in the upper panels of Fig. 7. These predictions are verified in Fig. 8, where we show scatter plots of the described quantities from 10 4 realizations for chains of length N = 8.
The expected correlations are already well established for moderate values of disorder W = 8.
Furthermore, nonzero modes should predominantly display entropies around S k = ln 2, which can be achieved by the widest variety of hybridization patterns, followed by S k = (ln 8)/2, whilst S k = ln 3 should occur relatively less frequently, as indeed observed in the lower panels of Fig. 5. The expected correlations with IPR k and Z max k are depicted in the lower panels of Fig. 7. These predictions are verified in Fig. 9.
Comparing the results in Figs. 8 and 9, we find that the zero modes and nonzero modes are most clearly discriminated by their distinct correlations between the inverse participation ratio IPR k and the leading Z-correlation eigenvalue Z max k .
Having confirmed these key predictions, we return to Fig. 6 to observe that the dominant hybridization patterns of the zero modes are appreciable over a larger range of detunings δ than for the nonzero mode. This verifies that the zero mode hybridizes more readily than the nonzero modes, and then exhibits more delocalized Fock-space configurations, which provides the general explanation for the numerical observation of this effect in the previous section.

VI. DISCUSSION AND CONCLUSIONS
In summary, in many-body systems zero modes protected by a chiral symmetry can localize, but then do so with distinctively different characteristics than nonzero modes. In particular, the zero modes are more delocalized both in terms of their real-space and Fock-space signatures. We explained these differences by the characteristic symmetry-restricted mechanisms allowing the localized basis states to hybridize. These symmetry constraints can be extended to all disorder strengths by considering the fragmentation of real-space correlations.
We developed and demonstrated these effects for the example of a disordered spin-1 Ising chain. For spin 1/2-chains, the chiral symmetry is already present, but symmetry-protected zero modes do not occur as the Hilbert space dimension is even in both parity sectors. In spin-1/2 chains, the nonzero modes are known to delocalize by a single dominant hybridization pattern, involving dimers with bipartite entanglement entropy S k = ln 2 [33]. In contrast, in the spin-1 chain, the delocalization mechanism of zero modes involves a dominant hybridization pattern with entropy S 0 = ln 3, whilst nonzero modes involve a competition of various hybridization patterns, including such with entropy S k = (ln 8)/2. Even though the underlying hybridizations differ, these entanglement values are reminiscent to those encountered in the fragmented ground state of the spin-1 system by Affleck, Kennedy, Lieb, and Tasaki (AKLT model) [34], as well as for other spin-1 systems, where such entanglement entropy values can be found between a single spin and the remainder of the system [35,36]. Furthermore, for the studied system the fragmentation of real-space correlations occurs both with respect to the spin orientation as well as with respect to the even and odd sublattices, where the latter is particularly noteworthy as statistically the system is translationally invariant.
A fundamental tenet for disordered interacting quantum systems is the expectation that many-body states close in energy share the same statistical signatures. The symmetry-protected zero-modes discussed here provide a mechanism to equip individual states with their own characteristic signatures. It would be interesting to explore whether the remarkable differences between zero modes and nonzero modes become further accentuated for larger integer spins, and whether these observations also extend to appropriately designed itinerant fermionic systems. the spectrum (note that the chiral symmetry furthermore implies r k = r −k ). In Fig. 10, the disorder-averaged spacing ratios are shown as a function of disorder strength for three system sizes N = 6, 8, 10. The statistical fluctuations for r 1 are large as only a single value is obtained for each realization. Nonetheless, both figures consistently point towards states becoming localized at about the same strength of disorder, with the averaged ratios of different system size crossing near a point of inflection at around W 5.
In Fig. 11, we show the full statistical distribution of the spacing ratios for the smallest system size N = 6, where enough data can be collected, for representative values of the disorder strength W = 1, 5, 20. The results for r 1 and r >1 resemble each other closely in all three cases, being consistent with ergodic behavior for W = 1 as well as many-body localized behavior for W = 20, and displaying similar intermediate statistics for W = 5.