Disorder in dissipation-induced topological states: Evidence for a different type of localization transition
DDisorder in dissipation-induced topological states: Evidence for novel localizationtransition
Alon Beck and Moshe Goldstein
Raymond and Beverly Sackler School of Physics and Astronomy, Tel Aviv University, Tel Aviv 6997801, Israel (Dated: November 20, 2020)The quest for nonequilibrium quantum phase transitions is often hampered by the tendency ofdriving and dissipation to give rise to an effective temperature, resulting in classical behavior. Inthis work we shed new light on this issue by studying the effect of disorder on recently-introduceddissipation-induced Chern topological states, and examining the eigenmodes of the Hermitian steadystate density matrix or entanglement Hamiltonian. We find that, similarly to equilibrium, eachLandau band has a single delocalized level near its center. However, using three different finite sizescaling methods we show that the critical exponent ν describing the divergence of the localizationlength upon approaching the delocalized state is significantly different from equilibrium if disorderis introduced into the non-dissipative part of the dynamics. This indicates a new nonequilibriumquantum critical universality class accessible in cold-atom experiments. Introduction.—
Recent years have seen a surge of inter-est in the driven-dissipative dynamics of quantum manybody systems [1, 2]. Of particular interest is the possibil-ity of new nonequilibrium quantum critical phenomena.However, typically far-from-equilibrium conditions giverise to an effective temperature governing the long timephysics, and leading to classical criticality. This stands inline with the usual perception of driving and dissipationas causing decoherence and destroying subtle quantumphenomena. However, this point of view has been chal-lenged by recent works showing how coupling to an envi-ronment could be engineered to drive a system towardsdesired steady states displaying quantum correlations [3–9], such as nonequilibrium topological states [10–23]. Inparticular, Refs. [18, 19] introduced a protocol, realizablewith cold atoms, for purely dissipative dynamics whichapproaches at a finite rate a mixed steady state as closeas desired to a pure topological state. Yet, the resultingtopology is encoded in the
Hermitian steady state den-sity matrix, giving rise to te same topological classes asin equilibrium [10–12, 15, 18, 21, 23–34]. Could this stilllead to new quantum nonequilibrium criticality?Every natural system exhibits imperfections and dis-order. In equilibrium, it has long been recognized thatdisorder is actually essential for stabilizing the most ba-sic topological phase, the integer quantum Hall state [35].Disorder localizes all states in a Landau level except oneat energy E c . The localization length diverges as oneapproaches it as [36, 37] ξ ( E ) ∼ | E − E c | − ν , (1)with a critical exponent ν governing the plateau tran-sition. Lately, a debate arose regarding the theoreticalvalue of ν [38–48], and its relation to experiment [49–51];the currently accepted value is 2.5–2.6.In this work we study the interplay between disorderand the recipe of Refs. [18, 19] for dissipatively-inducingChern-insulator states, through the effects of disorder onthe eigenmodes of the steady-state density matrix, whichis experimentally measurable in cold atoms [52–56]. This is thus a Hermitian localization problem, unrelated todisordered nonhermitian Hamiltonians [57, 58]. We showthat disorder in the system-bath coupling leads to thesame universality class as in equilibrium, but disorderperturbing the system Hamiltonian is not. We employthree different finite size scaling (FSS) methods, based on(a) the number of conducting states [48, 59]; (b) the lo-cal Chern marker [60]; (c) the transfer matrix Lyapunovexponent [37, 47]. The final results are presented in Ta-ble II; all methods show that the out-of-equilibrium ν is larger by 0.5–0.6 than equilibrium, hinting at a newuniversality class. Recipe.—
We now briefly recall the recipe for the dissi-pative creation of topological states, which is comprehen-sively described in Ref. [18]. Suppose we have a “refer-ence Hamiltonian”, H ref = (cid:80) i,j h ref ij c † i c j = (cid:80) λ ε ref λ c † λ c λ ( i, j being real space indexes in 2D, and λ an eigenvalueindex, which, in the clean case, would correspond to theband number and lattice momentum), with some de-sired (e.g., topologically-nontrivial) gapped ground statewhere only low-lying states ( λ ≤ λ ) are filled. Ratherthan implementing H ref as the system Hamiltonian, onemay set the system Hamiltonian to zero and employ dis-sipation to drive the system into a steady-state which isclose to the ground state of H ref . For this one takes a sys-tem consisting of two types of fermions (e.g., cold atomhyperfine states), with respective creation operators a † i (system) and b † i (bath). Both fermion species feel a lat-tice potential in the xy plane, but the bath b -fermionscould also escape in the z direction. Besides that, theHamiltonian of the a -fermions is trivial, ideally featur-ing no hopping; deviations from this will be describedby a system Hamiltonian H S = (cid:80) i,j h S,ij a † i a j . Rather,the dynamics originates from the system-bath couplingHamiltonian, which is built out of the matrix elements of a r X i v : . [ c ond - m a t . d i s - nn ] N ov Equilibrium Out of equilibriumMethod W Geometry
L L x N g M L eff x W µ ∗ γ in γ Geometry
L L x p N g M L eff x I 0.2 L × L a — 2 − . L × L b a —II 0.2 L × L − . L × L c —III 0.2 L × L x · — 5 10 − . L × L x . · M depends on L [61]. b N g = 25 for L ≤
49 and N g = 31 for L = 56 , c M = 3000 for L ≤
63 and M = 1500 for L = 70 , TABLE I. Methods parameters: W is the disorder strength, L and L x the system size in the y and x directions, respectively, N g the grid size (method I), M the number of disorder realizations, L eff x the effective x -length (method III), p the hopping rangecutoff (method III, nonequilibrium), µ ∗ the effective chemical potential, and γ in /γ the refilling rate in units of γ = 2 πν t . the reference Hamiltonian: H SB = (cid:88) i,j (cid:0) h ref ij − µ ∗ δ ij (cid:1) a † i b j + h . c . = (cid:88) λ (cid:0) ε ref λ − µ ∗ (cid:1) b † λ a λ + h . c . (2)The utility of the construction now becomes apparent:Suppose the lowest band of eigenstates of the referenceHamiltonian is almost flat (dispersionless). By tuning µ ∗ to its center ( ε ref λ ≈ µ ∗ for all λ ≤ λ ), its statesbecomes weakly coupled to the bath compared to statesin the other bands, λ > λ . Thus, all states are evapo-rated rapidly, except those belonging to the lowest band.One may then introduce another similar reservoir whichrefills all trapped states at a uniform rate to stabilize asteady state close to the ground state of the referenceHamiltonian, as we now explain.Integrating out the baths one gets a Lindblad [62] mas-ter equation, from which the Gaussian steady-state ρ canbe obtained. The latter is completely characterized bythe single-particle density matrix G ij ≡ tr( ρa † i a j ), whichobeys a continuous Lyapunov equation [18, 19, 63]: i [ G, h ∗ S ] + 12 (cid:8) G, γ out + γ in (cid:9) = γ in , (3)where γ in , γ out are nonnegative Hermitian matrices thatdescribes the rates which particles enter/escape of thesystem, and h ∗ S is the complex conjugate of the matrix h S . By Fermi’s golden rule, γ out λ = 2 πν ( ε ref λ − µ ∗ ) isdiagonal in the eigenbasis of the reference Hamiltonian[more generally, as a matrix γ out = 2 πν ( h ref − µ ∗ ) ],with ν the density of states of the b -species (assumedconstant), while γ in is taken as state independent (pro-portional to the unit matrix). For h S = 0, G is diag-onal in the eigenbasis of H ref , with eigenvalues n λ = γ in / ( γ in + γ out λ ) representing their mean occupation. Formax λ ≤ λ ( ε ref λ ) (cid:28) γ in (cid:28) min λ>λ ( ε ref λ ) we get n λ ≤ λ ≈ n λ>λ ≈
0, as desired. G [or, equivalently, the system-bath entanglement Hamiltonian − ln( ρ )] therefore has asimilar band structure to H ref , which is amenable totopological classification [10–12, 15, 18, 19, 28]. This mo-tivates the study of its eigenmodes and their localizationproperties in the presence of disorder. Localization transition.—
This work compares the lo-calization quantum phase transition of two systems. Thefirst is the equilibrium Hofstadter model [64] for the in-teger quantum Hall effect on a square lattice: H H = t (cid:88) r x ,r y e πiαr y a † r x +1 ,r y a r x ,r y + a † r x ,r y +1 a r x ,r y +h . c ., (4)where we take t = 1, α = 1 /
7. The second sys-tem is the out of equilibrium analogue, built using therecipe described above [18, 19]: The Hofstadter Hamil-tonian (whose lowest band is naturally almost-flat) istaken as the reference
Hamiltonian H ref = H H , while H S = 0. To study the localization phase transitionwe introduce disorder. In equilibrium we add a term H D = (cid:80) r x ,r y w r x ,r y a † r x ,r y a r x ,r y , where w r x ,r y ∈ [ − W, W ]are independent and uniformly distributed. Out of equi-librium, there are two options for introducing the samedisorder term: One may either add H D to H ref whilekeeping H S = 0, or keep H ref = H H and set H S = H D .We find that in both cases the disorder causes a nonequi-librium localization phase transition of the eigenmodesof G . Similarly to Eq. (1), ξ ( n ) ∝ | n − n c | − ν , where n is the occupation (eigenvalue of G ), and n c its criticalvalue. Does ν takes the same value as in equilibrium?In the first case the answer is yes; since H S = 0 we cansolve Eq. (3) explicitly: G − = 1 + 2 πν γ in ( h ref − µ ∗ ) . (5)Thus, even in the presence of disorder h ref and G stillshares the same eigenvectors , hence the same ν [61].This argument does not hold in the second scenario(disorder in H S ), since G and h ref have different eigen-vectors. Here we need to resort to numerical solutionof Eq. (3). We will investigate ν using three FSS meth-ods. For each we first calculate ν in equilibrium (dis-ordered Hofstadter model), and then out of equilibrium( H ref = H H and H S = H D ). Again, while in equilib-rium we examine the properties of the Hamiltonian (e.g.,eigenvector localization length, Chern number), out ofequilibrium we investigate G . The parameter values aresummarized in Table I, and the final results in Table II. (a) 𝑁 𝑁 𝑐 / 𝑁 𝑏 𝑁 R e s i du a l 𝜈 = 2.58 ± 0.04 R e s i du a l 𝑁 𝑁 𝑐 / 𝑁 𝑏 𝑁 (b) 𝜈 = 2.99 ± 0.10 Equilibrium Out of equilibrium
FIG. 1. Log-log plot of (cid:104) N c /N b (cid:105) as function L (systemsize), with N c the number of conducting states and N b = αL ( α = 1 /
7) the total number of states per band. Dashed linesrepresent linear fits with L ≥
28 in equilibrium and L ≥ Method I.—
Following Ref. [48, 61], we calculate thecritical exponent in equilibrium by the scaling of thenumber of conducting states, N c , N c ( L ) ∝ L − /ν , (6)where L is the system size and ν is the critical expo-nent. Working with a L × L Hofstadter model with pe-riodic boundary conditions, we calculate N c by countingthe number of single-particle states with nonzero Chernnumber, and average the result over M different disor-der realizations. In the presence of disorder, the Chernnumber can be defined as [65]: C = − π (cid:90) dθ x dθ y Im (cid:10) ∂ θ x ψ | ∂ θ y ψ (cid:11) dθ x dθ y , (7)where ψ ( θ x , θ y ) is the single-particle state and the inte-gral is over the space of twisted periodic boundary con-ditions, defined by the phases 0 ≤ θ x , θ y ≤ π . Forefficient calculation, we use the method suggested inRef. [66], employing grid size N g × N g [61]. Correctionsto the scaling (6) fade quickly with increasing the sys-tem size, hence may be ignored by excluding low sys-tem sizes. The nonequilibrium generalization is straight-forward: We calculate the Chern number of eigenstatesof G (instead of H ) by introducing the twisted boundaryconditions θ x , θ y into H ref . Then, we count the conduct-ing states within the band of highest occupations (whichcorresponds to the lowest Landau band of H ). Resultsare presented in Fig. 1. Method II.—
Here we study FSS of the topological in-dex [36, 61]. In equilibrium, in the vicinity of the criticalenergy E c , the total Chern number C L ( E ) of the statesbelow energy E scales as C L ( E ) = f (cid:0) ( E − E c ) L /ν (cid:1) ,with a scaling function f . To calculate C L ( E ) we use thelocal Chern marker [60, 67] with open boundary condi-tions, C r x ,r y = − πi (cid:104) r x , r y | ˜ X ˜ Y − ˜ Y ˜ X | r x , r y (cid:105) , (8)where ˜ X, ˜ Y are the projected lattice position operators:˜ X = P XP, ˜ Y = P Y P , P being a projection onto (a) 𝐸 𝐶 (b) 𝑛 𝐶 𝐸 − 𝐸 𝑐 𝐿 𝑛 − 𝑛 𝑐 𝐿 𝜈 = 2.26 ± 0.04 𝜈 = 2.91 ± 0.06 FIG. 2. Averaged local Chern number (a) in and (b) out ofequilibrium. Insets: scaling data collapse. the occupied states. The local Chern marker fluctuatesaround the value of the Chern number in the bulk ofthe system, but takes different values on the edges, since (cid:80) r x ,r y C r x ,r y = 0. Thus, we average C r x ,r y over the bulk,while excluding 1/4 of the sample length from each side, C = (4 /L ) × (cid:80) L/ ≤ r x ,r y ≤ L/ C r x ,r y , and also evaluatethe statistical error ∆ C . As in method I, irrelevant cor-rections exist, but their influence decreases rapidly withincreasing system size. We then search for ν, E c , and thecoefficients of a polynomial approximating f [61], whichminimize the chi-squared: χ = (cid:88) E,L (cid:34) f (cid:0) ( E − E c ) L /ν (cid:1) − C L ( E )∆ C L ( E ) (cid:35) , (9)Out of equilibrium, we calculate C as a function of n (eigenvalue of G ) instead of E . To obtain C ( n ) fromEq. (8) we use a projection P ( n ) to states with occupa-tion higher than n . The results are presented in Fig. 2. Method III.—
Here we perform FSS of the localiza-tion length ξ . Following Ref. [61, 68], we calculate thelocalization length with the transfer-matrix method: Weconsider a long cylinder of size L x × L , L x (cid:29) L . Let ψ bean eigenvalue of the Hamiltonian with energy E . Fromthe equation Hψ = Eψ we can construct the 2 L × L transfer-matrix T r x , defined as: (cid:32) ψ r x +1 ψ r x (cid:33) = T r x (cid:32) ψ r x ψ r x − (cid:33) , (10)where ψ r x is a vector with elements ψ r x ,r y =1 ··· L . Beingsymplectic, the eigenvalues of each transfer matrix comein reciprocal pairs { λ, λ − } . The same applies to theirproduct, T = (cid:81) L x r x =1 T r x . The Lyapunov exponent (in-verse localization length) is defined as:˜Λ ≡ ξ − = lim L x →∞ ln( λ min ) L x , (11)where λ min is the smallest eigenvalue of T that is largerthan unity. We have applied the Gram-Schmidt processto the columns of T every 7 multiplications to reducenumerical error. The results are presented in Fig. 3(a).As in the previous method, ν can be extracted by chi-squared minimization (9) for the dimensionless Lyapunov (a) 𝐸 Λ (b) 𝑛 −1 Λ Number of excluded points 𝜈 (d) 𝐿 eff 𝐺 −1 … 𝐺 −1 𝐿 𝑥 𝐿 𝐺 −1 𝐺 −1 𝐺 −1 (c) Equilibrium Out of equilibrium
FIG. 3. Dimensionless Lyapunov exponent (a) in and (b) outof equilibrium. Insets: Zoom into the vicinity of the criticalpoint. (c) Illustration of the nonequilibrium transfer matri-ces construction. (d) Comparison of the critical exponent inand out of equilibrium, without corrections to scaling. Thehorizontal axis represents the number data points that wereexcluded from each side of the critical point in the chi-squaredminimization. exponent Λ ≡ L ˜Λ. However, since the data containsstrong corrections to scaling (typical for the long cylin-der geometry), we account for a single irrelevant scalingfield [61].The nonequilibrium generalization from Λ( E ) to Λ( n )is more complicated compared to the previous methods.First, unlike H , G has non-local hopping terms whichprevent us from constructing a transfer matrix. This re-quires introducing a cutoff p on the hopping range inthe x direction, and setting terms of range larger than p to zero. From this perspective it is advantageous toconstruct the transfer matrices using G − , since Eq. (5)shows that for H S = 0 its elements have a finite range p = 2. We have verified numerically that the elements of G − decay exponentially with range for H S = H D , mak-ing truncation at p = 5 a very good approximation [61].A second issue is that in the presence of disorder thestructure of G − can only be obtained numerically, bysolving Eq. (3). Thus, we cannot analytically obtain thetransfer matrix at a specific x -position, and instead, wecan only generate the entire G − matrix, which is im-practical for L x (cid:29)
1. As a solution, we use the schemedepicted in Fig. 3(c): We generate a G − matrix of size L x × L for some large but practical L x (with periodicboundary conditions) [61]. We repeat this with M dis-order realizations, and denote the resulting matrices as (cid:8)(cid:0) G − (cid:1) m (cid:9) Mm =1 . From each (cid:0) G − (cid:1) m we extract K trans-fer matrices ( K = L x − c , excluding the c = 7 matricesclosest to each end) by putting a cutoff p on the hoppingrange, as explained above. We then define the sequence { T n } MKn =1 , with T ( n − K +1 , ..., T nK the transfer matrices Method I II IIIEquilibrium 2 . ± .
04 2 . ± .
04 2 . ± . . ± .
10 2 . ± .
06 not convergent, higherthan equilibriumTABLE II. Summary of the results for the critical exponent ν , in and out of equilibrium. extracted from (cid:0) G − (cid:1) n . The effective system lengthis thus L eff = M K . The mismatch between transfer-matrices that originate from different G − (for example, T K and T K +1 ) introduces an error, but it can be reducedby increasing L x .The results are shown in Fig. 3(b). The numericaleffort per sample is still much higher in the nonequilib-rium case, limiting our ability to reduce statistical errorby either sample averaging or using large system sizes.Hence, we can neither implement corrections to scalingnor drop small systems, and therefore cannot determine ν as accurately as before. We thus resolve to extractingthe uncorrected nonequilibrium exponent and comparingit with a similarly obtained equilibrium value, to appre-ciate the significance of their difference, see Fig. 3(d). Results and Discussion.—
The results are summarizedin Table II. In equilibrium they are generally in line withprevious studies. For method I, the obtained ν = 2 . ± .
04 is somewhat higher than the value ν = 2 . ± . α = 1 / ν = 2 . ± .
04) is smaller than recentestimates of the critical exponent, which seems to be ageneral feature of FSS of a topological index [46, 69]. Inmethod III, upon including corrections to scaling we get ν = 2 . ± . y = 0 . ± .
01, Λ( E c ) = 0 . ± .
01, with y the leading irrelevant exponent. This is slightly smallerbut still in agreement with ν = 2 . ± .
03 obtained inRef. [47] for the Hofstadter model.Out of equilibrium, methods I and II give rise to valueswhich are significantly higher than in equilibrium. Theresults of method III are not convergent, but they stillstrongly suggest that ν is higher than equilibrium by 0.5–0.6, in agreement with the other methods.Finally, let us reiterate that the single-particle densitymatrix G is Hermitian. Furthermore, G − is local inspace. The locality is exact for h S = 0, where G − isessentially the square of h ref , see Eq. (5). We have foundthat for disorder in H S the elements of G − have dis-tributions without fat tails, and with averages and cor-relations which decay exponentially with distance [61].Thus, our results indicate a new universality class of thelocal Hermitian disordered G − , which is rooted in thenonequilibrium nature of the system. Conclusions.—
In this work we have investigatedthe effects of disorder on dissipation-induced topologi-cal states. We demonstrated the existence of an out ofequilibrium localization phase transition similar to the in-teger quantum Hall plateau transition. Using three FSSmethods, we found a significant difference between thevalue of the critical exponent ν in and out of equilibriumwhen disorder is introduced into the non-dissipative partof the Lindbladian. This indicates a novel nonequilibriumquantum universality class, despite the steady state den-sity matrix being Hermitian and local. These findingscould be tested in cold-atom experiments, which allow toimplement the dissipative protocol [18, 19], and to mea-sure the single-particle density matrix [52–56]. In thefuture it would be interesting to try to attack the prob-lem using field theoretical methods [2, 37], and to study the relation between the steady state and the nonher-mitian [57, 58, 70] decay towards it (a relation which isnontrivial out of equilibrium [19]), as well as the possi-bility of new many-body localization transition [71, 72]. ACKNOWLEDGMENTS
We thank I. Burmistrov, R. Ilan, and E. Shimshoni foruseful discussions. Support by the Israel Science Foun-dation (Grant No. 227/15) and the US-Israel BinationalScience Foundation (Grant No. 2016224) is gratefullyacknowledged. [1] A. Kamenev,
Field Theory of Non-Equilibrium Systems (Cambridge University Press, 2009).[2] L. M. Sieberer, M. Buchhold, and S. Diehl, Keldysh fieldtheory for driven open quantum systems, Reports onProgress in Physics , 096001 (2016).[3] S. Diehl, A. Micheli, A. Kantian, B. Kraus, H. P. Buchler,and P. Zoller, Quantum states and phases in driven openquantum systems with cold atoms, Nature Physics , 878(2008).[4] B. Kraus, H. P. Buchler, S. Diehl, A. Kantian, A. Micheli,and P. Zoller, Preparation of entangled states by quan-tum markov processes, Physical Review A (2008).[5] F. Verstraete, M. M. Wolf, and J. I. Cirac, Quantumcomputation and quantum-state engineering driven bydissipation, Nature Physics , 633 (2009).[6] H. Weimer, M. Muller, I. Lesanovsky, P. Zoller, and H. P.Buchler, A Rydberg quantum simulator, Nature Physics , 382 (2010).[7] J. Otterbach and M. Lemeshko, Dissipative preparationof spatial order in rydberg-dressed bose-einstein conden-sates, Physical Review Letters (2014).[8] N. Lang and H. P. Buchler, Exploring quantum phasesby driven dissipation, Physical Review A (2015).[9] L. Zhou, S. Choi, and M. D. Lukin, Symmetry-protected dissipative preparation of matrix productstates, arXiv:1706.01995 [quant-ph] (2017).[10] S. Diehl, E. Rico, M. A. Baranov, and P. Zoller, Topologyby dissipation in atomic quantum wires, Nature Physics , 971 (2011).[11] C.-E. Bardyn, M. A. Baranov, E. Rico, A. ˙Imamo˘glu,P. Zoller, and S. Diehl, Majorana modes in driven-dissipative atomic superfluids with a zero Chern number,Physical Review Letters (2012).[12] C.-E. Bardyn, M. A. Baranov, C. V. Kraus, E. Rico,A. Imamoglu, P. Zoller, and S. Diehl, Topology by dissi-pation, New Journal of Physics , 085001 (2013).[13] R. Konig and F. Pastawski, Generating topological order:No speedup by dissipation, Physical Review B (2014).[14] E. Kapit, M. Hafezi, and S. H. Simon, Induced self-stabilization in fractional quantum Hall states of light,Physical Review X (2014).[15] J. C. Budich, P. Zoller, and S. Diehl, Dissipative prepa-ration of Chern insulators, Physical Review A (2015).[16] F. Iemini, D. Rossini, R. Fazio, S. Diehl, and L. Mazza,Dissipative topological superconductors in number- conserving systems, Physical Review B (2016).[17] Z. Gong, S. Higashikawa, and M. Ueda, Zeno Hall effect,Physical Review Letters (2017).[18] M. Goldstein, Dissipation-induced topological insulators:A no-go theorem and a recipe, SciPost Physics (2019).[19] G. Shavit and M. Goldstein, Topology by dissipation:Transport properties, Physical Review B (2020).[20] F. Tonielli, J. C. Budich, A. Altland, and S. Diehl, Topo-logical field theory far from equilibrium, Phys. Rev. Lett. , 240404 (2020).[21] T. Yoshida, K. Kudo, H. Katsura, and Y. Hatsugai, Fateof fractional quantum Hall states in open quantum sys-tems: Characterization of correlated topological statesfor the full Liouvillian, Phys. Rev. Research , 033428(2020).[22] S. Bandyopadhyay and A. Dutta, Dissipative prepa-ration of many-body Floquet Chern insulators,arXiv:2005.09972 [cond-mat.stat-mech] (2020).[23] A. Altland, M. Fleischhauer, and S. Diehl, Sym-metry classes of open fermionic quantum matter,arXiv:2007.10448 [cond-mat.str-el] (2020).[24] A. Rivas, O. Viyuela, and M. A. Martin-Delgado,Density-matrix Chern insulators: Finite-temperaturegeneralization of topological insulators, Phys. Rev. B ,155141 (2013).[25] Z. Huang and D. P. Arovas, Topological indices for openand thermal systems via Uhlmann’s phase, Phys. Rev.Lett. , 076407 (2014).[26] O. Viyuela, A. Rivas, and M. A. Martin-Delgado, Two-dimensional density-matrix topological fermionic phases:Topological Uhlmann numbers, Phys. Rev. Lett. ,076408 (2014).[27] E. P. L. van Nieuwenburg and S. D. Huber, Classificationof mixed-state topology in one dimension, Phys. Rev. B , 075141 (2014).[28] J. C. Budich and S. Diehl, Topology of density matrices,Phys. Rev. B , 165140 (2015).[29] F. Grusdt, Topological order of mixed states in correlatedquantum many-body systems, Phys. Rev. B , 075106(2017).[30] C.-E. Bardyn, A recipe for topological observables ofdensity matrices, arXiv:1711.09735 [cond-mat.quant-gas](2017).[31] C.-E. Bardyn, L. Wawer, A. Altland, M. Fleischhauer,and S. Diehl, Probing the topology of density matrices, Phys. Rev. X , 011035 (2018).[32] D.-J. Zhang and J. Gong, Topological characterizationof one-dimensional open fermionic systems, Phys. Rev.A , 052101 (2018).[33] A. Coser and D. P´erez-Garc´ıa, Classification of phasesfor mixed states via fast dissipative evolution, Quantum , 174 (2019).[34] S. Lieu, M. McGinley, and N. R. Cooper, Tenfold wayfor quadratic Lindbladians, Phys. Rev. Lett. , 040401(2020).[35] K. v. Klitzing, G. Dorda, and M. Pepper, New methodfor high-accuracy determination of the fine-structure con-stant based on quantized Hall resistance, Phys. Rev. Lett. , 494 (1980).[36] B. Huckestein, Scaling theory of the integer quantum Halleffect, Reviews of Modern Physics , 357 (1995).[37] F. Evers and A. D. Mirlin, Anderson transitions, Reviewsof Modern Physics , 1355 (2008).[38] K. Slevin and T. Ohtsuki, Critical exponent for the quan-tum Hall transition, Physical Review B (2009).[39] H. Obuse, A. R. Subramaniam, A. Furusaki, I. A.Gruzberg, and A. W. W. Ludwig, Conformal invariance,multifractality, and finite-size scaling at anderson local-ization transitions in two dimensions, Physical Review B (2010).[40] M. Amado, A. V. Malyshev, A. Sedrakyan, andF. Dom´ınguez-Adame, Numerical study of the localiza-tion length critical index in a network model of plateau-plateau transitions in the quantum Hall effect, PhysicalReview Letters (2011).[41] I. C. Fulga, F. Hassler, A. R. Akhmerov, and C. W. J.Beenakker, Topological quantum number and critical ex-ponent from conductance fluctuations at the quantumHall plateau transition, Physical Review B (2011).[42] K. Slevin and T. Ohtsuki, Finite size scaling of theChalker-Coddington model, International Journal ofModern Physics: Conference Series , 60 (2012).[43] H. Obuse, I. A. Gruzberg, and F. Evers, Finite-size ef-fects and irrelevant corrections to scaling near the integerquantum Hall transition, Physical Review Letters (2012).[44] W. Nuding, A. Klumper, and A. Sedrakyan, Localizationlength index and subleading corrections in a Chalker-Coddington model: A numerical study, Physical ReviewB (2015).[45] I. A. Gruzberg, A. Klumper, W. Nuding, and A. Se-drakyan, Geometrically disordered network models,quenched quantum gravity, and critical behavior at quan-tum Hall plateau transitions, Physical Review B (2017).[46] M. Ippoliti, S. D. Geraedts, and R. N. Bhatt, Integerquantum Hall transition in a fraction of a landau level,Physical Review B (2018).[47] M. Puschmann, P. Cain, M. Schreiber, and T. Vojta, In-teger quantum Hall transition on a tight-binding lattice,Physical Review B (2019).[48] Q. Zhu, P. Wu, R. N. Bhatt, and X. Wan, Localization-length exponent in two models of quantum Hall plateautransitions, Physical Review B (2019).[49] W. Li, G. A. Cs´athy, D. C. Tsui, L. N. Pfeiffer, andK. W. West, Scaling and universality of integer quantumHall plateau-to-plateau transitions, Phys. Rev. Lett. ,206807 (2005).[50] W. Li, C. L. Vicente, J. S. Xia, W. Pan, D. C. Tsui, L. N. Pfeiffer, and K. W. West, Scaling in plateau-to-plateautransition: A direct connection of quantum Hall systemswith the anderson localization model, Phys. Rev. Lett. , 216801 (2009).[51] A. J. M. Giesbers, U. Zeitler, L. A. Ponomarenko,R. Yang, K. S. Novoselov, A. K. Geim, and J. C. Maan,Scaling of the quantum Hall plateau-plateau transitionin graphene, Phys. Rev. B , 241411 (2009).[52] P. Hauke, M. Lewenstein, and A. Eckardt, Tomography ofband insulators from quench dynamics, Phys. Rev. Lett. , 045303 (2014).[53] N. Fl¨aschner, B. S. Rem, M. Tarnowski, D. Vogel, D.-S. L¨uhmann, K. Sengstock, and C. Weitenberg, Experi-mental reconstruction of the berry curvature in a FloquetBloch band, Science , 1091 (2016).[54] M. Tarnowski, M. Nuske, N. Fl¨aschner, B. Rem, D. Vo-gel, L. Freystatzky, K. Sengstock, L. Mathey, andC. Weitenberg, Observation of topological Bloch-statedefects and their merging transition, Phys. Rev. Lett. , 240403 (2017).[55] L. A. Pe˜na Ardila, M. Heyl, and A. Eckardt, Measuringthe single-particle density matrix for fermions and hard-core bosons in an optical lattice, Phys. Rev. Lett. ,260401 (2018).[56] J.-H. Zheng, B. Irsigler, L. Jiang, C. Weitenberg, andW. Hofstetter, Measuring an interaction-induced topo-logical phase transition via the single-particle density ma-trix, Phys. Rev. A , 013631 (2020).[57] N. Hatano and D. R. Nelson, Localization transitions innon-Hermitian quantum mechanics, Phys. Rev. Lett. ,570 (1996).[58] Y. Ashida, Z. Gong, and M. Ueda, Non-Hermitianphysics, arXiv:2006.01837 [cond-mat.mes-hall] (2020).[59] K. Yang and R. N. Bhatt, Floating of extended states andlocalization transition in a weak magnetic field, PhysicalReview Letters , 1316 (1996).[60] R. Bianco and R. Resta, Mapping topological order incoordinate space, Physical Review B (2011).[61] See the supplmental material [URL] for technical details.[62] P. Z. Crispin Gardiner, Quantum Noise (Springer BerlinHeidelberg, 2004).[63] F. Schwarz, M. Goldstein, A. Dorda, E. Arrigoni, A. We-ichselbaum, and J. von Delft, Lindblad-driven discretizedleads for nonequilibrium steady-state transport in quan-tum impurity models: Recovering the continuum limit,Phys. Rev. B , 155142 (2016).[64] D. R. Hofstadter, Energy levels and wave functions ofBloch electrons in rational and irrational magnetic fields,Physical Review B , 2239 (1976).[65] Q. Niu, D. J. Thouless, and Y.-S. Wu, Quantized Hallconductance as a topological invariant, Physical ReviewB , 3372 (1985).[66] T. Fukui, Y. Hatsugai, and H. Suzuki, Chern numbers indiscretized Brillouin zone: Efficient method of computing(spin) Hall conductances, Journal of the Physical Societyof Japan , 1674 (2005).[67] M. D. Caio, G. Moller, N. R. Cooper, and M. J. Bhaseen,Topological marker currents in Chern insulators, NaturePhysics , 257 (2019).[68] A. MacKinnon and B. Kramer, The scaling theory of elec-trons in disordered solids: Additional numerical results,Zeitschrift fur Physik B Condensed Matter , 1 (1983).[69] T. A. Loring and M. B. Hastings, Disordered topologi-cal insulators via C ∗ -algebras, EPL (Europhysics Letters) , 67004 (2010).[70] N. Silberstein, J. Behrends, M. Goldstein, and R. Ilan,Berry connection induced anomalous wave-packet dy-namics in non-Hermitian systems, arXiv:2004.13746[cond-mat.mes-hall] (2020).[71] R. Nandkishore and D. A. Huse, Many-body localization and thermalization in quantum statistical mechanics, An-nual Review of Condensed Matter Physics , 15 (2015).[72] E. Altman and R. Vosk, Universal dynamics and renor-malization in many-body-localized systems, Annual Re-view of Condensed Matter Physics , 383 (2015). SUPPLEMENTAL MATERIAL FOR: “DISORDER IN DISSIPATION-INDUCED TOPOLOGICALSTATES: EVIDENCE FOR NOVEL LOCALIZATION TRANSITION”
In this Supplemental Material we provide additional technical details and results. In Sec. S.I we display an argumentto the fact that when the disorder appears in the system-bath coupling Hamiltonian, the localization phase transitionis in the same universality class as in equilibrium. In Secs. S.II–S.IV we present additional details regarding the threemethods which were used to calculate the critical exponent ν out of equilibrium. Finally, in Sec. S.V we characterizethe distribution and correlation of the elements of the matrix G − . S.I. DISORDER IN THE DISSIPATIVE DYNAMICS
In the main text we have stated that if we consider disorder only in the reference Hamiltonian h ref (that is, thedynamics is purely-dissipative, H S = 0), then the critical exponent of the phase transition is the same as in equilibrium(that is, in the same universality class). We now present a more detailed argument for this. In fact, it is a specialcase of the following claim: Claim.
Let h be a nondegenerate Hamiltonian with a property ξ h ( E ), which is a function of the eigenstates of h (with energy E ). Suppose we have a phase transition described by a scaling law, ξ h ( E ) ∝ | E − E c | − ν , where E c isthe critical energy and ν is the critical exponent. Then, any analytic function G = g ( h ) of h that satisfies0 < dgdh (cid:12)(cid:12)(cid:12)(cid:12) E c < ∞ , (S1)will display a phase transition with the same critical exponent. That is, ξ G ( n ) ∝ | n − n c | − ν , where n represent aneigenvalue of G . Proof.
We notice that G has the same eigenstates as h , but with different eigenvalues described by the relation n ( E ) = g ( E ), where n ( E ) is the eigenvalue of G corresponding to the eigenvalue E of h . Since ξ is determined onlyby the eigenstates, we have ξ G ( n ( E )) = ξ h ( E ) , (S2)and thus ξ G ( n ) ∝ = | E − E c | − ν = | g − ( n ) − g − ( n c ) | − ν , (S3)where g − is the inverse of the function g . We note that g − is well defined around n c since 0 < dgdh (cid:12)(cid:12)(cid:12) E c < ∞ .Expanding g − to first order around n c , we get ξ G ( n ) ≈ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) dg − dn (cid:12)(cid:12)(cid:12)(cid:12) n c ( n − n c ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − ν , (S4)hence ξ ( n ) ∝ | n − n c | − ν , as expected.Going back to our case, Eq. (5) of the main text shows that for H S = 0 the single-particle density matrix G can beexpressed as a function of h ref , whose derivative is dGdh ref = − G πν γ in h ref − µ ∗ ) , (S5)that is, condition (S1) will be satisfied for µ ∗ (cid:54) = E c . Equilibrium Out of equilibrium L min ν χ L min ν χ . ± .
01 46.4 7 3 . ± .
01 1814 2 . ± .
02 0.92 14 3 . ± .
03 1.221 2 . ± .
03 0.94 21 3 . ± .
05 1.228 2 . ± .
04 0.57 28 3 . ± .
06 1.035 2 . ± .
07 0.32 35 2 . ± .
10 0.2842 2 . ± .
12 0.46 42 2 . ± .
19 0.41TABLE ST1. The results of the critical exponent ν without correction to scaling, extracted from fits of the number of conductingstates, Eq. (S10). The smallest system size is taken as L min and the largest system size is always L max = 63. S.II. ADDITIONAL DETAILS FOR METHOD I
Calculation of the Chern number.
For efficient calculation of the Chern number we have followed the method ofRef. [66]. We divide the parameter space 0 ≤ θ x , θ y ≤ π into a grid of size N g × N g with equal spacing. For eachpoint, θ = ( θ x , θ y ) = 2 πN g ( r x , r y ) , r x , r y = 0 , · · · , N g − , (S6)we define: U ˆ µ ( θ ) ≡ (cid:104) ψ ( θ ) | ψ ( θ + ˆ µ ) (cid:105) /N ˆ µ ( θ ) , (S7)where N ˆ µ ( θ ) ≡ |(cid:104) ψ ( θ ) | ψ ( θ + ˆ µ ) (cid:105)| is a normalization factor, and ˆ µ = ˆ x, ˆ y is a grid lattice vector in the x or y direction,respectively. We define a discretized version of the Berry curvature, F ( θ ) ≡ ln (cid:2) U x ( θ ) U y ( θ + ˆ x ) U x ( θ + ˆ y ) − U y ( θ ) − (cid:3) , (S8)where the principal branch of the logarithm is used. Finally, the Chern number is defined as C ≡ πi (cid:88) θ F ( θ ) , (S9)which must result in an integer value. However, the result might contain an error if N g is not large enough. In ourcase we can detect errors by checking that the sum of the Chern numbers of all the single-particle states equals − −
1, we know that at least one errorhas occurred in the calculation and therefore reject it. We have taken values of N g which would keep the rejectionrate smaller than 2%: In equilibrium, we choose N g = 30 for all of the system sizes. Out of equilibrium, we choose N g = 25 for L = 7 , ...,
49, and N g = 31 for L = 56 , Calculation of the critical exponent.
Unlike an infinite system, in which an extended (conducting) state exists onlyat a single energy E c , in a finite-sized system there is a range of extended states, corresponding to the range of energieswith localization lengths ξ ( E ) > L . Thus, to obtain the finite-size mobility edges need to solve ξ ( E ) = ξ | E − E c | − ν = L , leading to E , = E c ± ( ξ /L ) /ν . Therefore, the number of conducting states would be N c = L (cid:82) E E ρ ( E ) dE, where ρ ( E ) is the density of states of the band. Approximating ρ ( E ) ≈ ρ ( E c ) and recalling that the total number of statesin a Landau band is N b = αL ( α = 1 / N c N b = aL − /ν , (S10)where a is some constant. ν can be extracted by numerical calculations of N c for different system sizes. The results(without corrections to scaling) in and out of equilibrium are presented in Table ST1. In equilibrium, the correctionsto scaling are significant only when L min = 7. The first correction can be included by considering the generalizedscaling form: N c N b = a (cid:0) bL − y (cid:1) L − /ν , (S11) Equilibrium Out of equilibrium L min ν χ L min ν χ . ± .
007 3.5 7 3 . ± .
01 2514 2 . ± .
01 3.2 14 3 . ± .
03 1.721 2 . ± .
01 3.8 21 3 . ± .
04 1.828 2 . ± .
02 4.4 28 3 . ± .
05 0.6735 2 . ± .
03 5.1 35 3 . ± .
09 0.4542 2 . ± .
06 6.4 42 3 . ± .
15 0.34TABLE ST2. The results of the critical exponent ν without correction to scaling, extracted from fits of the width of conductingstates density, ∆ E c [cf. Eq. (S12)]. The smallest system size is taken as L min and the largest system size is always L max = 63. where y > L = 7 , ..., ν = 2 . ± . y = 4 . ± . χ = 0 . L min = 28 in equilibrium and L min = 35 out of equilibrium. This choice ensures that the three lowest system sizesare excluded to avoid corrections to scaling, but also that the change in ν between the employed L min and using thefollowing value L min + 1 /α = L min + 7 is smaller than the uncertainty in ν .An additional way for extracting the critical exponent is by looking at the width of the density of the conductingstates, ρ c ( E ), defined as ∆ E c = L N c (cid:90) ρ c ( E ) E dE − (cid:20) L N c (cid:90) ρ c ( E ) EdE (cid:21) , (S12)which is expected to scale as ∆ E c ∼ L − /ν . The results without corrections are presented in Table ST2. While theyalso suggest a higher value of the critical exponent out of equilibrium, they seem to be less reliable than the resultswith the number of conducting states, as evidenced by the significantly large chi-squared values. This may imply thatthe width of the distribution is more sensitive to finite-size corrections than the number of conducting states (See thediscussion around Fig. 7 of Ref. [48]). S.III. ADDITIONAL DETAILS FOR METHOD II
The results of the local Chern marker calculations are presented in Fig. SF1. Corrections to scaling are presentat low system sizes. This can be seen in the insets, which show that the low systems sizes curves cross the largesystem curves away from E c . The corrections to scaling for small system sizes turn out be be difficult to fit accurately.However, they also decrease rapidly with increasing system sizes. We therefore resort to omitting smaller system sizesand ignoring corrections to scaling. We note that the result also depends on the range of energies (or occupations outof equilibrium) that were included in the fit. Therefore, we chose a range of energies (occupations) in which the resultis most stable (least sensitive to increasing or decreasing the number of included points). For example, for L min = 35 Equilibrium Out of equilibrium L min ν χ L min ν χ
14 2 . ± .
06 40 . . ± .
04 2521 2 . ± .
05 17 . . ± .
06 1228 2 . ± .
05 8 . . ± .
07 5.435 2 . ± .
04 5 . . ± .
06 3.542 2 . ± .
04 6 . . ± .
08 449 2 . ± .
05 5 . . ± .
10 4TABLE ST3. The results of the critical exponent ν , extracted from method II. The smallest system size is taken as L min andthe largest system size is always L max = 77. 𝐶 𝐸 Equilibrium Out of equilibrium(a) (b) 𝑛 𝐶 FIG. SF1. The local Chern marker for system sizes L = 14 , ...,
77, (a) in equilibrium (results averaged over 3 · disorderrealizations); (b) out of equilibrium (results averaged over 3 · disorder realizations). Inset: zoom-in onto the vicinity of thecritical point. we have estimated ν = 2 . ± .
04 from the results presented in Fig. SF2(a). As for the degree D of the polynomial f ( x ) = (cid:80) Dq =0 a q x q that was used to approximate the scaling function f , we have verified that it is large enough tocapture the behavior in the given range, but is not too large, so as to prevent overfitting. We have thus used D = 5 inequilibrium and D = 7 out of equilibrium. The results in and out of equilibrium are shown in Table ST3. In order toavoid correction to scaling effects, we have chosen to exclude the four first system sizes. That is, the lowest includedsystem size was chosen as L min = 35 in and out of equilibrium. As in method I, we also verified that the change in ν with respect to using the following value L min + 1 /α = L min + 7 is smaller than the uncertainty in ν .To verify these results, we have also extracted the critical exponent from the L -dependence of the derivative of C at the critical point. For example, in equilibrium, since C L ( E ) = f (cid:16) ( E − E c ) L ν (cid:17) , we haveln (cid:32) ∂C L ( E ) ∂E (cid:12)(cid:12)(cid:12)(cid:12) E c (cid:33) = ln ( f (cid:48) (0)) + 1 ν ln( L ) . (S13)By fitting a polynomial expansion to each C L ( E ) data set at fixed L , we can obtain the left hand side of the lastequation by approximating ∂C L ( E ) ∂E (cid:12)(cid:12)(cid:12) E c ≈ max (cid:16) ∂C L ( E ) ∂E (cid:17) . Then, ν can be extracted from a linear fit of the logarithmof the latter quantity as function of ln( L ). While we found this method to be less stable, its results were still inagreement with the chi-square minimization results: In equilibrium we got ν ≈ . − .
3, while out of equilibrium wegot ν ≈ . − . 𝜈 Points excluded from each side E l o c a l C h e r nnu m b e r (a) (b) FIG. SF2. Determination of ν using method II in equilibrium: (a) Values of ν extracted from system sizes L = 35 , , ..., E c for each data set inpanel (b). The critical energy was found to be E c ≈ − . 𝜈 Points excluded from each side 𝜈 Points excluded from each side (a) (b)
FIG. SF3. Values of ν extracted in equilibrium using method III and system sizes L = 42 , , , , , , , , , E c for each data set inFig. 3(a) of the main text. Panels (a) and (b) correspond, respectively, to approaches (a) and (b) described in the text. S.IV. ADDITIONAL DETAILS FOR METHOD III
Equilibrium.
We will present additional details regarding the corrections to scaling that were used in the fittingprocedure in equilibrium. We assume the existence of only one irrelevant exponent (including more than a singleirrelevant exponent would first, be a numerical challenge which in general requires data with much lower uncertainties,and second, seems not to be required in practice for the system sizes used). That is, we assume the following scalingform: Lξ ( E ) = f ( u r L /ν , u i L − y ) , where ξ ( E ) is the localization length, f is some scaling function, u r , u i are the relevant and irrelevant scaling fields,respectively, and y > E c as u r ( E − E c ) = (cid:80) m r n =1 a n ( E − E c ) n , u i ( E − E c ) = (cid:80) m i n =0 b n ( E − E c ) n (the term n = 0 is absent for the relevant field since it mustvanish at the critical point). In addition, we expand f to the first order in the irrelevant field: f (cid:16) u r L /ν , u i L − y (cid:17) ≈ f (cid:16) u r L /ν (cid:17) + u i L − y f (cid:16) u r L /ν (cid:17) , where f , f are some single-parameter functions. We will now present two approaches which lead to similar results:(i) We assume a simple form of the scaling fields: u r ( E ) = E − E c , u i ( E ) = 1, and take f , f as the followingpolynomials: f ( x ) = (cid:80) n n =0 a n x n , f ( x ) = (cid:80) n n =0 b n x n , with n = 5 , n = 4.(ii) Following Ref. [47], we consider only even terms in the scaling functions: f ( x ) = (cid:80) n =0 a n x n , f ( x ) = (cid:80) n =0 b n x n , where n = 3, n = 2. However, we include additional terms in the expansion of the scaling fields: u r ( E ) = (cid:80) m n =1 c n ( E − E c ) n , u i ( E ) = (cid:80) m n =0 d n ( E − E c ) n , with m = 3, m = 1. This can be motivated by the factthat our data is close to being an even function around E c , and the small asymmetry is reflected by the odd terms ofthe expansion of the scaling fields.We have found both approaches to have low sensitivity to the choice of the smallest system size to be included inthe fit, but are still somewhat affected by the range of energies taken around the critical energy E c . As in method II,we chose a range of energies in which the result is most stable (least sensitive to increasing or decreasing the numberof included points). We also verified that increasing n , n , m , m has small impact on the results. A comparison ofthe results is presented in Fig. SF3. We have found approach (b) to be slightly more stable. Out of equilibrium.
We will provide here additional details regarding the transfer-matrix method and the choiceof the parameters that were used to extract the localization length out of equilibrium. As described in the maintext, we have calculated the matrix G − (of size L x × L , with periodic boundary conditions in both directions), for M different disorder realizations. We then set the hopping terms in G − with range along the x -direction largerthan a cutoff p to be zero. Since the matrix now has only a finite hopping range, it is possible to extract thetransfer-matrices by a straightforward generalization of the conventional nearest-neighbor case [68]: Suppose thatwe have a L x × L Hamiltonian in a quasi-1D geometry with hopping range p in the x -direction. It can be writtenas H = (cid:80) L x r x =1 (cid:80) pq = − p H r x ,r x + q , where H r x ,r x + q is the Hamiltonian describing hopping from slab r x to slab r x + q . Λ Λ n −1 n −1 (d)(a) (e) (b) (f) (c) 𝑝 𝑝𝐿 𝑥 𝐿 𝑥 Λ Λ ΛΛ n −1 = 8.4 n −1 = 5n −1 = 7.8 n −1 = 11.1 FIG. SF4. (a) the dimensionless Lyapunov exponent Λ = ˜Λ L = L/ξ , Eq. (11) of the main text, as function of n − for L = 28,range cutoff p = 2, and different values of L x (the size in the x -direction of the G − matrices). (b), (c) Λ as function of L x , fortwo specific values of n − from panel (a). (d) Λ as function of n − for L = 28 and for different values of p (the range cutoff).(e),(f) Λ as function of p , for two specific values of n − from panel (d). We also define ψ r x as a length- L vector containing the wavefunction amplitudes associated with the r x th slab. Theeigenvalue equation can be written as: (cid:80) pq = − p H r x ,r x + q ψ r x + q = Eψ r x for each r x , where E is the energy. We canisolate ψ r x + p and arrive to a recursion relation ψ r x + p = ( H r x ,r x + p ) − · (cid:32) Eψ r x − p − (cid:88) q = − p H r x ,r x + q ψ r x + q (cid:33) . (S14)We can then define the transfer matrix T r x (of size 2 pL × pL ) as: ψ r x + p ... ψ r x − p +1 = · · · · · · ( S · · · · · · I · · · I · · · · · · I ψ r x + p − ... ψ r x − p ≡ T r x ψ r x + p − ... ψ r x − p , (S15)where in the first line appear the corresponding components of equation (S14), and I is the L × L identity matrix.Therefore, from each G − matrix we can extract K = L x − c transfer matrices, where c = 7 is a “safety margin”,which was chosen to be larger than p in order to avoid mixing between the first and the last transfer matrices of G − .The effective length would then be L eff = M K . As for the choice of values for p and L x , a priori it seems that thebigger L x and p are, the smaller the resulting error (since the approximation becomes more accurate). While this isindeed the case for L x , for p the situation is more subtle: To use Eq. (S14) we are required to calculate the inverseof H r x ,r x + p , which become exponentially small as p becomes larger. Therefore, p that is too large will lead to largenumerical uncertainties in the inverse matrix.In order to examine the effects of the value L x on the calculation, we first set p = 2 [which is the exact hoppingrange of G − for the case of no disorder in the system Hamiltonian, see Eq. (5) of the main text] and investigatethe dimensionless Lyapunov exponent Λ for different L x , see Fig. SF4(a)–(c). We can see that L x = 105 is alreadyclose to the limiting value. Then, we set L x = 105 and investigate Λ for different values of p , see Fig. SF4(d)-(f). Inaddition, we have performed a calculation of the critical exponent for L x = 105 and several values of p , and verifiedthat p = 4 and p = 5 already give similar results. Based on this information, we chose L x = 105 and p = 5 in thecalculations that are presented in the main text. (c) (d) 𝑝 p r o b a b ili t y d e n s i t y scaled matrix element scaled matrix element p r o b a b ili t y d e n s i t y (b) 𝑞 (a) a v e r a g e c o rr e l a t i o n FIG. SF5. (a) The distributions of absolute values of the matrix elements ( G − ) ij connecting sites along the x -direction,where i = ( r x , r y ) and j = ( r x + p x , r y ) with p x = 0 , ..., p x = 0 is the onsite term), shifted by their respective expectationvalues and normalized by their respective standard deviations. (b) The distributions of the absolute values of the matrixelements of ( G − ) ij connecting sites along the π/ xy -plane, where i = ( r x , r y ) and j = ( r x + p xy , r y + p xy ) with p xy = 0 , ..., p xy = 0 is the onsite term), shifted by their respective expectation values and normalized by their respectivestandard deviations. (c) Semi-log plot of the expectation values and standard deviations of the absolute values of the elements( G − ) ij connecting sites along the x -direction (blue), and connecting sites along the π/ xy -plane (red).(d) Semi-log plot of the correlation C p x ( q ) [Eq. (S16, normalized by its value at q = 0] of the absolute values of matrix elements( G − ) ij connecting sites along the x -direction with p x = 0 , ..., q is the distance between the correlated terms (along eitherthe x - or the y -direction). For the calculations, M = 2000 realizations of G − were generated, each with L x = 105, L y = 28(with periodic boundary conditions), γ in = 0 . γ , µ ∗ = − .