Mesoscopic transport signatures of disorder-induced non-Hermitian phases

We investigate the impact on basic quantum transport properties of disorder-induced exceptional points (EPs) that emerge in a disorder-averaged Green's function description of two-dimensional (2D) Dirac semimetals with spin- or orbital-dependent potential scattering. Remarkably, we find that EPs may promote the nearly vanishing conductance of a finite sample at the Dirac point to a sizable value that increases with disorder strength. This striking behavior exhibits a strong directional anisotropy that is closely related to the Fermi arcs connecting the EPs. We corroborate our results by numerically exact simulations, thus revealing the fingerprints of characteristic non-Hermitian spectral features also on the localization properties of the considered systems. Finally, several candidates for the experimental verification of our theoretical analysis are discussed, including disordered electronic square-net materials and cold atoms in spin-dependent optical lattices.


I. INTRODUCTION
The study of disorder in quantum materials has led to the discovery of numerous intriguing phenomena, prominently including metal-insulator transitions and disorder induced topological states [1][2][3][4][5][6]. Quite generally, scattering off random impurity potentials tends to inhibit bulk transport by localizing the current-carrying Bloch states of a clean system [1,2]. In this context, the recently proposed framework of many-body localization has aimed at generalizing the basic picture of Anderson localization beyond the independent electron approximation [6][7][8].
Despite the Hermitian character of the microscopic Hamiltonian, non-Hermitian (NH) physics naturally emerges in the description of disordered systems due to the presence of a self-energy Σ = Σ H +iΣ A that accounts for impurity scattering in the disorder averaged Green's function. As a simplest scenario, the anti-Hermitian part Σ A of the self-energy just gives a scattering-induced finite lifetime to the stationary Bloch states of the free system, while the Hermitian part Σ H amounts to a correction of the free band structure described by the Bloch Hamiltonian H 0 (k) with the lattice momentum k. However, as soon as Σ A does not commute with H 0 (k), remarkable phenomena unique to NH matrices may occur, such as the formation of topologically stable exceptional points (EPs) [48][49][50] representing the NH counterpart of topological semimetals known from the Hermitian realm, and thus an example of NH topological phases  (see Ref. [37] for an overview).
In this work, we study the occurrence of disorderinduced exceptional phases in two-dimensional Dirac semimetals with spin-or orbital-dependent impurity scattering (see Fig. 1). To this end, we compute the effective NH Hamiltonian H e (k) = H 0 (k) + Σ(k, ω = 0), from exact numerical data on the disorder averaged Green's function [65] at the Fermi Energy (ω = 0) (cf. Section III below). While previous work on disorder-induced NH phases has mostly concentrated on spectral properties, our main focus is on the mesoscopic quantum transport properties of exceptional NH phases in a simple two-terminal setting [66]. Interestingly, we find that the fast exponential decay of the zero-energy conductance with system size, which is caused by the vanishing densitiy of states (DOS) at the nodal points in the clean regime, can be prolonged extensively in the presence of disorder-induced EPs. As a result, the nearly vanishing conductance of a clean mesoscopic sample with the Fermi energy at the Dirac point may be promoted to a sizable value that increases with disorder strength (see Fig. 1(c), and Fig. 6 below, respectively). Both the directional anisotropy and the value of the transmission are found to be closely related to the orientation and length of the Fermi arcs connecting the EPs as well as the associated quasiparticles with arbitrarily prolonged lifetime that emerge in a certain parameter range. Comparing our results to a conventional disordered phase without EPs, where transport is overall strongly damped and the exponential decay of the zeroenergy conductance occurs even faster than in the clean regime, our findings presented in Section IV provide clear fingerprints of disorder-induced EPs on basic transport properties. We find that the aforementioned quasiparticles with largely extended lifetime in the exceptional phase are tied to exact Bloch wave eigenstates that emerge in the limit of orbital-selective disorder and persist for any perturbation amplitude. In section V B, this scenario is generalized to a theorem that guarantees the presence of Bloch eigenstates, whenever only a sufficiently constrained subset of the orbitals are affected by disorder in an arbitrary tight binding model.
Additionally, we investigate the localization properties of the different phases in section V, and find hallmarks of anomalous localization (i.e. slower than exponential decay of the eigenfunctions) in the exceptional phase for system sizes that are amenable to numerical studies.
Finally, in Section VI, we outline several candidate platforms for the experimental verification of our theoretical analysis such as magnetically disordered 2D electron systems, square net materials, which refer to a class of three-dimensional (3D) multiatomic crystals with squareshaped 2D sublattices [72,73] hosting a Dirac-like dispersion, and cold atoms in spin-dependent optical lattices.

