Edge mode locality in perturbed symmetry protected topological order

Spin chains with symmetry-protected edge zero modes can be seen as prototypical systems for exploring topological signatures in quantum systems. These are useful for robustly encoding quantum information. However in an experimental realization of such a system, spurious interactions may cause the edge zero modes to delocalize. To stabilize against this influence beyond simply increasing the bulk gap, it has been proposed to harness suitable notions of disorder. Equipped with numerical tools for constructing locally conserved operators that we introduce, we comprehensively explore the interplay of local interactions and disorder on localized edge modes in the XZX cluster Hamiltonian. This puts us in a position to challenge the narrative that disorder necessarily stabilizes topological order. Contrary to heuristic reasoning, we find that disorder has no effect on the edge modes in the Anderson localized regime. Moreover, disorder helps localize only a subset of edge modes in the many-body interacting regime. We identify one edge mode operator that behaves as if subjected to a non-interacting perturbation, i.e., shows no disorder dependence. This implies that in finite systems, edge mode operators effectively delocalize at distinct interaction strengths. In essence, our findings suggest that the ability to identify and control the best localized edge mode trumps any gains from introducing disorder.


Introduction
Topological states of matter have been the focus of intense research over the past 30 years. Within systems of condensed matter physics, topological effects are known to occur in quantum Hall systems of the electron gas [1] and topological insulators [2]. Experiments on wires with proximity induced superconductors gave compelling evidence for Majorana zero modes [3]. Cold atomic gases and photonic devices offer possibilities of creating synthetic topological properties [4]. These new phases of matter by definition have no local order parameter but can be detected via their entanglement properties [5] and are classified by topological invariants [6]. When considering one-dimensional spin chains, these invariants give rise to protected gapless edge modes [7][8][9][10][11][12][13][14], which survive only if perturbations do not break the symmetries of the Hamiltonian.
Such edge modes are interesting from the perspective of quantum information science as well: They are one of many proposed candidates to encode quantum information robustly using topology [6,9,15]. However, the localization of the edge modes can be compromised by the onset of interactions allowing them to hybridize via transport. To counteract this effect, it has been suggested that topological phases can be stabilized by disorder [16][17][18], which is supposed to inhibit transport by localizing the bulk. These works have given rise to the narrative that disorder is expected to always be beneficial when it comes to protecting topological signatures.
The interplay of topological features, interactions and disorder is far from being fully understood; the narrative of stabilization by disorder is at present not backed up by comprehensive studies, numerical or otherwise. What is known rigorously is that, for topologically ordered systems, sufficiently weak local perturbations do not lift the ground-state degeneracy [19][20][21] -but this kind of statement rather shows that small noise levels are not detrimental than making explicit constructive use of them. While this implicitly defines a coherence length for the edge modes, it is far from clear how the local structure of these operators is deformed in the presence of interactions and disorder. These question, though basic, must be addressed before more sophisticated scenarios can be meaningfully studied. This work sets out to tackle the issue of the deformation of edge modes under disorder in a comprehensive fashion, laying the groundwork for a more general picture of the interplay of topological features and disorder.
In the following text, we will study the XZX cluster Hamiltonian, a topological chain which hosts one qubit at each edge and analyze if disorder can help stabilize them [17,18] in the presence of weak interactions. We devise algorithms capable of calculating the support of the edge operators of the disordered XZX cluster Hamiltonian perturbed by either XX or XXZ type interactions. Equipped with these tools, we are in a position to determine the sensitivity of the edge mode localization length to each perturbation type. Contrary to the established narrative, we find that disorder only aids localization slightly, and only in the presence of interaction terms which are non-quadratic in the fermionic dual. Furthermore (and surprisingly), we also find that some edge modes are completely insensitive to disorder. Building on these findings, we elaborate on the lessons to be learned on the interplay of disorder and topological features.

