Detrimental effects of disorder in two-dimensional time-reversal invariant topological superconductors

The robustness against local perturbations, as long as the symmetry of the system is preserved, is a distinctive feature of topological quantum states. Magnetic impurities and defects break time-reversal invariance and, consequently, time-reversal invariant (TRI) topological superconductors are fragile against this type of disorder. Non-magnetic impurities, however, preserve time-reversal symmetry and one naively expects a TRI topological superconductor to persist in the presence of non-magnetic impurities. In this work, we study the effect of non-magnetic disorder on a TRI topological superconductor with extended $s$-wave pairing, which can be engineered at the interface of an Fe-based superconductor and a strongly spin-orbit coupled Rashba layer. We model two different types of non-magnetic random disorder and analyze both the bulk density of states and edge state spectrum. Contrary to naive expectations, we find that the disorder strongly affects the topological phase by closing the energy gap, while trivial superconducting phases remain stable and fully gapped. The disorder phase diagram reveals a strong expansion of a nodal phase with increasing disorder. We further show the decay of the helical Majorana edge states in the topological phase and how they eventually disappear with increasing disorder. These results alter our understanding of effects of impurities and disorder on TRI topological phases and may help explain the difficulty of experimental observation of TRI topological superconductors.


I. INTRODUCTION
Much of the recent developments in condensed matter physics have focused on theoretical prediction and experimental verification of topological phases of matter [1][2][3]. Topological superconductivity constitutes a particular exciting case, thanks to the presence of exotic Majorana quasiparticles and their potential for technological applications [4][5][6]. With topological phases being generally robust and protected against local perturbations, such as disorder and decoherence, topological superconductors have become particularly interesting candidates for topological quantum computing [7,8].
Generally, a topological phase is expected to be robust against any perturbation that does not break its classifying symmetries. In particular, a time-reversal invariant (TRI) topological phase would be expected to be robust against any perturbation that does not break time-reversal symmetry. As a consequence, several topological TRI systems have been shown to prevail in the presence of non-magnetic impurities and disorder [29][30][31][32]. Despite this strong disorder protection, TRI topological superconductivity has so far been challenging to observe experimentally.
One promising proposal for experimental realization of a TRI topological superconductor is based on proximity-induced superconductivity from an Fe-based superconductor into a Rashba spin-orbit coupled (SOC) layer [28]. Belonging to class DIII, the effective Hamiltonian for this hybrid structure in two dimensions (2D) is topologically characterized by a Z 2 index [9,10]. In the topological phase, a Kramers pair of helical Majorana edge states propagate along the boundary of the system. By simply tuning the chemical potential it is also possible to tune from the topological to a trivial phase, with an intersecting nodal phase [28]. The effects of single magnetic and non-magnetic impurities have previously been studied in this system [33]. In particular, it was shown that non-magnetic impurities behave very differently in the topological and trivial phases. For any material realization and its functionality it is however more relevant to consider random disorder and its effects.
In this work, we perform numerical calculations to investigate the effect of random non-magnetic disorder on TRI topological superconductors by studying hybrid structures composed of an Fe-based superconductor and an effective 2D material with Rashba SOC. In particular, we consider both Anderson disorder, where the chemical potential is randomly fluctuating, and concentration disorder, where randomly placed dilute but strong non-magnetic scatterers are present. By evaluating the bulk density of states (DOS), we find that even moderately weak non-magnetic disorder induces subgap states that quickly fill the entire energy gap in the topological phase.
Some aspects of an disorder-induced metallic phase has previously been discussed [34,35] and attributed to sign-changing potential fluctuations, while we demonstrate fragility of topological phase for both small disorder and without any need for sign-changing scattering. Furthermore, we are able to trace this disorder sensitivity of the topological phase to the existence of subgap states in the single-and few-impurity limit, although a single impurity never generates states near zero energy, such that only random bulk disorder can fully destroy the gap. We find that the critical disorder strength that fully closes the gap is almost exactly proportional to the size of the topological gap in the clean limit. This leads to the topological phase being particularly fragile for small gaps. We further show that this disorder behavior of the topological phase stand in sharp contrast to that of trivial conventional and extended s-wave superconductors, which both remain fully gapped even for strong disorder. Consequently, by increasing the disorder, we observe a strong expansion of the nodal region in the phase diagram, essentially only at the expense of the topological phase. We also study the evolution of the helical Majorana edge modes with increasing disorder and show how they quickly delocalize into the bulk and are eventually destroyed when the bulk gap closes. These results clearly demonstrate that the naive expectation of robustness against symmetry-preserving perturbations are not accurate for TRI topological superconductors; instead they are surprisingly sensitive to non-magnetic disorder. With non-magnetic disorder being very common in most superconductors and hybrid structures, our results indicates that experimental observation of TRI topological superconductivity might be challenging.
The remainder of this work is organized as follows. In Sec. II, we introduce a 2D lattice model to study the hybrid structure of an Fe-based superconductor and a Rashba SOC layer. In Sec. III, we present our results where we first analyze the disorder dependence of the bulk DOS in Sec. III A. Next, in Sec. III B we discuss the phase diagram in the presence of disorder, and finally in Sec. III C we show the effect of disorder on the topological edge states. We summarize our results in Sec. IV.