II. MODEL AND METHODS
We consider a conceptually simple square lattice model of a Dirac semimetal in two spatial dimensions (2D), specified by the Hamiltonian H 0 in second-quantized tight-binding form as with j = (j x , j y ) as well as δ x = (1, 0), δ y = (0, 1), and σ µ , µ = x, y, z denoting the standard Pauli matrices. Length is measured in units of the lattice constant, energy in terms of the spin-dependent nearest neighbor hopping chosen as t x = t z = 0.5, and the spinors of fieldoperators ψ † j = (ψ † j,↑ , ψ † j,↓ ) create a fermion in unit cell j with spin ↑ (↓). In reciprocal space, the model is characterized by the Bloch Hamiltonian H 0 (k) = d R (k) · σ with d R (k) = (cos(k y ), 0, cos(k x )), the spectrum of which exhibits four nodal points in the first Brillouin zone (see Fig. 2).
Random disorder with a spin-dependent structure is introduced to the system by means of the operator V , given in real space representation as with s 0 , γ, φ ∈ R. The uncorrelated random amplitudes {a j } are drawn from the box distribution on the real interval [−α, α], so the overall disorder strength scales with α. Depending on the choice of disorder parameters, an exceptional or an ordinary phase may emerge.
In Section VI, we outline several possible platforms for the experimental implementation of our model. We again stress that the total Hamiltonian H = H 0 + V is Hermitian, and that non-Hermitian physics only enters via the self-energy in the Green's function description discussed in the following.
To describe the electronic excitations close to the Fermi level in disordered systems, we employ an effective NH Hamiltonian H e (k) = H 0 (k) + Σ(k, ω = 0) (cf. Eq. 1), which we obtain from numerically exact data on the disorder-averaged Green's function. Starting from the is the Heisenberg picture spinor of creation operators in reciprocal space, a Fourier transform with respect to (t − t ) yields the retarded Green's function in frequency space. Averaging over the random disorder terms restores translational invariance and renders the the resulting disorderaveraged Green's function G R,av k block diagonal in momentum space (see Appendix A for details). The blocks G R,av where η > 0 is an infinitesimal regularization, determine the self-energy correction Σ(k, ω) that enters the effective Hamiltonian H e (k). In the following, we calculate G R,av k in a numerically exact fashion for a system size of 100 × 100 sites.

III. SPECTRAL PROPERTIES
In this section, we briefly discuss the phenomenon of EPs and show how a conventional disordered phase without EPs or an exceptional phase featuring EPs emerges in our model system specified by Eqs. (2, 3). (2)) exhibiting 4 Dirac cones with the parameters set to tx = 0.5, tz = 0.5.

