Scaling of critical wavefunctions at topological Anderson transitions in 1D

Topological Anderson transitions, which are direct phase transitions between topologically distinct Anderson localised phases, allow for criticality in 1D disordered systems. We analyse the statistical properties of an emsemble of critical wavefunctions at such transitions. We find that the local moments are strongly inhomogeneous, with significant amplification towards the edges of the system. In particular, we obtain an analytic expression for the spatial profile of the local moments which is valid at all topological Anderson transitions in 1D, as we verify by direct comparison with numerical simulations of various lattice models.

Topological Anderson transitions, which are direct phase transitions between topologically distinct Anderson localised phases, allow for criticality in 1D disordered systems. We analyse the statistical properties of an emsemble of critical wavefunctions at such transitions. We find that the local moments are strongly inhomogeneous, with significant amplification towards the edges of the system. In particular, we obtain an analytic expression for the spatial profile of the local moments which is valid at all topological Anderson transitions in 1D, as we verify by direct comparison with numerical simulations of various lattice models. The interplay of topology and disorder in electronic systems can give rise to critical transitions between localised phases [1][2][3][4][5][6]. This is of particular interest in 1D, as here there are good reasons to expect a system of electrons to realise an insulating phase, even in the absence of electronic correlations. On the one hand, a 1D chain is generically unstable to dimerisation via Peierls distortion, while on the other, any amount of extensive disorder leads to Anderson localisation. Remarkably however, these two effects can be tuned in combination to induce critical behaviour. This occurs as distinct dimerisation patterns realise distinct topological phases, whereas strong disorder washes out the detail of the dimerisation, which drives the system between the phases. In general, such critical transitions occur in 1D for disordered systems which realise symmetry-protected topological phases, and are described by two-parameter scaling [6]. We shall refer to the transitions as topological Anderson transitions, and it is the resulting critical states that are of interest to us in this paper.
Whether a system may realise distinct topological phases is dictated by both its dimension and the antiunitary symmetries of its Hamiltonian [7][8][9]. In particular, for each dimension a certain five of the ten symmetry classes of the Altland-Zirnbauer classification [10] allow for non-trivial topological phases. In 1D these are the three chiral classes (AIII, BDI, CII) and two of the Bogoliubov-de Gennes classes (D, DIII). The situation in 1D is however more special as here, irrespective of symmetry class, all topological Anderson transitions share the same critical properties, a phenomenon dubbed superuniversality [4]. This was established by identifying exact mappings between representative single-channel models for each of the five classes, and by showing that for a long chain the properties of a system with many channels are determined by a dominant effective single channel.
It is of great interest to understand the properties of the critical wavefuntion at this superuniversal transition. This has been partially investigated by Balents and Fisher [11] who calculated the system-size scaling of the moments of the normalised wavefunction in the bulk and their spatial correlations where . . . denotes an averaging over a disorder ensemble, N is the system size and q > 0. The 1/N -scaling of the moments is characteristic of localisation, while the algebraic decay of the spatial correlations reflects the critical nature of the wavefunction [12].
Here we analytically show that the local moments are in fact inhomogeneous throughout the system. In the bulk they scale as and directly at the edges they scale as 1/ √ N . Equation (2) describes how the 1/N -scaling in the bulk changes to 1/ √ N towards the edges. The strong amplification of the disorder-averaged wavefunction near the boundaries |ψ(1)| 2q / |ψ(N/2)| 2q ∝ √ N distinguishes it significantly from all other critical wavefunctions, and can be understood as a remnant of the topological edge mode that is created/destroyed by the transition.
These results are derived from a specific single-channel model in class BDI for which the critical wavefunction at the transition is exact and analytically tractable. Based on superuniversality, we argue that the results are generally valid at all 1D topological Anderson transitions. We verify this by directly comparing them with numerics for both a class BDI model with no E = 0 wavefunctions away from the N → ∞ limit, and for a multi-channel model belonging to symmetry class D.
To explore the critical physics of 1D topological Anderson transitions we focus on a disordered dimerized tightbinding chain. This model appears as the static limit of the Su-Schrieffer-Heeger model [13,14], which is used for example to describe polymers such as polyacetylene. For reasons that will be made clear shortly, we primarily focus on a chain of odd length 2N + 1, and label the sites from 0 to 2N : Disorder appears through the hopping amplitudes where ω j , ω j are independent random variables uniformly distributed on the interval [−0.5, 0.5]. The model is invariant under the time reversal symmetry T = K, where K is complex conjugation, and has a chiral symmetry The parity operator P = T S squares to the identity when acting on single-particle states, and so the model belongs to the BDI symmetry class.
The odd length of the chain guarantees the existence of a zero-energy eigenstate for any set of parameters t j , t j [15,16]. This follows from the chiral symmetry of the model which restricts the spectrum to be symmetric about zero. The Schrödinger equation at zero energy reduces to a simple recursion relation with solutioñ where the ∼ indicates that this wavefunction is unnormalised. This eigenstate describes a topological boundary mode that is localised at one edge of the chain. For the clean case (t k = t, t k = t ) the mode resides at the left edge if |t| < |t | and at the right edge if |t| > |t |, while precisely at |t| = |t | there is a topological transition at which the wavefunction is uniform. In the presence of disorder this transition becomes a topological Anderson transition, and it can be driven purely by a change in the disorder strength.
It is worthwhile to contrast this behaviour of the edge mode with that for an even length chain, for which there are boundary modes at both edges of the chain on one side of the transition, and no boundary modes on the other side. The connection can be made by considering adding an extra site to the odd length chain. On its own the extra site has a single localised zero energy mode. On coupling it to an end of the chain this localised mode will either remain (protected by the chiral symmetry) if there is no zero energy state at that edge, or it will combine with and annihilate the boundary mode if there is one present. More formally, for the even length chain the two sides of the transition can be distinguished by the Zak-Berry phase [17,18], which is a Z 2 bulk invariant, and the extra site (which breaks the 2-site unit cell) can be thought of as giving rise to a defect mode [19]. We remark that for the odd length chain there are two choices for the unit cell, and so there is no natural way to call one phase trivial and the other topological.
We now calculate the local moments per unit cell of the normalized zero-energy wavefunction: where n indexes the unit cells of the underlying lattice, and Λ n is the set of sites in the n th unit cell.
In the present case I q (n) = |ψ 2n | 2q , with ψ 2n ≡ ψ 2n /( N m=0 |ψ 2m | 2 ) 1/2 , as the wavefunction has support on only the even sites of the chain. The normalisation factor contains non-trivial dependence on the disorder parameters ω k , ω k , and to handle this we use a technique developed for continuum models that maps the calculation to one in Liouville quantum mechanics [11,20]. In the present case, working explicitly with a lattice model allows for a controlled description of the problem in the strong disorder regime. We outline here only the main steps of the calculation and defer the details to the Appendix.
The disorder averaged local moments of the wavefunction I q (n) can be expressed as a functional integral over random walk configurations of the variables z k = 2 ln |t k /t k |. In the thermodynamic limit N → ∞, the behavior of this integral is determined by the first two moments of the distribution function of the z k . As a result, evaluating I q (n) reduces to solving a discrete-time version of Liouville quantum mechanics -the quantum mechanics of a particle moving in an exponential potential. This latter problem can be solved by approximating the exponential potential by a hard wall.
As shown in the Appendix, the averaged local moments can be expressed in terms of a propagator P (x, x , m) as where A q is a normalisation constant and V (x) = e −2q √ ax with a = ( z 2 k − z k 2 )/2. Within the hard wall approximation, the propagator P (x, x , m) is given by a solution of the classical drift-diffusion equation with an absorbing boundary condition at x = 0: which satisfies the initial condition P (x, √ a corresponds to the insulating phases, which are distinguished by the sign of γ, while zero drift corresponds to criticality, and thus z k = 0 is the condition on the distribution of t j , t j that the wavefunction of Eq. (5) is critical.
We briefly summarise the behaviour of the moments away from the critical regime. Evaluating the integrals in Eq. (7) for γ = 0, we obtain where C = 1 − e −2 √ aγ , Γ is the gamma function, and we assume that γn √ n 1 for γ > 0, and |γ|(N − n) √ N − n 1 for γ < 0. Therefore the wavefunction is exponentially localized at the left edge for γ > 0 and at the right edge for γ < 0. Thus the two topologically distinct insulating phases in the model can be distinguished by probing the local statistics of the wavefunctions. In terms of the underlying diffusion process, the spatial asymmetry of the wavefunctions can be interpreted through the interplay of the non-zero drift with the absorbing boundary condition at x = 0, see the Appendix.
For critical wavefunctions, corresponding to γ = 0, integrating Eq. (7) yields The scaling of the moments is I q (n) ∝ 1/N deep in the bulk of the system, where n ∝ N , which is the scaling law typical for a localized wavefunction [12]. At the same time, the moments get strongly amplified approaching the edges of the system. Exactly at the edge, the moments are given by The 1/N τ -scaling of I q (0) with a non-integer exponent τ = 1/2 is the scaling law specific for fractal wavefunctions [12]. The fact that τ < 1 is however unusual, as typically amplitudes of the critical wavefunctions are suppressed rather than enhanced near the boundaries and τ > 1 [21]. The inhomogeneous nature of the local moments and the amplification of the critical wavefunction at the edges can be traced to the presence of the localized boundary modes in the insulating phase. In terms of the classical diffusion process it can be traced to the role of the potential V (x), which gives an increased weight to trajectories passing near the absorbing boundary condition at x = 0. Viewed in this way it is natural that The spatial correlations of the critical wavefunctions can be characterised by the two-point correlation function C mn = |ψ 2m ψ 2n | q . Using the same approach as above, we obtain that C mn ∝ |m − n| −3/2 N −1 in the bulk of the system, and C 0n ∝ n −3/2 N −1/2 at the edge. The power-law scaling of the correlation function is another manifestation of its fractal nature [11].
To investigate the range of validity of the above results, which are derived in the limit of large N , we analyse the critical wavefunction numerically. For this we examine a critical point of the Hamiltonian of Eq. (3) given by t = 1, t = 0, W = 2, W = 4, and compute the local moments by averaging the zero-energy wavefunction over 10 6 −10 8 disorder realisations. The system-size scaling of the local moments, both directly at the edge and in the centre of the chain, are presented in Fig. 1. It is seen that the edge scaling is more robust than the bulk scaling for  small system sizes, while for increasing N an excellent fit to Eq. (11) and (10), for the edge and bulk respectively, emerges.
In Fig. 2 the spatial profile of the q = 1/4, 1/2, 1, 2, 4 local moments for N = 50, 100, 200 are plotted, and compared with a fit to Eq. (10). It is seen that for q > 1 the formula captures the moments right out to the edge,  3), with the disorder average taken over wavefunctions whose energy lies within a window about 0, as described in the text. Plot (i) shows the system-size scaling of the edge moments for q = 1/4, 1/2, 1, 2, 4 and the dashed lines give fits to 1/ √ N scaling. The data is terminated at the black squares on the q = 1, 2, 4 lines as beyond these system sizes numerical inaccuracy arises due to machine precision. In (ii) the spatial profile of the q = 1 local moments are plotted for N = 25, 50, 100, 200 and the dashed lines show fits of the data to Eq. (10) where the normalisation is set by fixing the value at the centre of the chain.
while for q < 1 there are deviations near the edges of the chain which get smaller as the system size is increased. Since the moments become uniform in the limit q → 0 one would need diverging system sizes in order to see the predicted form in this regime.
We now perform checks on the generality of our analytic expressions, to support the claim that they are valid at all topological Anderson transition in 1D. First we consider a dimerised tight-binding chain with an even number of sites, as then there is no wavefunction exactly at E = 0 for any finite length. To examine the critical properties in this case we take the disorder average . . . to run over all wavefunctions in the disorder ensemble whose energy is within a window [−E typ , E typ ], where E typ = exp( ln(min α |E α |) ) with {E α } the set of eigenenergies for a given realisation. The even length chain is realised by setting t 1 = 0 in Eq. (3), which decouples the site 0, and it is analysed using exact diagonalisation. In Fig. 3(i) the system-size scaling of the local moments at the edge of the chain are plotted and compared with the predicted 1/ √ N scaling, while in Fig. 3(ii) the spatial profile of the q = 1 local moments are presented and compared with the prediction of Eq. (10). In both cases a very good fit is seen, which indicates that the results are valid in general and not just to exactly critical wavefunctions. The second check on the generality of our results that we perform is to examine a multi-channel model in the Bogoliubov-de Gennes symmetry class D. Specifically we consider the following Hamiltonian where t D is a hopping parameter, C j is a column vector {c j,1 , c j,2 , . . . , c j,N ch } of N ch fermions at site j, B j is an N ch × N ch random matrix satisfying B j = −B T j , and N ch is the number of channels. In practice we set B j = A j − A T j where A j are real matrices whose elements are independent random variables uniformly distributed on the interval [−0.5, 0.5]. The Hamiltonian belongs to symmetry class D when N ch ≥ 2, as complex conjugation K provides an anti-unitary symmetry that anti-commutes with the Hamiltonian, whereas in the absence of a sublattice symmetry there is generically no anti-unitary symmetry that commutes with the Hamiltonian. The model sits exactly at a topological Anderson transition for any t D , and can be taken into a topological phase by adding a dimerisation to the hopping. Indeed, to view the topological properties of the model one must regard it as having a two-site unit cell, and it is for this reason that we set the sum to run to 2N in Eq. (12).
Again we numerically examine the disorder averaged local moments of the critical wavefunction. For this we take N ch = 5 and obtain the critical wavefunctions through exact diagonalisation. In Fig. 4 we plot again the system-size scaling of the local moments at the edge of the chain and the spatial profile of the q = 1 local moments, as were presented in Fig. 3. We find once more an excellent fit, and view this as confirmation that our results are generally applicable at topological Anderson transition in 1D.
In summary, from the known exact zero-energy wavefunction of the Su-Schrieffer-Heeger model we have derived analytical expressions for the system-size scaling of the local moments, which are valid at all topological Anderson transitions in 1D. The generality of the expressions follow from the superuniversality of these transitions, and is also evidenced by our numerical simulations of other models, which show excellent agreement despite having distinct features from the original model. The results show that the moments are inhomogneous, with strong amplification towards the edges of the system, and this distinguishes critical wavefunctions at topological Anderson transitions in 1D from all other critical wavefunctions studied so far. It is an interesting open problem whether similar features of critical wavefunctions can be found in higher dimensions.