II. MODEL
A TRI topological superconductor can be constructed at the interface of a 2D Rashba SOC material and an Fe-based s ± -wave superconductor, as originally proposed by Zhang et al. [28]. In this work, the authors considered a 2D square lattice with nearest-neighbor hopping t and Rashba spinorbit interaction λ R from the 2D Rashba SOC material and with superconducting pairing, consisting of on-site ∆ 0 and isotropic nearest neighbor ∆ 1 pairing terms, induced by the extended s-wave symmetry of the Fe-based superconductor. The tight-binding Hamiltonian for this system within the standard mean-field framework for superconductivity reads: with ∆ ij = ∆ 0 and t ij = µ for i = j, where i = (i x , i y ) represents a site in a 2D square lattice and µ is the chemical potential. We restrict the range of hopping and superconducting pairing to nearest-neighbors, i, j , with t i,j = t and ∆ i,j = ∆ 1 . For clarity, we scale all parameters to the nearest neighbor hopping by setting t = 1. The combination of ∆ 0 and ∆ 1 results in an extended s-wave pairing symmetry that retains the symmetry of the normal-state Hamiltonian but changes sign between the Γ and M points of the first Brillouin zone. The Rashba SOC also splits the normal-state Fermi surface, such that by adjusting the chemical potential one can make the inner and outer Fermi surfaces experience pairing gaps with either same or opposite signs. This has direct consequences for the topology of the system. It is noteworthy that in Hamiltonian Eq. (1), no scattering takes place between the two helical Fermi surfaces. Furthermore, there is no interband scattering due to impurities, as the topological features protected by symmetry prevent any such interband scattering. Time-reversal symmetry is preserved for the Hamiltonian in Eq. (1), which means it belongs to class DIII. For this class in 2D the topological classification is through a Z 2 index and not a Chern number [9,10]. The Z 2 invariant is set by the relative sign of the superconducting gap for each of two SOC split Fermi surfaces. Opposite signs of the gap yields a topological phase with ν = 1 and s ± -wave superconducting pairing, while the same sign gives a trivial phase with ν = 0 and s ++ -wave pairing [36,37]. The topological phase is known to have a Kramers pair of 1D helical Majorana states at any boundary toward a trivial region [28].
In Fig. 1 we present the phase diagram of the Hamiltonian Eq. (1) as a function of the Rashba SOC λ R and chemical potential µ. Depending on the size of λ R and µ compared to the ratio of the superconducting gaps, r = t∆ 0 /∆ 1 , the system has three different phases: topologically trivial (ν = 0), nontrivial (ν = 1), and nodal phase [28]. The topological phase occurs for |µ − r| < 2λ R |r| − r 2 /4 (red in Fig. 1), while the trivial phase requires |µ − r| > 2λ R 2 − r 2 /8 (blue in Fig. 1). The white region in between these two phases represents the nodal phase, where at least one of the SOC split Fermi surfaces crosses the nodal line of the superconducting gap function. Interestingly, the trivial, topological, and nodal phases all meet at the origin µ = r, λ R = 0, forming a triple point in the phase diagram.
To include the effect of non-magnetic, or potential, disorder, we independently implement two different disorder models, namely local random chemical potential fluctuations, or Anderson disorder, and dilute but strong non-magnetic scatterers, or concentration disorder. Anderson disorder is modeled by adding an on-site random chemical potential i to every site, To simulate concentration disorder, we add a constant strong potential to a randomly chosen, small, subset of lattice sites, where Λ * is a small subset of all sites of the square lattice Λ.
The disorder concentration c is then given by the ratio of the Phase diagram of a clean hybrid structure consisting of an Fe-based s ± -wave superconductor and a 2D Rashba SOC layer, described by Eq. (1), as a function of µ − r and λ R . Topologically trivial s ++ -wave (ν = 0), non-trivial s ± -wave (ν = 1), and nodal phases are depicted by blue, red, and white, respectively. The color intensity coding corresponds to the energy gap size. Blue and red solid lines mark the phase boundaries.
dimensions of Λ * and Λ, which becomes a chosen input parameter, while the position of the sites contained in Λ * is random. A notable difference between these two disorder types is that, for Anderson disorder the added on-site terms average to zero, while for concentration disorder, the effective chemical potential is shifted, on average, by cV.
We solve H = H 0 + H (i) dis using exact diagonalization within the TBTK code [38][39][40][41]. We use large lattices with both periodic boundary conditions (PBC) and open boundary conditions (OBC) in order to both assess the effect of disorder on bulk properties and its influence on the topological edge states appearing at the boundary in the topological phase. To be able to compare our results for Eq. (1) with those of a conventional s-wave superconductor, we also perform the same calculations but setting ∆ 1 = 0. This results in modeling a trivial and conventional s-wave superconductor with only onsite pairing, which has no topological features and benefits from strong protection against disorder as established by Anderson's theorem [42].
In this work, for simplicity, we assume ∆ 0 = ∆ 1 = 0.2 and λ R = 0.5. We also choose µ = 1.0 for representing the topological s ± -wave phase and µ = 2.9 for representing the trivial s ++ -wave phase, unless otherwise stated. For these two values of the chemical potential, the size of the excitation gap in the clean case is almost the same, 2∆ ≈ 0.24, and we can thus directly compare the influence of disorder in the subgap regime between these two topologically distinct cases. For the conventional s-wave superconductor we set ∆ 0 = 0.15, which gives rise to an excitation gap similar to the two aforementioned superconducting states. Note that we here keep the superconducting order parameter constant throughout the sample. Technically, a more accurate treatment would require us to calculate ∆ 0 and ∆ 1 in a self-consistent manner. However, previous tests on similar proximity-induced hybrid structures have revealed only minor quantitative corrections to the gap size and hardly any changes to phase diagram [31,43].
Finally, finite size effects are usually unavoidable in solid state simulations. The inclusion of random disorder makes it even more challenging to obtain configuration-independent results. Throughout this work we report exact diagonalization results using a square lattice with 51 × 51 sites (unless otherwise stated), but we have checked our results also using smaller lattice sizes. Moreover, when we perform disorder averaging we choose the number of disorder realizations n such that σ/ √ n < 0.015 is fulfilled, where σ is the standard deviation. Typical values here are n = 40 and σ = 0.02. This procedure guarantees that all our disorder averaged results are sufficiently independent of a specific disorder realization. Furthermore, in order to calculate the bulk DOS we use a Gaussian kernel for smoothing the data, setting the standard deviation for the Gaussian kernel to Γ = 0.02.