A. Exceptional points
For parameter-dependent NH matrices, EPs refer to points in parameter space at which not only the eigenvalues but also the corresponding eigenvectors coalesce, thereby rendering the matrix non-diagonalizable [61][62][63]. Since all Hermitian matrices possess a complete set of eigenvectors, this is a genuinely NH phenomenon. For an overview on the role of EPs in contemporary physics, see, e.g., Refs. [37,64].
As an example for EPs occuring in a NH Bloch band setting such as the disorder averaged description of our disordered model system, here we focus on the case of a two-banded model with a generic NH Hamiltonian of the form where d = d R + id I with d R , d I ∈ R 3 and d 0 ∈ C. The eigenenergies of H e are E ± = d 0 ± d 2 R − d 2 I + 2id R · d I and at degeneracy points where the real and the imaginary part under the square root vanish simultaneously, i.e. d 2 R − d 2 I = 0 and d R · d I = 0, the eigenvectors coalesce, leading to an EP. In conclusion, EPs require the tuning of two separate conditions, which is also the case for more than two bands, and thus are a generic and stable phenomenon in systems with at least two spatial dimensions.
From a mathematical point of view, EPs are the end points of a branch cut in the complex plane, which manifests itself as a line of purely real or purely imaginary energy gap that connects the EPs and corresponds to the conditions d R · d I = 0 and d 2 R − d 2 I > 0 or d 2 R − d 2 I < 0 in the two-banded case. In the context of NH physics, these lines are referred to as (NH) Fermi arcs and they play a central role in the interpretation of the anomalous transport properties of the model studied in the following.
B. Complex spectrum of the model For the system at hand, the simplest possible perturbation, consisting only of an on-site disorder potential, is reflected by the parameter choice s 0 = 1, γ = 0, and leads to a trivial self energy Σ, such that the real part of the spectrum remains the same as for H 0 (cf. Fig. 2) and the imaginary part of the spectrum is flat and twofold degenerate with a value of about −0.4i.
Upon imbalancing the amplitude of the disorder potential between the orbitals, the system enters an exceptional phase. We start by considering the extremal case of the disorder only affecting the A-sublattice through the parameter choice s 0 = 1, γ = 1 and φ = 0, noting that deviations from this limiting case will be discussed further below. The Dirac cones are split into EPs connected by a Fermi arc of purely imaginary energy gap (marked by red lines in Fig. 3a), which increases in length with the disorder strength α and is perpendicular to the k x -direction. By tuning the parameter φ in Eq. (3), the angle between the Fermi arcs and the k x -axis may be adjusted, where the limiting cases φ = 0 (disorder on the A-sublattice) and φ = π (disorder on the B-sublattice) afford a simple physical interpretation. Other values of φ amount to a rotation around the σ y axis, thus resulting in orbital-selectivity in a rotated orbital basis. More generally, the parameter φ allows us to investigate the influence on the transmission of the angle between the Fermi arcs and the transport direction (cf. 5) in our setup. In an actual experiment, the angle dependence could be readily probed by measuring the transmission in different directions through the sample.
Another notable feature of the exceptional phase with orbital selective disorder is the presence of states with a vanishing imaginary part at all energies (marked by green lines in Fig. 3a). As we demonstrate in section V B, these states are Bloch modes of the clean system that remain unaffected by the disorder potential and are responsible for the enhanced and anisotropic transport capabilities of the exceptional phase.
Since the self-energy Σ shows no dependence on k (see Appendix A), the splitting looks similar around all nodal points. It is therefore sufficient to closer inspect the two Fermi arcs in the (+, +) and (+, −) quadrant and a cut along these two Fermi arcs is shown in Fig. 3b. The nodal points in the spectrum of H 0 from Fig. 2 have been inflated into nodal lines in the real part of the spectrum of H e . By stark contrast to the conventional phase, the DOS at zero energy does not vanish in the quasiparticle picture. To gauge the sensitivity of the exeptional phase towards deviations from the restriction of the disorder potential to one orbital, we also show the result for γ = 0.9 and γ = 0.8 in Fig. 3b, which corresponds to distributing 5% and 10% of the disorder amplitude on the second orbital. With decreasing value of γ, the Fermi arcs shrink and the top of the "Fermi bubbles" starts to move away from zero, i.e. the long-lived states acquire a long but finite life-time. . Lines of vanishing imaginary part (marked in green) entail propagating quasiparticles with infinite lifetime. The group velocity along the green lines is parallel to the kx-axis, which leads to a directional dependence of transport signatures. Disorder parameters are α = 1.5, s0 = 1.0, γ = 1, φ = 0 (b) Cut along the line where kx = π 2 , showing the profile of the two Fermi arcs in the (+, +) and (+, −) quadrants of Fig. 3a. Note that the imaginary "Fermi bubbles" touch the real axis. Disorder parameters are α = 1.5, s0 = 1.0, φ = 0. Different colors correspond to γ = 1 (blue), γ = 0.9 (orange), and γ = 0.8 (red), respectively, where for γ < 1 the disorder potential is no longer exactly restricted to one orbital.
For all sets of disorder parameters, we find only a very weak if any ω dependence of H e (k, ω), which allows us to draw conclusions on the quasiparticle-dispersion through the complete energy range from H e (k, ω = 0) or Fig. 3a, respectively.

IV. TRANSPORT SIGNATURES
Here, we investigate the influence of the quasiparticle dispersion inside the bulk on the transport properties of a finite-sized sample for the different NH phases. To this end, we generate a random rectangular instance of the system with N x = 100 sites in x-direction and N y = 200 sites in y-direction. After attaching two spin-independent metallic leads to both ends of the sample, we calculate the two-terminal linear-response conductance in x-direction. The dispersion in the leads is given by the Bloch Hamiltonian H L = (cos(k x ) + cos(k y ))σ 0 .