SPT chains with spatial disorder
Our main focus is on the interplay of disorder, interactions and SPT order. Take for example a spin chain hosting a XZX cluster Hamiltonian where the h j are drawn uniformly from − ∆ 2 , ∆ 2 and X j , Y j , Z j are the Pauli operators acting at site j. This system is known to be in a symmetry protected topological (SPT) phase, and thus supports localized, spin 1/2, edge zero modes. The choice of disorder model here may seem unphysical to readers familiar with many-body localization, where a disordered local magnetic field is commonplace. Such a field competes with the SPT order, driving a transition to a topologically trivial phase. By disordering the cluster terms directly, we are implementing an ideal version of disorder, in that it breaks the degeneracies in the excited state spectrum while preserving the ground state manifold. Thus, if we should fail to observe increased localization in this circumstance we do not expect any improvement by moving to a more physical model of disorder.
Without any extra interaction terms, the two edge zero modes enforce a four-fold degeneracy at every energy level in the spectrum, and are perfectly localized on the two sites nearest to the boundary. By inspection, we can find local operators which exactly commute with the Hamiltonian (1) that are located at the left and right edge. Apart from these local conserved quantities, the Hamiltonian in (1) also commutes with the time reversal operator where K is complex conjugation. Note that all local edge operators fail to commute with time reversal and thus cannot be used to split the degeneracy without breaking the symmetry. Each set of these operators describe a spin-1/2 Hilbert space. Although T 2 = 1 on the full Hilbert space, T 2 = −1 when restricted to these spin-1/2 Hilbert spaces, i.e. these edge states transform projectively under the global time reversal (see Appendix A).
To perturb H 0 , we will introduce a translationally invariant XXZ-coupling to our spin chain of the form which, for J 1, is representative of interactions typically found in real solid-state material where Heisenberg-type interactions are ubiquitously present as a result of exchange interactions. Since we will only be investigating the effects of either an XX or a Heisenberg perturbation, we have included a parameter η which interpolates between the two, and leave J as the overall interaction strength. Specifically, we consider the following Hamiltonian defined on N lattice sites where we choose η = 0 or η = 1. By choosing ∆ = 0 we can switch on the presence of local disorder that can have the effect of diminishing the influence of the perturbation added to the exact Hamiltonian.
Note that H int (J, η) commutes with the time reversal operator for any values of J and η, so if it is sufficiently weak it will only lift the degeneracy by an amount exponentially suppressed in system size [7][8][9]23]. This occurs because, as soon as J = 0, the edge mode will no longer be perfectly localized at the edges and are in fact expected to be smeared within an exponential envelope [7][8][9][10][11][12][13][14]. With the degeneracy lifted, the edge mode operators will no longer commute exactly with the Hamiltonian, since the existence of operators which anticommute with T (a feature of the edge modes in (2)) and commute with Hamiltonian would require exact degeneracies due to Kramer's theorem. The failure of the edge mode operators to commute exactly with play an important role in informing our algorithm in Section 3.2.
From here, we set out to understand this interplay of topology, interactions and disorder by explicitly constructing edge modes for this perturbed XZX cluster Hamiltonian. For Anderson and many-body localized systems, the localization length of all local conserved quantities depends on the disorder strength. Thus, one would expect that localizing the bulk of the SPT chain should stabilize the edge modes as excitations cannot traverse the full system to allow the hybridization of opposite edges [17,18]. We devise methods capable of computing the edge mode support in presence of both non-interacting and many-body interacting perturbations. Contrary to the previously stated heuristic argument, we find numerous cases where the edge modes are completely insensitive to disorder.

Edge mode construction
In this section, we describe two methods employed to construct the edge mode operators E in the perturbed XZX cluster Hamiltonian. Computing the broadening of the edge modes is a particularly daunting task, precisely because the operator encodes information about states throughout the entire spectrum, and thus cannot be studied using low energy techniques such as DMRG. The first method assumes that the perturbed Hamiltonian represents non-interacting fermions and yields an efficient solution in terms of Majorana eigenmodes, which allows for a direct computation of the edge modes. The second approach tackles the generic interacting case, where no Bogoliubov transformation will suffice and hence the construction of the edge modes becomes more intricate. In this case we rely on a method developed to construct conserved operators called "l-bits" for a many-body localized system [24,25]. The construction of these conserved quantities from first principles is difficult, but algorithms which can construct them using various methods do exist [26][27][28][29][30][31][32][33][34][35][36][37][38].