III. RESULTS
Having established the phase diagram and the overall properties of our system in the clean limit we now turn to the effects of non-magnetic disorder in this TRI topological superconductor. First, we study the properties of the bulk DOS in the presence of different disorder types and strengths. Despite naive expectations of disorder robustness for topological systems as long as their symmetries are preserved, we show that the TRI topological superconducting phase is actually very fragile against non-magnetic disorder. This leads us to derive a very different phase diagram at finite disorder, where especially the nodal phase is heavily expanded. Finally, we explore the impact of disorder on the helical Majorana edge states where we capture their delocalization and eventual destruction with increasing disorder.

A. Bulk density of states
To illustrate the impact of non-magnetic disorder on the bulk properties of a TRI topological superconductor, we first evaluate the bulk DOS for the Hamiltonian Eq. (1) in the presence of both Anderson and concentration disorder, H (i) dis (i = 1, 2), with PBCs imposed. In Fig. 2 we show the DOS for a single disorder configuration for three different disorder strengths W (top) and disorder concentrations c (bottom). Each panel in Fig. 2 shows the DOS for three different superconducting states with topological s ± -wave (red), trivial s ++wave (blue), and trivial conventional onsite s-wave (green) symmetries. In particular, in Figs. 2 (a-c), we consider Anderson disorder using W = 1.0, 1.8 and 2.4, respectively. Initially, in the absence of the disorder and for only small disorder, the size of the energy gap is equal for all three superconducting states, see Fig. 2 (a). However, when increasing the strength of the Anderson disorder, we find the DOS to be different between the topological and trivial phases in a number of respects. Most importantly, the gap size in the topological s ± -wave superconducting phase is strongly reduced. Surprisingly, even for W ≈ 1.8, which is still small in comparison to Next, we investigate concentration disorder in Fig. 2 (d-f). Here, we set the impurity scattering strength to V = 2 and vary the disorder concentration between 1%, 4%, and 7%, respectively. Again, we observe that very small concentrations of disorder strongly affect the topological phase. As seen in Fig. 2(e), an impurity concentration of only 4% fully closes the energy gap of the topological phase, producing a gapless state. In contrast, the conventional s-wave state remains completely unaffected even for large concentrations. Also the trivial s ++ -wave state survives large impurity concentrations, albeit with a slightly reduced gap size. To conclude, the results in Fig. 2 show that for both Anderson and concentration disorder the topological phase is very fragile against disorder, in sharp contrast to the robustness of trivial s-wave phases. Given that the disorder does not break any of the symmetries protecting the topological phase, this is an unexpected finding.
We can start to understand the strong effect of non-magnetic disorder on the topological phase and its gap by first considering the single-impurity case. As shown in Ref. [33], a single potential impurity added to the TRI superconductor in Eq. (1) produces a single pair of subgap states (symmetric around zero energy and two-fold degenerate due to timereversal symmetry) in the energy spectrum in the topological phase, while the trivial phase never hosts any subgap states for non-magnetic impurities. However, a single impurity was also shown to never induce states at or near zero energy. In our results, we clearly see that we generate a finite bulk DOS also at zero energy for a finite concentration of disorder.
To connect the earlier single impurity results [33] to our bulk results, we next analyze the energy spectrum when we add N randomly positioned impurities, letting N = 1, . . . , 20. The results are summarized in Fig. 3 where we plot the subgap spectrum as a function of the impurity strength V. The case of a single impurity N = 1 is shown in Fig. 3 (a), where we observe subgap state clearly forming for V 1.5 but notably always staying rather far from zero energy. With increasing number of impurities, N, we observe an increasing number of subgap states that also move closer to zero energy. In Fig. 3 (b,c) we plot the cases of N = 10 and N = 20, respectively. Already for N = 20, which corresponds to an impurity concentration of only c = 0.008, we observe a considerable contraction of the energy gap. From the plotted sequence of disorder sites, N = 1, 10, 20, we can easily extrapolate to the results for larger concentrations obtained in Figs. 2(a-f) using V = 2, marked by a vertical dashed blue line in Fig. 3. In this way we obtain a smooth connection between the singleimpurity results and bulk results at finite disorder concentration, noting that the difference is that multiple impurities pushes the subgap impurity states closer and closer to zero energy, eventually filling the whole energy gap.
To contrast the results for the topological phase, we also plot results for N = 10, 20 impurities in the trivial s ++ -wave phase in Fig. 3 (e,f). For a few impurities, we hardly observe any subgap states, only for N ∼ 10 do subgap states start to slowly leak into the bulk gap, but they always remain extremely close to the bulk band edge. This explains the slightly reduced bulk gap in Figs. 2(a-f) in the trivial phase. Finally, in Fig. 3 (d) we consider the special case of an "impurity dimer", i.e., two impurities next to each other and of identical strength.
The dimer generally produces lower-lying impurity subgap states than the single impurity, but the exact values vary depending on the dimer orientation on the lattice. In fact, horizontally bonded dimers even achieve zero-energy states at certain disorder strengths. This further supports the fact that for random finite disorder concentration, zero and near zero energy states will always be present. If we instead redo the same calculation as in Fig. 3 but use a disorder strength that on each disorder site varies between [−V, V], we arrive at qualitatively the same results, which eliminates any possible influence of the shift of the overall chemical potential. We also note that all plots shown in Fig. 3 represent only a single disorder realization, but we do not expect disorder-averaging over multiple realizations to qualitatively change the results.
To shed further light on the effects of disorder we investigate more closely when the topological gap is lost. To this end, we introduce W c , the ensemble-averaged disorder strength for which the transition from gapped topological to gapless phase takes place for Anderson disorder. To be computationally achievable, we choose to evaluate W c for points in the topological phase along the dashed line in the phase diagram in Fig. 1 and for λ R ∈ [0, 1], which is the most physically relevant range. In the clean system, the phase diagram is geometrically symmetric along µ = r and the gap size in the topolog- ical phase, ∆ Topo , grows with increasing λ R as shown by the blue curve in Fig. 4 (a). By plotting W c as a function of λ R in red and in the same panel, we find that W c has almost exactly the same dependence as the gap ∆ Topo on λ R . In more detail, we observe an algebraic relation, W c ∝ λ 1/2 R , as indicated by the fit (solid red line) in Fig. 4(a). This strongly suggests that ∆ Topo is the relevant energy scale for the transition to the gapless phase. Hence, the larger the topological gap in the clean limit, the more robust is the topological phase, as also intuitively expected. This means that in regions close to the nodal phase in the clean limit in Fig. 1, even very weak disorder will destroy the topological gap and render a nodal spectrum.
We can alternatively evaluate the critical density c * , for which concentration disorder completely fills the topological gap and generates a nodal phase. In Fig. 4(b) we plot c * as a function of λ R using several different values of V. We find again that stronger disorder, here represented by higher disorder densities, are needed to reach the critical density for larger λ R , or, equivalently, a larger topological gap ∆ Topo in the clean limit. Instead keeping λ R constant, for larger V, we find smaller c * , since stronger scatterers pushes the impurity subgap states further into the gap as seen in Fig. 3, which facilities the transition to the nodal phase.