A. Energy-dependent conductance
First, we present the energy-dependent conductance through the clean system in Fig. 4a. Resonance effects lead to small fluctuations of the transmission anplitude over the energy spectrum and the vanishing DOS at zero energy completely suppresses transport there. Fig. 4b shows the energy-dependent conductance for different disorder configurations. The green line represents the conventional disordered phase, where almost no transport occurs at any energy, which is a consequence of the large imaginary part and therefore damping of modes in the corresponding quasiparticle dispersion. The orange line belongs to the exceptional disordered phase with φ = 0, where the Fermi arcs are perpendicular to the direction of transport along k x . There, an enhanced conductance at all energies when compared to the conventional phase is visible, carried by the states with an extended lifetime in the vicinity of the green lines in Fig. 3a. Strinkingly, the transmission at zero energy surpasses that of the clean system by far. We interpret that as a signature of the increased DOS in the quasiparticle picture caused by the inflation of the nodal point into a Fermi arc (compare to Fig 3b). Finally, the blue line indictates the exceptional disordered phase with φ = π/2, i.e. with the Fermi arcs parallel to direction of transport. Here, transport is suppressed heavily in comparison to the case with φ = 0.
B. Effect of φ on the conductance As Fig. 4b shows, the angle of the Fermi arc relative to the direction of transport affects transport at all energies in the exceptional phase. More specifically, the transmission decreases with increasing φ, which can be intuitively understood from the quasiparticle dispersion. In the imaginary part of the spectrum for φ = 0 in Fig. 3a, ridges with vanishing imaginary part are visible and accentuated by green lines. Due to their prolonged or even infinite (at γ = 1) lifetime, it is reasonable to expect that these quasiparticles and the ones close to them mostly carry the transport. Additionally, all of them have a group velocity parallel to the k x -axis, i.e. the direction of transport. With increasing φ, the Fermi arc and also the ridges rotate in the k x -k y -plane, while the direction of transport remains along the k x -axis. The group velocity of the states associated with the ridges rotates as well and always encloses the angle φ with the direction of transport. An elementary trigonometric construction shows that the wave fronts will have to travel a longer distance of N x 1 + tan 2 (φ) through the sample (see Appendix B), which causes the decay of the transmission amplitude.
To further substantiate this influence of the parameter φ on the transmission, we present the conductance in the exceptional phase at zero energy in dependence on φ in the main plot of Fig. 5 for a fixed system length of N x = 100 sites, along with the dependence on system length in the insets (a) -(c). In order to mitigate the influence of the coupling to the leads and varying degrees of reflection with φ, we work close to the wideband limit by flattening the dispersion in the leads to H L = 0.1(cos(k x ) + cos(k y ))σ 0 for the calculations in We start by discussing the insets (a), (b) and (c) of Fig. 5, which show the decay of the transmission in dependence of system length for φ = 0 along with a fit of the law for decreasing values of p. Specifically, we have p = 1 in (a), p = 1 2 in (b), and p = 1 5 in (c). Apart from small resonance-induced fluctuations at short system lengths, the fit in inset (c) with p = 1 5 is closest to the actual data and yields the parameters A = 4.63, τ = 1.06.
This fit to the length-dependent transmission is motivated by previous work on disordered one-dimensional systems suggesting that such a slower-than-exponential decay of the transmission amplitude is tied to anomalous localization properties [67][68][69]. There, the value of p that determines the length dependent transmission decay corresponds to the value obtained from the localization of the wave-function as Ψ(x) ∝ exp [−λ|x| p ], where p = 1 for standard Anderson localization. Since our fits agree best with the if we choose p = 1 5 or smaller, it seems plausible that the exceptional phase features very weakly localized microscopic eigentstates, which is also indicated by a more careful analysis of the localization properties in section V.
Moving to the main plot of Fig. 5, a decay of the conductance with φ is clearly visible. As mentioned before, the distance that the modes travel through the sample increases with φ as N x 1 + tan 2 (φ). Taking into account the findings from the insets, the φ-dependent conductance should roughly follow a law of the form T (φ) = e 2 h exp A − τ N x 1 + tan 2 (φ) p with p = 1 5 . A fit reveals a good agreement and basically recovers the parameters from the fit to the length-dependent decay with A = 4.67, τ = 1.08.