Edge modes under free fermion perturbations
After a Jordan-Wigner transformation, the choice of η = 0 is equivalent to a non-interacting fermionic problem whereas η = 0 maps to an interacting fermionic model with quartic interactions. Introducing the Majorana operatorsγ j , γ j for j = 1, . . . , N as with the Hamiltonian (5) becomes which is non-interacting if and only if η = 0. Written in terms of the Majorana operators, the edge modes for J = 0 take the form with P = Z 1 · · · Z N being the global parity operator which commutes with H. For this we note that (8) can be written as with the coupling matrix As C ∈ R N ×N is real, the singular value decomposition of C takes the specialized form C = Q T ΣQ with two orthogonal Q,Q ∈ O(N ) and Σ ∈ R N ×N a diagonal matrix with real non-negative entries. Using the two orthogonal matrices Q,Q ∈ O(N ) we introduce new modes which again fulfill the Majorana anti-commutation relations (7) and the Hamiltonian (10) becomes diagonal taking the form where we have defined the single particle energies σ j = Σ j,j and assume without loss of generality that they are in increasing order. For J = 0, we find that σ 1 = σ 2 = 0 with m 1 = γ 1 , m 2 = γ 2 ,m 1 =γ n ,m 2 =γ n−1 being the corresponding localized edge mode operators. At finite J > 0, σ 1 ∼ σ 2 ∼ e −n/n 0 are not exactly zero anymore but decay exponentially with increasing system size [9] and hence much smaller compared to the next largest value σ 3 . Thus an approximate four-fold degenerate ground-state sector remains well defined. The operators m 1 , m 2 ,m 1 ,m 2 hence correspond to perturbed edge mode operators which can be individually studied in the free fermionic setting via their single particle wavefunctions. Note however, that only their products m 1m1 and m 2m2 , which are supported at both ends of the chain, are exact constants of motions of the Hamiltonian. This will also be a main difference to the interacting relaxation algorithm which from the outset seeks operators that exactly commute with the Hamiltonian.