B. Phase diagram
Having established that the energy gap in the topological phase is quickly diminished and even disappears in the presence of disorder, we next establish the phase diagrams for finite disorder. By using extensive numerical simulations extracting both the topology and gap size, we produce in Fig. 5(a-c) the phase diagrams analogous to the clean system in Fig. 1 for different strengths W of Anderson disorder, averaged over 50 independent disorder configurations. As expected, stronger disorder affects the phase diagram drastically. First of all, the nodal region (white color) grows substantially and especially occupies regions that were topological in the clean limit. Second, the trivial phase is much less affected by disorder than the topological phase. Additionally, while the nodal phase shrinks in the clean limit into a triple point at µ = r and λ R = 0, in the presence of disorder we find the nodal phase expanding over a large area for small λ R . We also note that disorder affects the phase diagram in an asymmetric way with respect to µ − r. This is not surprising as the gap size in the clean limit is not symmetric with respect to µ − r, see color intensity in Fig. 1. As a consequence, the filling of the gap through disorder-induced states occurs differently for positive and negative values of µ − r.
To further illustrate the effect of disorder on the phase diagram, we plot in Fig. 6 phase diagrams for fixed Rashba SOC λ R , while we tune the chemical potential and disorder strength. As shown in Fig. 6(a), there is no topological phase for λ R = 0. For moderate disorder strength, the trivial phase also remains gapped. By further increasing the disorder strength, the gap is however reduced and eventually, for disorder of the order of the bandgap, we lose the gap even in the trivial phase. Next setting λ R = 0.5 in Fig. 6(b), the clean system is topological for |µ − r| 1. By adding disorder, we clearly see how the topological phase is much more fragile to disorder than the trivial phase. Finally, in Fig. 6(c), we set λ R = 1.0 which in the clean limit generates a topological phase for the whole range of |µ − r| 2. Again, we observe that disorder leads to strong suppression of the gap in the topological phase, but due to the initially larger topological gap thanks to the larger λ R , the transition to the nodal phase requires a larger W c compared to Fig. 6(b). Somewhat surprisingly, we find that further increasing the disorder strength W beyond W c in Fig. 6 (b-c) eventually drives the system from the gapless nodal phase into the trivially gapped phase for a range of µ − r. By using OBC and calculating the local DOS we have verified that there are no edge states in this regime, consistent with the trivial topology. This behavior of a re-entrant gapped phase is reminiscent of disorder-induced topological phase transitions found earlier in superconductors, where Anderson disorder has been shown to lead to topological phase transitions between phases with different Chern numbers [44].