C. Effect of the disorder amplitude on the conductance
In Fig. 6, we analyze the dependence of the conductance in the exceptional phase on the disorder strength α. The parameters s 0 = 1 and φ = 0 are fixed and we explore the behaviour for the "ideal" case γ = 1 with infinitely lived quasiparticles as well as γ = 0.9, and γ = 0.8, where the disorder potential is no longer restricted to precisely one orbital such that the Fermi bubbles no longer touch zero (cf. Fig. 3b). Quite remarkably, both for γ = 1 and γ = 0.9, the conductance increases with disorder strength, together with the length of the Fermi arc which is plotted alongside. For γ = 0.8, the conductance exhibits a similar tendency at small disorder amplitudes, but drops again for strong disorder, where the increasing overall damping exceeds the inflation of the Fermi arcs.
The conductance enhancing effect of disorder in the exceptional phase is a consequence of the growing quasiparticle DOS and size of the imaginary Fermi bubbles (cf. Fig. 3b), which increases the amount of states with pro- longed life-time (suppressed damping). In this sense, the quite counter-intuitive observation of a conductance that increases with disorder strength affords a simple explanation within the framework of our effective NH Hamiltonian analysis. Our transport simulations are performed using the Kwant library [75].

V. MICROSCOPIC LOCALIZATION PROPERTIES
Here, we turn to the microscopic properties of the system. We calculate the inverse participation ratio (IPR) and the actual DOS obtained from a full diagonalization of the microscopic Hamiltonian. The IPR is a measure for the localization of states, which is given by dr|Ψ(r)| 4 . Fig. 4c shows the microscopic DOS and the averaged IPR at each energy for the disordered system in the conventional phase with the disorder parameters set to s 0 = 1, γ = 0, α = 1.5. The localization strength varies with energy and is weakest at zero energy with an IPR of about 0.008. It continuously increases when moving away from this point. The DOS at zero energy does not vanish (in contrast to the clean system) as a consequence of the σ 0 -term in the disorder (cf. Eq. (3)) that breaks the chiral symmetry of H 0 .

A. Microscopic DOS and average IPR
For the exceptional phase in Fig. 4d with the disorder parameters set to s 0 = 1, γ = 1, φ = 0, α = 1.5, the states are uniformly localized for energies between -1 and 1, with an IPR of about 0.007. Outside of that range, the localization grows stronger rapidly. The interval of relatively weak localization in the exceptional phase roughly coincides with the energy range in which transmission occurs (c.f. Fig. 4b) and also with the range that is covered by the real part of the spectrum (c.f. Fig. 3a). Although the microscopic DOS covers a larger energy interval than in the conventional phase, no qualitative difference is visible.

B. Surviving Bloch modes in the exceptional phase
While the average IPR around zero energy is only marginally smaller in the exceptional phase (cf.  Fig 4c, 4d), we find a clear distinction in the distribution of the individual IPRs of the microscopic eigenstates between the two phases. Fig. 7 compares the individual IPR values in both phases as a scatter plot over the energies of their respective eigenstates. The conventional phase in 7a features eigenstates with varying degree of localization such that the minimum value of the IPR at zero energy is about 0.003 or higher. By contrast, the exceptional phase depicted in 7b possesses a set of completely delocalized states that come in clusters of four throughout the energy range from -1 to 1. These delocalized states are Bloch waves of the clean system that live entirely on the sublattice unaffected by the random disorder and thus survive under any perturbation amplitude. This phenomenon is not a fine-tuned peculiarity of our model but rather a manifestation of a more general theorem: Consider an arbitrary tight-binding model H 0 with translational invariance and a number of n orbitals. For a random disorder term V that only affects a subset of n α < n orbitals, there can be Bloch modes of H 0 that live on the remaining n − n α orbitals and thus are still eigenstates of the disordered system H = H 0 + V regardless of the concrete nature of V . Such modes require the tuning of n α complex conditions in momentum space, which can be made real or reduced in number by the presence of certain symmetries (see Appendix C for a proof ).
These insights complement the analysis of Ref. [70], where the existence of zero-energy (localized) eigenstates that are restriceted to one sublattice of a disordered bipartite lattice has been derived. The surviving Bloch modes make up the quasiparticles with vanishing imaginary energy that we observe in the spectrum of the effective Hamiltonian (cf. Fig. 3a). To illustrate this, we represent the inifinitely-lived eigenstate of H e with momentum k x = π/2, k y = π/2 (in the center of the Fermi arc) in the basis of the microscopic eigenstates of the system and visualize the result in Fig. 7b by the color and size of the plot points. The quasiparticle is comprised entirely of the four delocalized Bloch eigenstates from the cluster at zero energy. The other three eigenstates of H e with zero (real) energy and infinite lifetime belonging to the other three Fermi arcs can also be represented within the same cluster. A similar outcome is observed for the non-decaying quasiparticles at other energies.
On the other hand, the quasiparticles with a large imaginary energy part are composed of a wide range of microscopic states with high IPR-values. Exemplarily, we show the overlap of one of the quasiparticles in the conventional phase (where all energies have an imaginary part of about −0.4i) with the microscopic eigenstates of the system in Fig. 7a.
The symmetry of our model system reduces the one complex condition, that arises from a disorder term restricted to one orbital, to a real one (see Appendix C). In conclusion, the surviving Bloch states have codimension one and should form one-dimensional submanifolds in the momentum space of our two-dimensional system, which are precisely the green lines in Fig. 3a.
Even though there is no direct connection between the surviving Bloch modes and the occurrence of EPs in the spectrum of the effective Hamiltonian, the two effects play together to create the observed transport phenomenology. While the Bloch states themselves only provide a single transport channel for a given energy regard- less of the system size, they enforce the lines of vanishing imaginary energy in the spectrum of the effective Hamiltonian (cf. Fig. 3a) and thereby a population of quasiparticles with a very long lifetime in the immediate vicinity of these lines. The transmission amplitude carried by these regions of the spectrum does also scale with system size. Since the distance between the EPs grows with disorder strength, the associated Fermi bubbles broaden and bring the imaginary part of the spectrum surrounding the Bloch states closer to zero (cf . Fig 3b), thereby increasing the transmission amplitude (cf. Fig. 6).