Edge modes under perturbative many-body interactions
The intuition behind our approach to constructing edge modes of a system with many-body interactions is as follows: the edge modes E 0 of the unperturbed model H 0 should deform smoothly to those of the full interacting Hamiltonian. Indeed, they turn out to be good starting points to obtain the actual edge mode operator E which commutes with H exponentially well in the system size whilst remaining local to some degree.
The method requires an ansatz which is expected to resemble the conserved operators. Since we are perturbing away from a solvable point, we employ a natural choice, the exact edge modes obtained in the unperturbed fixed-point model. One might be inclined to use a single edge mode as an ansatz for finding the perturbed operators, but this approach will fail in general as the operators produced with this method necessarily commute exactly with the Hamiltonian by construction. This is not the case for single edge modes as they cannot commute with Hamiltonian, as discussed in Section 2. We can circumvent this problem by instead using products of edge modes supported on both the left and right ends of the chain which we call B 0 = E L 0 ⊗ E R 0 . Such products respect the time reversal symmetry and thus are not prevented from commuting exactly with the Hamiltonian. Due to the topological degeneracies present in our model, we have to make sure that the basis in any subspace also diagonalizes our edge mode guess E 0 if we want to obtain the form in Eq. (15). This is reminiscent of standard degenerate perturbation theory and in fact requires by far the most ressources of the total algorithm as we need to rediagonalize 2 N −2 many 4 × 4 matrices.
We will now detail how to construct the quasi-local conserved operators. By definition, these have a compact representation in the energy basis. This basis, which we label by {|k }, is obtained from full exact diagonalization. Consider a basis for diagonal operators in energy space fulfilling the Pauli-algebra. The minimal elements of this basis may take the following form where Z is a Pauli-operator. This might look complicated at first glance, but it is really nothing but Pauli-operators defined in energy space. The full basis can be obtained by calculating all products of these N -many operators. These operators by construction exactly commute with the Hamiltonian. Hence, any constant of motion can be brought into the form of the Ξ operators. Their real space representation U D Ξ i U † D will however in general not be local. Their locality is completely dependent on the ordering of the eigenstates in the unitary U D as this is the only freedom left. Due to the immense number of possible permutations -the order of the symmetric group S 2 N is 2 N ! -a brute force approach is out of scope. We instead rely on a heuristic algorithm which dynamically relaxes an ansatz operator to obtain a good permutation. Abstractly speaking, we bank on the time independent or equilibrium part of B 0 to resemble the product of the perturbed edge modes B already quite well. In cases where this is not given, the algorithm will fail to produce a local edge mode.
The unperturbed operator B 0 and a diagonalization unitary U D serve as inputs to our method. This unitary has an arbitrary ordering of eigenvectors at the start (for ED, usually determined by the size of the energies of the Hamiltonian). Upon mapping B 0 to its equilibrium representation, we use it to obtain an ordering of the eigenstates that resembles the Pauli structure well. This representation is obtained by calculating the infinite time average of B 0 , which stems from equilibration theory the time average will hence be diagonal in the energy eigenbasis for non-degenerate spectra and is thus a constant of motion. However, since the infinite time average is not trace preserving, it in general causes B 0 to lose its algebraic structure. We would like to point out that while localizing systems are in general not expected to thermalize, they do equilibrate which makes this ansatz meaningful [39].
Due to the topological degeneracies present in our model, we have to make sure that the basis in any subspace also diagonalizes our edge mode guess E 0 if we want to obtain the form in Eq. (15). This is reminiscent of standard degenerate perturbation theory and in fact requires by far the most resources of the total algorithm as we need to rediagonalize 2 N −2 many 4 × 4 matrices.
We then set out to find a permutation of the eigenvectors of H such that the time average of the B 0 best resembles Ξ 1 which can be done by a sorting of the eigenvalues of E(B 0 ). This permutation P also gives rise to a new diagonalization unitary U D = PU D . Upon conjugating Ξ 1 with U D , we obtain an edge mode which fulfills all algebraic properties and commutes with the Hamiltonian. Since our sorting method is heuristic, we cannot rule out the existence of better localized edge modes. Nevertheless, the support that we find serves as a robust upper bound. Because of this, we note that a breakdown of our method, i.e. finding a non-local operator, does not necessarily imply that there are no localized edge modes.
The following pseudocode describes a possible way to implement this procedure numerically. We use a notation close to python. This algorithm builds on a previously introduced method used in the context of many-body localized systems [38]. In this problem, the authors designed operators which commute exactly with the given Hamiltonian and are quasi-local. In a many-body localized system, one searches for extensively many quasi-local constants of motion and the system features a fully non-degenerate spectrum caused by the disordered potential landscape. In contrast, the SPT model is characterized by only constantly many edge mode operators which enforce degeneracies throughout the spectrum. These differences necessitated major modifications to the method from Ref. [38].

Measure of locality
In the following analysis we set out to assess the locality of the constructed edge mode operators. Therefore, we want to compare the action of the full operator to a itself truncated to a local region only. As a first step, we will need to specify a reduction map, which reduces our operators to such operators with local support in a region S where support is defined using site indices. This map truncates an operator down to its local support on S and afterwards embeds it into the full real space again by tensoring identities on S C . This operator can now be compared to the original operator supported on the full system. The difference between the two will be a measure of the support Due to the interacting procedure yielding products of edge modes, we expect their support to be mainly on both edges of the system. To assess their locality, we hence use an S which is centered in the middle of the chain and extends by increasing this block on its both ends by one site. We note that the norm used here is most sensitiveand in many other applications operators which are expected to be local are so only in weaker norms than in operator norm [37,38,40,41].
In the non-interacting case η = 0 we have access to the individual edge mode operators m 1 , m 2 , m 1 ,m 2 and we hence consider the sets S L,k = [k] and S R,k = [N − k] C oriented at the left and Here, color encodes the three modes and markers again encode disorder strength ∆. The data for Z L X R and Z R X L overlaps completely which is why it is hard to spot the orange markers indicated in the legend. The grey line is taken from the non-interacting results as a comparison. right boundary of the system. The larger system size considered in this case, prohibits to use the full Hilbert-space representation of the operators. However, as we show in [22], one can exploit the algebraic properties of the non-interacting fermions in order to directly compute the reduction and norm of it within the fermionic picture which yields for p = 1, 2