Moments of the normalized wavefunctions a. Moments of the normalized wavefunctions and classical diffusion
The explicit form of the wavefunction given by Eq. (5) allows us to straightforwardly calculate its statistical properties, if the normalisation condition is not taken into account. However, for the normalized wavefunction the problem becomes much more complicated, as the normalisation constant depends on all variables t k , t k in a non-trivial way. In order to overcome this problem we map the calculation to one in Liouville quantum mechanics [A1, A2].
First, we define the random variables x n = n k=1 t k t k 2 and we set x 0 = 1 . Then the disorder averaged local moments of the normalized wavefunction I q (n) = |ψ 2n | 2q are given by I q (n) = where Γ is the gamma function. With the help of the above identity we obtain where J(α) = x q n e −α N p=1 xp and we used the fact that x 0 = 1. The disorder averaging can be performed using the statistical independence of the random variable z k = ln t k t k 2 : where P z (z k ) is the distribution function of z k . The structure of the above integrand suggests to introduce new variables in terms of which J(α) takes the form where t 0 ≡ 0, g p (t, α) = e −αe t for p = n and g n (t, α) = e qt−αe t . Using the expression for P z (z) in terms of its Fourier transform F (s) = ∞ −∞ dze isz P z (z), we obtain One can notice that hence the main contribution to the integral in the limit N → ∞ comes from s p → 0. Therefore we can use the Taylor expansion of F (s) about s = 0. Keeping only the terms up to the quadratic order we can write F (s) ≈ e −ibs−as 2 , where b = − z k , a = ( z 2 k − z k 2 )/2. Then the integrals over s p become Gaussian and can be evaluated explicitly: This expression can be viewed as a discrete in time path integral for a particle moving in an exponential potential -Liouville quantum mechanics. One can approximate such a potential by a hard wall potential: where θ(t) is the Heaviside step function, αe z(α) = C and C is a constant of order 1. Within this approximation J(α) is given by where r p = (z − t p )/2 √ a, r 0 = z/2 √ a and γ = b/2 √ a. In order to calculate the above integral, we introduce an integral operatorK defined as The propagator corresponding to the path integral (A10) is related toK as If P (x, x , m) is known, then J(α) can be found as where V (x) = e −2q √ ax , l 1 = N − n, l 2 = n. In the hard wall approximation P (x, x , m) = 0 for x < 0, and hence the integration over α in Eq. (A2) must be restricted to α < C. Taking into account that α q−1 e qz(α) = C q /α, the expression for the moments reads Changing the integration variable α by z(α) and noting that e −Ce −z ≈ 1 in the hard wall approximation for z > 0, we obtain where we finally changed z by x 3 = r 0 = z/2 √ a. This equation is equivalent to Eq. (7) in the main text with A q = 2 √ aC q /Γ(q). In order to find the propagator P (x, x , m), we notice that, as follows from its definition, it satisfies an integral equation Assuming that P (y, x , m) is a slowly varying function compared to e −(x−y+γ) 2 (this assumption is justified by the solution that we obtain, so is self-consistent), one can transform the above integral equation into the drift-diffusion equation with the absorbing boundary condition at x = 0. The solution of the drift-diffusion equation obtained by the standard eigenfunction method reads One can check by direct calculation that the above expression for P (x, x , m) solves the integral equation (A16) provided that m, x 1, i.e. it works well for large times and far away from the boundary x = 0. Eq. (A15) and Eq. (A17) give a complete description of the local moments of the wavefunctions in terms of classical diffusion.