C. Edge states
To further investigate the detrimental effect of disorder on topological s ± -wave TRI superconductors, we study the disorder impact on the helical Majorana edge modes by performing OBC calculations for a disc geometry with radius r = 47 lattice points. Such a disc size is sufficient to avoid overlap- ping of the edge wave functions at opposite sides of the disc, which would lead to their hybridization. By tuning the Anderson disorder strength W, we observe a clear destructive influence of disorder on the edge modes. As an illustration, we plot in Fig. 7 the accumulated local DOS for subgap energies |E| < ∆ Topo /4 using a specific disorder configuration. Fig. 7 (a-c) and (d-f) contain identical data but plotted in 3D and 2D, respectively. In Fig. 7 (a,d), in the absence of disorder, we find well-localized helical Majorana edge modes and a fully gapped interior. By introducing moderately small Anderson disorder, as shown in Fig. 7 (b,e), the helical Majorana edge modes already begin to fade away and near-zero subgap states begin to be visible throughout the disc interior. When the disorder strength is increased further, see Fig. 7 (c,f), the edge modes are strongly suppressed and a clear subgap DOS develops in the interior. Similar calculations for concentration disorder confirm that the helical Majorana edge modes in the topological phase disappear with increasing impurity concentration. As a complement, we also plot line cuts through the disc in Fig. 7 (g). Here it becomes clear how the delocalization of the edge states goes hand in hand with the emergence of subgap states in the bulk. This result shows how the loss of the topological gap in the bulk directly affects the edge states, which fade into the bulk and eventually disappear with increasing disorder. In other words, the edge of the topological phase is affected by disorder in the same way as the bulk, a manifestation of the bulk-boundary correspondence.