Numerical Results
In this section, we show and discuss the resulting operators for both models. We have worked with at least four interaction strengths J ∈ {10 −2 , 10 −3 , 10 −4 , 10 −5 } and three disorder strengths ∆ ∈ {0.1, 0.3, 0.5}. For a clearer presentation, we only picked a subset of these results but the calculations not shown behave analogously.

Free fermionic perturbation
For η = 0, as described above, we show the support of the single edge modes supported on the left part of the chain. The system has a total size of N = 32 sites. While much larger systems are treatable and have been investigated by us with this algorithm, we find that this system size suffices to properly display the edge mode decay. For a system size scaling of this method, we refer interested reader to the Appendix D. Results can be found in Fig. 1. Left and center plots show the logarithmic support log 10 supp(E, S) of the edge modes m 1 , m 2 . We use the data of these plots to extract a localization length ξ, shown in the right plot. The errorbars for small interactions stem from the fact that in these systems, the support of the edge mode falls of very strongly yielding only few non-zero points and therefore less accurate fits. We find that the inverse localization length ξ depends logarithmically on the interaction strength J. The different modes m 1 and m 2 show the same qualitative behaviour. The support falls off in exponential fashion with the size of the support region. This aligns nicely with the intuition that additional interaction terms should only dress the original modes. Furthermore, the observed plateaus can be derived for the infinite system size limit as shown in Appendix C. With increasing interaction strength J, the edge modes become less local, as expected from perturbation theory.
A feature of particular note in these results is the insensitivity of the edge mode locality to the disorder strength ∆. As a comparison, we also show data without any disorder. This is surprising when contrasted with the intuition that disorder should help localize the edge modes [17,18]. This suggests that the edge modes of this SPT do not couple to the bulk operators in circumstances of Anderson localization.