b. Localized wavefunctions
The insulating phase is characterized by a non-zero drift constant γ = 0. Substituting Eq. (A17) into Eq. (A15) and integrating over x 1 and x 3 we obtain is the error function and erfc (x) = 1 − erf (x). To evaluate this integral we may use the following approximation for erf (x) which is justified in the above integral provided that γl i √ l i 1. This approximation yields the following result for I q (n) The positive γ case corresponds to the wave function localized at the left edge, while the negative γ case corresponds to the wave function localized at the right edge. Thus the two different topological phases in our model can be distinguished by probing the local statistics of the wavefunctions.
In order to determine C we use the normalization condition N n=0 I 1 (n) = 1, which gives The spatial asymmetry of the result can be naturally understood in terms of the underlying diffusion process. Let us consider for simplicity γ > 0 and two sub-cases n = 0 and n = N . For positive γ, the propagator P (x, x , m) describes trajectories with a drift to the left (the mean value of x decreases in "time" m), as it follows from Eq. (A17). Now consider the wavefunction at n = 0, the moments are given by The propagator P (x 2 , x 1 , N ) describes trajectories of classical particles, which start from x = x 1 at time m = 0 and propagate to x = x 2 at m = N . The contribution of each trajectory is weighted with an exponential function V (x 2 ) = e −2q √ ax2 at the final point. For this reason the main contribution to the results comes from the trajectories with x 2 close to zero. Such trajectories can be found and, due to the drift term, the most probable starting point for them is x 1 = γN . The contribution of such trajectories yields the result of order of one.
If we now consider the opposite edge n = N , then The weight function V (x 1 ) = e −2q √ ax1 is associated with the starting point, which forces the initial coordinate x 1 to be close to zero. On the other hand, the absorbing boundary condition at x = 0, together with the presence of the drift term which tends to decrease the position of a particle in time, lead to the exponential small survival probability of such trajectories and to the exponentially small final result.