IV. SUMMARY
In this work, we study the effects of non-magnetic disorder in a time-reversal invariant superconductor, which can be realized in a hybrid structure consisting of an Fe-based superconductor and a spin-orbit coupled Rashba layer. The hybrid structure hosts three different phases, a topological s ± -wave and a trivial s ++ -wave, with an intersecting nodal phase. We perform intensive numerical calculations in the presence of two different types of disorder, Anderson disorder, i.e., random chemical potential fluctuations, and concentration disorder, i.e., random dilute strong potential scatterers, both preserving time-reversal symmetry. The naive expectation is that the topological phase should be robust against this disorder as it preserves all symmetries protecting the topology. Instead we find that the disorder leads to subgap states quickly accumulating in the gap of the topological phase, leading to the closure of the topological gap even for moderately weak disorder. These findings are in contrast to the topologically trivial phase, as well as conventional s-wave superconductors, which are both exceptionally robust to disorder with no notable subgap states even for strong disorder. We are able to trace this disorder fragility of the topological phase to the existence of subgap states in the single-and few-impurity limit. Notably though, single impurities do not generate (near) zero-energy states, while random disorder fully closes the gap, generating a gapless phase. We derive the phase diagram in the presence of disorder, where we very generally observe the expansion of the nodal phase at the expense of the topological phase. In addition to a bulk analysis, we also investigate the helical Majorana edge modes associated with the topological phase. We find that even moderately weak disorder causes the edge states to decay into the bulk and eventually to disappear. Thus, we find that disorder affects both edge and bulk of the topological phase in the same detrimental way. Our results suggest that the lack of a hard energy gap and the extension of the nodal phase due to disorder can easily prevent experimental detection of time-reversal invariant topological superconducting phases.