Many-body interacting perturbation
Now we resort to the calculations performed for η = 1, corresponding to the interacting system. Due to the size of the Hilbert space and the effort of the re-diagonalization of the topologically degenerate subspaces, we had to resort to system size N = 12. Fig. 2 again shows the logarithmic support log 10 supp(E, S) for E ∈ {Y L Y R , X L Z R } on the left and center panel. The right panel shows the extracted localization length. A more detailed discussion on the fitting procedure and cross validation of the code can be found in Appendix D. The symmetry of the plot is due to the choice of the system S as laid out in section 3.3.
The support again shows an exponential decay of operator support into the bulk. Moreover, this localization length ξ grows with increasing the interaction strength J as expected. The localization length observed for all edge modes is of the same order as the one found in much larger system sizes for the non-interacting perturbation (cp. gray dash-dotted line). However, when increasing the interaction strength to J = 10 −2 there is a sharp drop in the localization length for some operators which might be ascribed to a transition towards a topologically trivial phase. This transition point is far lower than the expected value of J ∼ 1. This is an expected finite size effect as the edge modes are a priori closer together and thus able to hybridize more easily. Nevertheless, the compatibility between non-interacting and many-body interacting localization lengths away from this transition indicates that the signal of SPT behaviour can still be reliably observed in system sizes tractable by exact diagonalization.
For the Z L X R mode we find that the heuristic picture is recovered as increasing disorder strength aids localization. This contrasts strongly with our findings in the non-interacting case, indicating that many-body interactions are necessary to couple the edge modes to the bulk operators. Moreover, the localization length is generically longer than in the free fermion case with disorder strength pushing the length down towards the free fermion value. This suggests that the free fermion value represents the best localization of the edge modes for fixed interaction strength. The same behaviour is found for X L Z R (cp. [22]).
An exception to the behaviour reported above is displayed in the localization length of the Y L Y R edge mode. Despite the presence of many-body interactions, it shows a disorder insensitivity akin to that of the non-interacting regime. This goes beyond mere analogy as the value of the localization length of the Y L Y R operator overlaps perfectly with the non-interacting results. We computed the localization behaviour for all six possible edge mode hybridization patterns and found that only the Y L Y R operator displays free fermion localization behaviour. This suggests that this mode is subject to a selection rule which precludes the many-body interaction effects which delocalize the Z L X R mode. The source of this selection rule is at this point mysterious, but we note that the Y L Y R operator is unique among the choices of edge mode products in being local to the edges in both spin and fermionic variables, i.e. it does not feature a parity string across the whole chain. The absence of such a non-local feature in the fermionic picture may explain the reduced sensitivity to bulk localization behaviour.
Put succinctly, our results suggest that in the presence of many-body interactions, there may be a splitting of the modes into those which delocalize faster, i.e. Z L X R and X L Z R and are sensitive to disorder and one mode Y L Y R , which is insensitive to disorder and shows a stronger localization comparable to the one of non-interacting edge modes. We would like to point out, that since our method can only provide upper bounds to the localization behaviour, it is still conceivable that all three modes behave the same. Also, it is possible that the disorder sensitivity observed in all other products vanishes in larger systems than we are able to treat. However even if a finite size effect, this splitting constitutes an interesting result as it would be relevant for short synthetic chains or cold ion systems. In such circumstances where one seeks to improve edge mode locality in presence of manybody interactions to encode quantum information, the gains from disorder potentials are marginal compared to those from picking "better" edge modes.

Conclusions
In this work, we investigated the localization behaviour of topological edge mode operators upon introducing both non-interacting and many-body interacting perturbations as well as disorder. Specifically, we started out from the disordered XZX cluster Hamiltonian which as a fixed-point model is exactly soluble and added XX and XXZ interactions which are expected to drive the transition towards a topologically trivial model. We introduce different methods of finding the topological edge mode operators, one based on the Majorana description which yields the lowest lying eigenmodes for non-interacting systems and a second one, which uses the relaxation of the fixed-point edge modes as an ansatz to heuristically find local edge modes for many-body interacting chains. While the support of the obtained edge operators with the interacting method is only an upper bound, the commutation with the Hamiltonian is exact.
Both interactions considered delocalize as the interaction strength is increased. However, the noninteracting model displays no disorder dependence whereas the interacting system does. Curiously, a single edge mode combination which in the fermionic language corresponds to the two density operators at both ends, namely Y L Y R , shows no disorder dependence even when adding many-body interactions. Our results suggest that for a finite size chain, one might find different localization behaviour for different edge mode operators. Specifically, we find one edge mode that is most stable and completely insensitive to disorder, picking it out as the one best-suited to encode a logical qubit.
Since we fully diagonalize the Hamiltonian we are limited to small system sizes even for this onedimensional problem. We hope to extend the method to larger systems by truncating to the ground state sector, which would possibly allow a tensor network implementation as well. The interacting method used in this work relies only on guessing a suitable ansatz edge mode operator. Hence, we plan to apply it to more physical models and other types of perturbations such as open dynamics. Here, one might hope to overcome the thermal instability of topological systems [42] with the help of disorder [43].

Acknowledgements
This work has been supported by the ERC (TAQ), the DFG (CRC 183, FOR 2724, EI 519/7-1), and the Templeton Foundation. This work has also received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 817482 (PASQuanS).