C. Localization behaviour
Based on the methods used in an earlier study on anomalous localization [70], we measure the spatial decay of the probability amplitude relative to the maximum by defining the correlation function Here, we consider the probability density summed over the two internal degrees of freedom, i.e. |Ψ| 2 = |Ψ A | 2 + |Ψ B | 2 , let |Ψ max | 2 denote the maximum, |Ψ(r, ϕ)| 2 the probability density in relative polar coordinates to the location of the maximum, and ... the average over multiple eigenstates. For an average localization behaviour Ψ(r, ϕ) ∝ exp −λr p(ϕ) with 0 < p(ϕ) < 1), the correlation function should behave as g(r, ϕ) = λr p(ϕ) . In general, our data agrees with the above ansatz and suggests that there is no dependence of the value of p on ϕ for larger systems with similar dimension in x and y direction. To properly capture the asymptotic localization behaviour, we investigate systems with a size of N x = 5000, N y = 200 by using the Krylov solver of the sparse matrix library of Scipy [71] to obtain 100 eigenstates close to zero energy of a randomly generated instance of the system, where we exlude the surviving Bloch eigenstates in the exceptional phase. From those, we compute the correlation function as per Eq. (5) in the direction ϕ = 0, i.e. along the x direction. Fig. 8a shows g(r, ϕ = 0) for the conventional phase, which behaves roughly linear up to the point where Ψ(r, ϕ = 0) becomes smaller than machine precision and the correlation function saturates. A fit of the form g(r, ϕ = 0) = λr p yields λ = 0.5, p = 0.73. For other sample geometries and perturbation amplitudes α, given that the length of the sample in the direction of ϕ is big enough, we observe similar behaviour and obtain fit parameters p ranging between 0.7 and 1, which we deem compatible with regular Anderson localization.
For the exceptional phase in Fig. 8b, a much slower and non-linear decay of the correlation function is observed. At the maximum system length presented here, g(r, ϕ = 0) is far from saturating and the fit g(r, ϕ = 0) = λr p reveals the parameters λ = 6.82, p = 0.07. Again, we observe similar behaviour for other sample geometries and obtain values of p around 0.1 from the corresponding fits. In summary, our data provides numerical evidence that the exceptional phase of our model exhibits anomalous localization. This is also consistent with the subexponentially decaying transport amplitude (cf. Fig. 5).