c. Critical wavefunctions
The critical phase is characterised by a vanishing drift constant γ = 0. In this case integration over x 3 and x 1 in Eq. (A15) leads to the following result: Recalling that l 1 = N − n, l 2 = n and assuming that the site n is in the bulk of the system, we are allowed to replace two error functions in the thermodynamic limit N → ∞ by their asymptotic expressions for l 1 , Then the remaining integration over x yields In order to fix the value of the constant C we use again the normalization condition N n=0 I 1 (n) = 1, which gives C = a/2.
If n = βN with N -independent constant β, then I q (n) ∝ 1/N , which is the scaling law typical for a localized wavefunction.
However, if we are interested in the scaling of the local moments on a site which is close to the boundary, the result will be very different. Indeed, let us consider n = 0. Then using Eq. (A22) and the same steps as above, we obtain Thus I q (0) ∝ 1/ √ N , which is the scaling law typical for a fractal wave function. The same scaling is valid for the opposite edge n = N .
The two different results for the scaling of the moments at the edge and in the bulk of the system can be traced to the presence of the localized boundary modes in the insulating phase. It can be also understood in terms of the classical diffusion process as follows. The propagator P (x 2 , x 1 , N ) describes trajectories of classical particles, which start from x = x 1 at time m = 0 and propagate to x = x 2 at m = N . In Eq. (A22) the contribution of each trajectory is weighted with an exponential function V (x 2 ) = e −2q √ ax2 at the final point. For this reason the main contribution to the result comes from the trajectories with x 2 close to zero. Thus Eq. (A22) gives a probability of getting from any starting point x 1 to an x 2 in the vicinity of zero in time N , and according to our calculations this scales as 1/ √ N . Similarly, Eq. (A15) for n = N/2 counts the contributions from the trajectories starting at any point x 1 reaching an x 2 in the vicinity of zero in time N/2, and then reaching an arbitrary final point x 3 in time N/2. The probability for the first part of the process scales as 1/ √ N , but the same is true for the second part of the process due to the time reversal invariance of the problem. Hence the total probability must scale as (1/ √ N ) 2 = 1/N .