A Fractionalization
To see this, we must rewrite the time-reversal operator in a way that makes the edge action explicit. We can do this by re-expressing our time reversal operator using the cluster operators of our Hamiltonian, which lets us recast T as If we decompose our Hilbert space into edge and bulk tensor factors, we can identify the emergent edge action and so we can see that the localized operators of time reversal on the edge states are which, curiously, these operators do not square to 1. Instead, which demonstrates the symmetry fractionalization expected in an SPT phase.
B Operator norm decay of edge modes for η = 0 In the following we explicitly compute the norm of the edge mode operators m p and Pm p with p = 1, 2 obtained in the case η = 0. Here, it is important to keep in mind, that we want to study the effect of the reduction map on the level of the qubits in which the original Hamiltonian in (1) is defined. All following norms and traces are hence evaluated by reversing the Jordan-Wigner transformation and relating the fermionic operators to the Pauli operators. The results are then mapped again to the fermionic level as the expressions take a more compact form here.
Due to the non-interacting structure of the problem, the edge mode operators are linear combinations of the initial Majorana operators γ j ,γ j . In this case the result of the reduction map can be studied in detail. We find where, as explained above, the trace is evaluated on the level of the qubits. The differences A − Γ(A) are then given by and again essentially only linear combinations of γ j andγ j . We can however compute easily the norm of any linear combination of Majorana operators. Let S ⊂ [N ] and define A = j∈S a j γ j (30) to be any linear combination of the γ j operators with a j ∈ R for j ∈ S. One finds that the square of A is given by a j a k γ j γ k = j∈S a j a j 1 + j,k∈S:j<k a j a k γ j γ k + j,k∈S:k<j a j a k γ j γ k = j∈S a 2 j 1.
From this we can directly conclude, that A has only two degenerate eigenvalues ±( j∈S a 2 j ) 1/2 such that can be directly computed. The same argument holds for the operators P γ j . Hence, we obtain C Scaling behaviour of support results for η = 0 Some features of the edge mode support can be inferred from analytical results computed for N → ∞.
Starting with we can attempt to construct the edge modes iteratively. We can infer from the structure of the Hamiltonian that the left edge modes will be of the form In the thermodynamic limit, the edge zero modes are expected to be exact. This imposes a condition on the coefficients of the zero modes via the commutator from which we can construct a recurrence relation that willl allow us to produce two linearly independant edge mode operators. Since we assume that these should be smoothly related to the perfectly localised operators when J = 0, we choose to begin the recurrance relation with either α 1 = 1 and α 2 = 0, which we identify with m 1 , or vice versa, which we identify with m 2 . In the case of m 1 , a clear scaling behaviour emerges from j = 4 onward while the coefficients of m 2 exhibit a similar scaling behaviour from j = 3 onward which predicts a significant amount of structure in our support measure plots. Since we can compute our support measure exactly in the free fermion context, we see that

D Finite size scaling and cross validation
In this appendix, we discuss the finite size scaling of the non-interacting code and also explain the cross validation between the two methods. The data for the finite size scaling for the non-interacting (η = 0) code is shown in Fig. 3. The color encodes support data of m 0 for different different system sizes and the gray line is machine precision. Each point is an average over 100 realizations for interaction and disorder values similar to the main text. We find that the different system sizes agree quite well in the parameter regime investigated here. Moreover, the data indicates that at a size of 24 sites the support decays to the machine precision, which is why we settled for a system size of that order for the main text material despite the availability of even larger systems.
For the interacting code things are quite different, as here 12 sites is the maximal system size that can be reached due to the effort needed for the sorting procedure. Furthermore, the emergent plateau structure which can be demonstrated in the infinite system limit for the non-interacting case (see App.,C) still persists in the interacting model. Thus, we are forced to effectively fit the exponential envelope to only three points, which unfortunately renders a system size scaling towards smaller systems meaningless. Nevertheless, the interacting procedure naturally also works if the Hamiltonian is actually non-interacting so we used this to at least cross validate results between both algorithms for small systems where they agree.

E Additional numerical data
In this appendix, we show additional numerical data obtained for the X L Z R edge mode and a bulk operator. Fig. 4a shows the support of X L Z R edge mode on the same scale as in the main text. It shows the same disorder dependence and has the same localization length as Z L X R . Fig. 4b shows the localization behaviour of a bulk operator in a chain of length N = 11. Here again, the disorder strength decreases the localization length.