A. Basic requirements
Besides the intriguing algebraic aspects of EPs by themselves, a core ingredient for the physics discussed in this work is the possibility of quasiparticles with infinite lifetime in a dissipative environment caused by disorder scattering. Such extraordinarily long-lived modes are associated with Fermi-arc regions in between the EPs, where the imaginary part of the complex spectrum of H e (at least approximately) touches the real axis (cf. Fig. 3a). For this to occur, we have identified three crucial ingredients.
(i) First, a chiral-symmetric base model H 0 , which has symmetry-protected nodal points. The symmetry implies that the Bloch Hamiltonian H 0 (k) contains only two of the three Pauli matrices σ x,y,z (or two linearily independent combinations of them) and reduces the complex condition from the theorem in Section V B to a real one and thus ensures the survival of Bloch modes at all energies.
(ii) Second, a random disorder with on-site terms of the form with |s 0 | = |s|. The vector s must be chosen such that s · σ respects the chiral symmetry of the base model H 0 . This ensures that the total model can be brought back to a form where the disorder potential only affects one orbital and the Bloch Hamiltonian possesses eigenstates that live on the opposite orbital by a unitary transformation.
(iii) Third, a symmetric distribution of the random amplitudes {a j } around zero is required. Any deviation from a symmetric distribution reduces the length of the Fermi arcs and shifts them away from zero energy.

B. Square-net materials
As a first candidate, we consider a setup where the two internal degrees of freedom are given by different atomic orbitals. We may choose s = (0, 0, 1) to represent disorder on only one of the two orbitals, since A model of this type may occur within the class of multiatomic crystals with square sub-nets which have been found to exhibit a Dirac dispersion [72,73]. A remaining challenge then is to restrict the disorder potential (at least to good approximation) to one of the two orbitals of the atomic species that forms the square net. In many cases though, the orbitals are p-like as in LaCuSb 2 , where the Sb atoms form the square net that yields Dirac physics [74]. Assuming a crystal growth direction that results in a square parallel to the surface, positive and negative ions adsorbed to the surface may affect one of the orbitals much stronger, since they are club-shaped in different spatial directions, thus creating a disorder term similar to Eq. (7). Our calculations show that the disorder does not have to be restricted to one orbital perfectly. In particular, if the disorder strength on one of the orbitals is about 5 % of the other one, the quasiparticles along the green lines in Fig. 3a obtain a small imaginary part, but the characteristic transport signatures still remain (cf. Fig. 3b  and 6). Only at about 10 % deviation from the ideal configuration, the characteristic signatures start to become unrecognizable.
According to the above ingredient (ii), the Bloch Hamiltonian H 0 (k) of the base model should posses a chiral symmetry that allows σ z and one other Pauli matrix. However, this should not be problematic, since the setup with atomic orbitals exhibits real hoppings. Together with the square-shaped lattice, a chiral symmetry that rules out σ y and permits σ x , σ z is to be expected.
Along the same lines of reasoning, honeycomb materials such as graphene or silicene as well as other Dirac materials with non-rectangular lattices do not represent promising candidates to create an exceptional phase akin to the one discussed in this work. The chiral symmetry of such lattices will generally permit σ x and σ y in the Bloch Hamiltonian H 0 (k), which means that the condition for the existence of a disorder-insensitive Bloch mode is complex. Thus, their codimension is two and they will only appear at certain isolated momenta. This is confirmed by the numerical studies we conducted on disordered Honeycomb models, where Bloch modes survived at isolated momenta under the presence of lattice-selective disorder but no EPs where observed. However, lattice-selective disorder in graphene has been linked to other interesting phenomena by a previous study [51,52].
C. Ultracold atoms in optical lattices with spin-selective disorder In the context of atomic many-body systems, a promising platform for experimentally implementing our model system is provided by ultracold atoms in (spindependent) optical lattices [53]. In the quest for probing many-body localization, the experimental realization and control of disorder potentials in such synthetic material systems in both 1D and 2D has been impressively demonstrated [54,55]. Moreover, quantum gas microscopy methods with spin-selective single site resolution have enabled the observation of anti-ferromagnetic phases with cold atoms in optical lattices [56][57][58]. Combining these recent developments, our proposed 2D Dirac models with spin-selective disorder potential are well within the state of the art toolbox of atomic quantum simulators. Remarkably, the enormous experimental control over such systems does not only allow for the observation of spectral properties: Despite the electrically neutral character of the atoms, also two terminal transport setups have become amenable in the laboratory [59,60].