Two-point correlation function
Spatial correlations of the wavefunctions can be determined through the two-point correlation function Performing all the steps which led us to Eq. (A15) we obtain where U (x) = e −q √ ax , l 1 = N − n, l 2 = n − m, l 3 = m. Integration over x 1 and x 4 , and expansion of the resulting error functions gives The above integral can be computed explicitly in terms of the error functions. However, in order to find how it scales with |m − n|, it is more convenient to write it in the momentum space using the complete set of standing waves [A2] z|k = 2 π sin(kz), k > 0.
and substituting two equations above into Eq. (A30), we find Evaluating the integral over x and scaling the variable k → q √ ak we obtain the final result C mn = 64C q π 2 aq 3 Γ(q) √ l 1 l 3 ∞ 0 dk k 2 e −l2q 2 ak 2 /4 (k 2 + 1) 4 , which coincides up to a constant with the result by Balents and Fisher [A2]. For l 2 = 0, the integral over k yields π/32 [A3] and hence we reproduce the expression for the moments (A26). For l 2 1 and m and n of order of N 1 we find showing that the critical wavefunctions exhibit power law correlations in the bulk of the system. If one of the points is located at the edge of the system m = 0, then the above derivation should be modified. Eq. (A29) must be replaced by where l 1 = N − n, l 2 = n. Performing all the steps similar to the ones, which led to Eq. (A34), we find Computing the integrals over x 2 and x 3 we arrive at the final result C 0n = 16C q π 3/2 √ aq 2 Γ(q) √ l 1 ∞ 0 dk k 2 e −l2q 2 ak 2 /4 (k 2 + 1) 3 .
For l 2 = 0, the integral over k gives π/16 and we reproduce Eq. (A27). For l 2 1 and N − n of order of N we find which has the same scaling with m as in the bulk, but a different scaling with N . Finally, for the both points located at the opposite edges an analogue of Eq. (A29) reads which can be transformed to the momentum representation as Calculating the integral over x, we obtain C 0N = 4C q πqΓ(q) ∞ 0 dk k 2 e −N q 2 ak 2 /4 (k 2 + 1) 2 , which implies that C 0N ∝ N −3/2 .