D. Surface states of a topological insulator
Another experimentally well studied platform for 2D Dirac physics is based on the surface states of a timereversal symmetric topological insulator (TI) exposed to magnetic perturbations that affect the electron spin representing the internal degree of freedom in this setting. There, an exceptional NH phase has previously been proposed to be induced by the tunnel coupling to a ferromagnetic thermal reservoir [44]. These ideas may be adapted to our present context of NH physics induced by disorder scattering in a straight-forward way. Instead of a magnetic lead, the relevant perturbation is then given by magnetic impurity ions adsorbed to the surface such that the electro-static potential generates the s 0 σ 0 disorder and the magnetic moment the s·σ term (cf. Eq. (6)). However, it is fair to say that several issues remain with the disordered version of this setting. To satisfy ingredient (ii), fine-tuning the amplitude of s 0 and s is required. Furthermore, the above ingredient (iii) requires that the on-site potential s 0 σ 0 changes sign with the direction of the magnetic moment s · σ. This amounts to having pos-itive and negative magnetic ions with oppositely oriented magnetic moments.

E. Comparison to materials with tilted Dirac cones
Previous works have demonstrated that EPs can emerge from trivial disorder (i.e. affecting all degrees of freedom similarly) as well in two-dimensional Dirac materials with tilted cones [48][49][50]. While such systems are probably easier to fabricate than the ones listed so far, the resulting Fermi arcs are unfortunately buried under a finite imaginary part that increases with the size of the Fermi arcs, since the non-decaying quasiparticles can only be obtained through orbital-restricted disorder. This is likely to obscure transport signatures and renders this class of systems less interesting for our present purposes.

VII. CONCLUSION
In this work, we have demonstrated the occurence of disorder-induced EPs by means of exact numerical calculations and explored their impact on the transport properties of a finitely-sized sample. For this, we study a minimal model of a Dirac semimetal, which can enter a conventional phase (without EPs) or an exceptional phase through the addition of random disorder on both or predominantly on one of the sublattices. The exceptional phase exhibits quasiparticle excitations with prolonged lifetime that we tie to orbital restricted disorder through a rigorous theorem. Besides studying the spectral properties, we compute the two-terminal transmission, and find that a finite sample in the exceptional phase shows greatly enhanced transport properties at zero energy when compared to both the clean system and the conventional disordered phase. We demonstrate that this effect is carried by the aforementioned long-lived quasiparticles, and accompanied by an anomalously slow decay of the localized eigenfunctions. Finally, we outline possible platforms for the experimental realization of the proposed model systems, where the main challenge is represented by engineering orbital-selective disorder.
can self-insert the Dyson equation to obtain a perturbation series for the full retarded GF If the unperturbed system H 0 is a tight-binding model with translational invariance and n internal degrees of freedom on each site, the free GF is block-diagonal in the Bloch basis and reads G 0 We will consider the case of random disorder represented by similar impurities at each site, but with uncorrelated random amplitudes a j (with j ∈ {1, 2, ..., N x } × {1, 2, ..., N y }) distributed according to a distribution function f (a). The disorder term V may include on-site terms and hoppings between sites connected by a vector δ = (δ x , δ y ). The generic form in real space representation is then The Ψ † j are n-spinors of creators for on-site states and the c † k = 1 √ NxNy j e ij·k Ψ † j are n-spinors of creators for the Bloch-states. The series from Eq. (A1) leads to an expression for the blocks of the full GF To obtain the effective Hamiltonian, this expression is averaged over the impurity amplitudes {a j } by calculating G av k,k (iω) = G k,k (iω) f = da 1,1 f (a 1,1 ) da 1,2 f (a 1,2 )... da Nx,Ny f (a Nx,Ny )G k,k (iω).

Now we need to introduce a small error of the order 1
NxNy by letting the h sums run unrestricted, e.g.
The average doesn't act on the exponentials and the distribution of the amplitudes is uncorrelated, so we can pull the average of the amplitudes out of the sums. Performing the sums gives a delta function for all momenta connected to the same impurity and we arrive at For on-site impurities, i.e. impurities without hopping terms, Eq. (A2) shows that V k,k is a constant and thus the self-energy Σ(k, ω + iη) does not depend on k as well, since we integrate over all internal momenta.  Consider a generic tight binding Hamiltonian H 0 with n internal degrees of freedom and translational invariance perturbed by a random disorder term V that only affects a subset {α} of the orbitals in the sense of random on-site potentials as well as random hopping terms connecting {α} orbitals from different unit cells. We will denote the complement of {α} by {β}, the number of disorder affected orbitals by n α and the number of the remaining orbitals by n β = n − n α .
Assuming periodic boundary conditions or the thermodynamic limit, we may diagonalize H 0 by switching to the Bloch basis with the transformation Ψ j,α = 1 √ N k e ik(r j +rα) c k,α , Ψ j,β = 1 √ N k e ik(r j +r β ) c k,β