{'fraction_non_alphanumeric': 0.0442152466367713, 'fraction_numerical': 0.015485799701046337, 'mean_word_length': 4.5208780326786595, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bulk-boundary correspondence in non-Hermitian systems: stability analysis for generalized boundary conditions

The bulk-boundary correspondence (BBC), i.e. the direct relation between bulk topological invariants defined for infinite periodic systems and the occurrence of protected zero-energy surface states in finite samples, is a ubiquitous and widely observed phenomenon in topological matter. In non-Hermitian (NH) generalizations of topological systems, however, this fundamental correspondence has recently been found to be qualitatively altered, largely owing to the sensitivity of NH eigen-spectra to changing the boundary conditions. In this work, we report on two contributions towards comprehensively explaining this remarkable behavior unique to NH systems with theory. First, we analytically solve paradigmatic NH topological models for their zero-energy modes in the presence of generalized boundary conditions interpolating between open and periodic boundary conditions, thus explicitly following the breakdown of the conventional BBC. Second, addressing the aforementioned spectral fragility of NH matrices, we investigate as to what extent the modified NH BBC represents a robust and generically observable phenomenon.


I. INTRODUCTION
In a broad variety of physical situations ranging from classical settings to open quantum systems, non-Hermitian (NH) Hamiltonians have proven to be a powerful and conceptually simple tool for effectively describing dissipation. In the classical context, including optical setups [1][2][3][4][5][6][7][8][9][10][11], electric circuits [12][13][14][15][16][17], and mechanical meta-materials [18][19][20], the equations of motion are naturally determined by a NH matrix. This situation may in many cases be mapped to an effective (tightbinding) Hamiltonian [14,20] familiar from the quantummechanical modeling of electrons in crystalline solids within the independent particle approximation, but with additional NH terms. Conversely, as an effective description for dissipative quantum systems, which, at a fundamental level, are governed by Liouvillian dynamics, similar NH models can in several scenarios be directly derived [21,22].
the occurrence of protected edge-modes, the breakdown of this direct relation has been demonstrated for various NH extensions of topological insulators [49,50,52,54]. Motivated by this observation, several approaches to formulate a modified NH BBC have been put forward [34,44,[54][55][56][57][58][59][60][61][62]. Since the mentioned breakdown of the standard BBC is closely related to known instabilities of the eigenvalue spectrum of NH matrices [63], the physical robustness of the NH BBC has been questioned [34,54], and analyzing the more stable singular value spectrum has been proposed as and alternative diagnostic tool [57].
In this work, we address several remaining issues regarding the BBC in NH systems, focusing on the biorthogonal basis approach reported in Ref. [54]. First, arXiv:1912.07687v1 [cond-mat.mes-hall] 16 Dec 2019 we demonstrate the robustness of NH zero-energy boundary modes against physically relevant local perturbations, despite the aforementioned fragility of the eigenstate spectrum towards adding generic (random) matrices to the NH Hamiltonian. Specifically, the matrixelements of a NH tight-binding model are found to typically decay fast enough in real space to overrule the rapidly growing relevance with spatial distance of (non-local) perturbations. Second, considering onedimensional (1D) NH systems, we derive analytical expressions for the occurrence of exceptional point (EP) transitions and the formation of zero-energy edge modes as a function of a generalized boundary condition parameter Γ continuously interpolating between periodic (Γ = 1) and open (Γ = 0) boundaries (see Fig. 1 for an illustration). In this context, the discrepancy between the topological phase diagrams of systems with different boundary conditions is intuitively explained by the occurrence of topological phase transitions, in which the boundary condition parameter Γ plays the role of a control parameter. Our findings provide additional insights on both the analytical origin and experimental relevance of anomalous (from a Hermitian perspective) bulk boundary effects in topological NH systems.
The remainder of this article is organized as follows: We start by formally introducing generalized boundary conditions in topological NH one-and two-band models in Section II. We shortly discuss previous approaches to restore the BBC for NH systems, thereby focusing on the biorthogonal framework proposed in Ref. [54] in Section III. After deriving the analytical EP transitions for two paradigmatic NH models in Section IV, the stability of their eigenspectra and edge modes is discussed in Section V. Finally, a concluding discussion is presented in Section VI.

II. NH TIGHT-BINDING MODELS WITH GENERALIZED BOUNDARY CONDITIONS
We focus on one-dimensional tight-binding models with N unit cells and generally asymmetric hopping amplitudes between the sites that render the Hamiltonian NH. Compared to their Hermitian counterparts, some NH tight-binding models show qualitatively different eigenspectra depending on the imposed boundary conditions [49,50,54], i.e. periodic boundary conditions (PBC) or open boundary conditions (OBC). In fact, the whole eigenspectrum and with it all the eigenstates are affected by the boundary conditions: In the OBC case, not only zero-energy edge modes familiar from Hermitian topological systems but also a macroscopic number of bulk modes may be exponentially localized at one of the edges -a phenomenon coined the non-Hermitian skin-effect [55,64], whose connection to the failure of the conventional BBC has been widely discussed [17,36,51,54,59,65]. Indeed, the NH skin effect is found to occur alongside the discrepancy between the PBC and OBC eigenspectrum in these models [58]. To examine the extreme sensitivity towards the boundary conditions, one may introduce generalized boundary conditions by scaling the hopping between the last site n = N and the first site n = 1 with a parameter Γ ∈ [0, 1] such that Γ = 0 (Γ = 1) corresponds to OBC (PBC) [52,54] as depicted in Figure 1. It has been observed [52,54,60] (e.g. by using a complex flux to tune the boundary conditions [60]) that the transition between the periodic and the open system happens almost instantaneously. That is, the eigenspectrum exhibits a qualitative change of order one as soon as the boundary conditions is modified by the critical parameter Γ c that is exponentially small in system size. In other words, Γ c ∝ exp(−αN ) with α some positive real-valued number depending on the system's parameters is enough to trigger a bulk transition in the spectrum [54]. A minimal model with asymmetric hopping amplitudes is the Hatano-Nelson model [66]. The real space Hamiltonian reads as [34] where N is the number of sites, c j (c † j ) are the fermionic annihilation (creation) operators on site j, t L , t R ∈ C are the hopping amplitudes, and Γ ∈ [0, 1] determines the boundary conditions. Its Hamiltonian in reciprocal space with PBC (Γ = 1) is given by where the lattice momentum k is summed over the first Brioullin zone (BZ). While in the Hermitian context, the meaning of a bulk gap is lost for one-band models, the complex spectrum of the Hatano-Nelson model exhibits a point-gap [34], as long as there is no state at E = 0, i.e. |t R | = |t L |. This one-band model possesses pointsymmetric eigenspectrum in the complex plane if the total number of sites is even. Hence, the eigenenergies come in pairs {E, −E} such that the eigenenergy E = 0 is degenerate for an even number of sites. Moreover, H(k) = 0 determines the zero energy, i.e. point-gap closing points [34], for the periodic system in exact analogy to Hermitian systems. Another emblematic class of models is provided by noninteracting two-band models which are fully described by their Bloch Hamiltonian H(k) in reciprocal space. The NH version of H(k) reads as where σ= (σ x , σ y , σ z ) are the standard Pauli matrices, k is the lattice momentum and d Re , d Im ∈ R 3 . Its eigenvalues are given by Since the Hamiltonian is NH, we encounter points in parameter space that not only feature degenerate eigenvalues but also their corresponding eigenvectors coalesce. Such points are called exceptional points (EPs) [67] that render the Hamiltonian defective (non-diagonalizable).
In the case of generic NH Bloch-Hamiltonians (see equation (3) and (4)), both eigenvalues coincide if E + = E − = 0 which generally is an EP of order two apart from the trivial case d Re = d Im = 0 known as the diabolic point. Note, that the one-band equivalent H(k) for the Hatano-Nelson model has no EPs because a scalar cannot be defective. There, however, the NH real-space tight-binding Hamiltonian can still exhibit EPs.
The unit-cell of the considered two-band tight-binding models consists of two alternating sites A, B. Having imposed periodic boundary conditions (PBC), the position space Hamiltonian H is obtained upon Fourier transform with the fermionic creation (annihilation) operators where the momentum k is summed over the BZ. Thus, H still encounters an EP at E = 0 whenever H(k) does, which at the same time marks the bulk-gap closing point.
We illustrate our further analysis based on a NH version of the SSH-model [68] that has been widely discussed in Refs. [54,55,57,58,60,62], where are the components of d in equation (3). Since this Hamiltonian preserves the chiral symmetry σ z H(k)σ z = −H(k), the eigenenergies come in pairs {E, −E} such that the eigenenergy E = 0 is degenerate if the total number of sites is even. In Fig. 2 the absolute values of the eigenenergies are depicted for a) PBC and b) OBC to show the qualitative differences of both spectra for this model. The bulk-gap closing points as a function of the hopping amplitude t 1 are determined by the roots of (4) in the periodic system while it has been found [54,55] that give the bulk-gap closing points for the open system. The NH SSH-model has thus either two ( t 2 > γ/2) or four bulk gap closing points for OBC [55]. Alongside the whole eigenspectrum, the position of the bulkgap closing points in parameter space is altered as well when imposing open OBC. The resulting qualitatively differing eigenspectra of the periodic and the open system (see differences in a) and b) of Figure 2 as an example) cause the conventional bulk-boundary correspondence (BBC) to break down [49,54,55]. Usually, this celebrated correspondence connects the occurrence and number of protected surface states of a Hermitian Hamiltonian with open boundaries to a topological invariant calculated from the bulk Hamiltonian, i.e. the system with PBC. Then, the bulk-band touchings mark the borders between different topological phases [53]. Introducing again the generalized boundary conditions parameterized by Γ ∈ [0, 1], the real space NH SSH Hamiltonian explicitly reads as In Section IV, analytical solutions for the bulk  [34]. The blue shaded area shows where a half-integer winding number (complex extension of the Chern-number [49]) indicates a non-trivial topological phase. The lighter green area shows where the biorthogonal polarization P = 1 [54] and edge states depicted in red occur while the darker green shaded areas show where no edge forms even though the periodic spectrum suggest a topological non-trivial phase (compare to shading in a)).
transition-points as a function of Γ for both the Hatano-Nelson model (1) and the NH SSH model (8) will be presented.

III. BIORTHOGONAL BULK BOUNDARY CORRESPONDENCE APPROACH
Recently, several suggestions for a non-Hermitian version of the bulk-boundary correspondence have been proposed [34,[55][56][57][58][59][60][61][62], including an approach based on the biorthogonal polarization P [54] on which we elaborate below. Accounting for the aforementioned sensitivity of NH systems to boundary conditions, the crucial difference to conventional bulk topological invariants is that P is calculated from bulk states of a system with OBC, thus accurately predicting bulk-band touching points and surface states for the open system (light green shaded area in Figure 2b)). Similar results have been obtained in a complementary way by considering a generalized BZ for NH systems [55] that contains model-specific information on both the periodic and the open boundary system. In contrast, there have been some notable suggestions for non-Hermitian topological invariants that are deduced solely from the periodic spectrum [33,34,39,49]. Firstly, the non-Hermitian half-integer winding number [49] (blue shaded area in Figure 2a)). Secondly, one can construct a different type of winding number based on the encircling of complex eigenenergies around the origin E = 0 in the complex plane that is trivially zero in all Hermitian systems [34] (red shaded area in Figure 2a)). However, both in general do not correctly predict the occurrence and disappearance of edge states in the open system for some parameter regime (compare to Figure 2 b)). Biorthogonal Polarizaton. We now briefly review the construction of the biorthogonal polarization P [54] for later reference. If the Hamiltonian is non-Hermitian, one has to distinguish between right and left eigenstates |ψ R,i and |ψ L,i found from the right and left eigenvalue problem, respectively: Choosing the biorthogonal normalization . . } that spans the complete eigenspace unless the system is at an EP. One-dimensional Bloch-Hamiltonians, where the vector d in (3) is given by d x (k) ± id y (k) = f ± + g ± exp(±ika) with a the lattice vector between neighboring unit cells, always possess an exact eigenmode with energy E = d z if the total number of sites is odd (e.g. the last unit cell is broken) and OBC are imposed [70]. Assuming the chain starts and ends with an A-site, the corresponding exact left and right eigenstates read with N the number of unit cells and c † n,A(B) [c n,A(B) ] the creation [annihilation] operators of a particle on sublattice A(B) in unit cell n, and where N R(L) is a normalization factor according to the biorthogonalization condition. These states are exponentially localized to one of the edges of the chain. In Ref. [54], it has been found that a topological phase transition in the open system occurs if the biorthogonal analogue of the mode localization changes from one side to the other. That is, the biorthogonal expectation value of the projector Π j = c † j,A |0 0|c j,A + c † j,B |0 0|c j,B onto the jth unit cell is exponentially localized to the left (right) edge for |r * L r R | > 1 (|r * L r R | < 1). Furthermore, the bulk gap of the open system closes if Using this, one constructs the biorthogonal polarization which is quantized in the thermodynamic limit and moreover jumps precisely at the bulk band touching points between 0 and 1. P thus plays the role of a bulk invariant characterizing systems with OBC by predicting edge states in the following sense: odd # of sites even # of sites P = 0 left localized edge state no edge states P = 1 right localized edge state two edges states

IV. ANALYTICAL SOLUTION WITH GENERALIZED BC
Going beyond previous literature, we now analytically derive the bulk transition points at the "critical" generalized boundary conditions Γ c = exp(−αN ) [54] for the NH models introduced in Section II, thereby specifying the dependence of α on model-parameters. To this end, we solve the right and left eigenvalue problems of the position space Hamiltonian H for the specific eigenvalue E = 0 with generalized boundary conditions. If PBC are imposed (Γ = 1), this ansatz gives the eigenstates associated to the gap-closing points of the Bloch (or Bloch-like) Hamiltonian. By following these eigenstates while Γ is continuously interpolated between 1 and 0, the exact parameter dependence of Γ c is revealed. We furthermore find an EP at zero energy that moves through parameter space as a function of the boundary condition parameter Γ for an even number of sites. In particular, this EP coincides with the gap-closing point whenever there is no edge mode in the OBC spectrum.
a. Nelson-Hatano model. Solving the right eigenvalue problem H HN |ψ R 0 = E|ψ R 0 (see eq. (7)) with the ansatz |ψ R 0 = N n=1 ψ R n c † n |0 and assuming that the eigen-energy is given by E = 0, we arrive at the following set of equations 0 = t L ψ n−1 + t R ψ n+1 , for n = 2, 3, . . . , N − 1, (15) with the boundary conditions Thus, the Hatano-Nelson model also possesses one exact zero-energy eigenstate throughout the whole parameter range with vanishing amplitudes on all even numbered sites if the total number of sites N is odd and OBC (Γ = 0) are imposed. If the number of sites N is even, we solve the system of equations for Γ and we find The corresponding right eigenvectors read where equation (17) [ (18)] ensures the existence of (19) [ (20)] and thus is a consistency relation. Since (17) and (18) cannot be fulfilled at the same time (apart from the special case |t L | = |t R |), the amplitude ψ n on every second site has to vanish and we only find one linearly independent eigenvector. However, the eigenvalue E = 0 is at least two-fold degenerate as we have pointed out in Section II and we thus encounter an EP. In order to arrive at the left eigenvectors t L ↔ t R have to be interchanged in (19) and (20). In Ref. [34] it was argued that the Nelson-Hatano model features a topologically non-trivial phase, in which an edge state forms for semi-infinite BC. We however will motivate in the following, that the Nelson-Hatano model does not have edge states for OBC which is confirmed by the numerical calculation of the eigenspectrum with exact diagonalization (see Fig. 3a)). Since the Nelson-Hatano model still has one exact zero-energy eigenmode |ψ R(L) for an odd number of sites, which explicitly reads as the biorthogonal polarization P [54] (cf. (14)) can still be constructed with the simplified projection operator Π n = c † n |0 0|c n . One already sees that the condition (13) is always met. Hence, the biorthogonal projection of these exact eigenstates |Ψ L , |Ψ R is never localized and P = 0. The only exception is if at least one of the hoppings t L , t R equal zero. At this point, one of the exact eigenstates is the zero vector and the other becomes ill defined. Hence, P is not defined. In fact, here we encounter a higher order (real-space) EP with eigenenergy E = 0 (see Figure 3a)). Note that the effect of higher order EPs on the failure of the conventional BBC has been discussed [52,71]. Since we find no edge states in the Hatano-Nelson model for an even number of sites and generalized boundary conditions, the exact gap closing points for a given t L , t R are found for Γ a(b) c (t L , t R ) (see Figure 3b) for an example). Figure 1b) thus visualizes equations (17) and (18) describing how the gap closing points move through parameter space when tuning Γ from PBC to OBC.
b. Non-Hermitian SSH-Model. Solving the right eigenvalue problem H SSH |ψ R 0 = E|ψ R 0 with the ansatz |ψ R 0 = N n=1 ψ R n,A c † n,A + ψ R n,B c † n,B |0 and assuming that the eigenenergy is given by E=0, we arrive at the following set of equations with n = 1, . . . , N − 1 in 23 and n = 2, . . . , N in 24.
The system of equations is decoupled with respect to the amplitudes on A(B)-sites, such that solving for Γ gives two solutions The left eigenvalue problem H † |ψ L interchanges t 1 − γ/2 ↔ t 1 + γ/2 in equation (23), (24) and (25). The existence of the eigenvalue E = 0 is guaranteed by equation (25) (23) and (24). The remaining non-zero amplitudes are defined apart from a normalizing factor, leading to only one linearly independent right eigenvector |Ψ R 0 and we encounter an EP. The order of this EP, i.e. the difference between algebraic and geometric multiplicity plus one, is at least two and E = 0 is at least two-fold degenerate. With exact diagonalization techniques, in can be observed that the EP at E = 0 in this model is in most cases indeed of order two. One exception can be found at the fine-tuned point t 2 = γ/2 and t 1 = 0. Several remarks are in order. First, we note that the EP found here is not necessarily the only one but others might occur at different eigenenergies E = 0. One known example is the (N − 1)-fold degenerate bulk-EP at E = ±1 that appears for OBC (Γ = 0) when |t 1 | = γ/2 [52,71] that is accompanied by an EP of order 2 with energy E = 0, which can be seen from (25). Second, we stress that the occurrence of the EP depends on the exact number of sites: If it is odd for instance, the eigenenergy E = 0 is no longer degenerate. Unlike the Hatano-Nelson model, the NH SSH model features edge states in the open system [54,55]. With our method alone however, we cannot decide whether the exact zero energy eigenmodes occuring at Γ c are bulk gap-closings (scenario (i)), or whether they are isolated edge states forming in the gap (scenario (ii)). To distinguish both scenarios we additionally compute the biorthogonal polarization P [54] (see also Section III). Scenario (i) (P = 0): Even though the PBC-Hamiltonian suggests a topologically non-trivial phase, there are no edge states found in the OBC eigenspectrum (dark green shaded areas in Figure 2b)). The EP at E = 0 is then a bulk-gap closing point, since it is supported by other modes. This means, that the gap-closing point shifts through parameter space as a function of the boundary conditions Γ c , which is the case for the Hatano-Nelson model. This happens whenever PBC and OBC differ in their topological phase. Scenario (ii) (P = 1): the model is expected to be in a topologically non-trivial phase judged from the periodic spectrum and also features edge states in the open system. Then, indeed, the EP at E = 0 that we find somewhere between OBC and PBC corresponds to an edge mode. For true OBC, the eigenvalue of this edge mode is not exactly zero but exponentially small. It is also no longer exactly degenerate and the eigenvectors are distinct but they only differ by exponentially small parameters.
To further stress the similarity to the construction of the biorthogonal polarization P , we rewrite (25) and arrive at where the argument of the logarithm must be ≤ 1, such that the logarithm is negative, and ξ A/B equals the localization length of the exact boundary modes (10,11) used for the biorthogonal polarization P . Note that this form stresses the great sensitivity of the eigenspectrum towards the boundary conditions (compare also to the corresponding Eq. (17) for the Hatano-Nelson model). Because of this similarity, a reformulation of (13) is Γ A(a) c · Γ B(b) c < 1 which gives the parameter regime for which edge states occur for OBC for both models. Analogously, Γ A(a) c ·Γ B(b) c = 1 matches with the condition that the biorthogonal Polarization needs to jump between 0 and 1. The case in which the left and right eigenvectors |ψ R , |ψ L used for the construction of P equal each other is reflected by Γ A c = Γ B c . One example is given by tuning γ = 0 in the NH SSH model (that is either the Hermitian limit if t 1 , t 2 ∈ R or if t 1 , t 2 ∈ C a version of the model discussed in Ref. [72]), and one generally can construct two linearly independent eigenvectors except from the cases Γ = 1 and Γ = 0. Then, the eigenspectra for OBC and PBC coincide apart from edge modes [54,72]. Additionally, by inserting parameter sets {t 1 , t 2 , γ} satisfying the conditions (7), i. e. parameter sets for which P jumps between 0 and 1, into Γ A(B) c , we define Γ P . Then, the bulk-gap closing point stays fixed for Γ ≤ Γ P . Furthermore, Γ P separates the scenarios (i) and (ii) distinguished above.
Quite remarkably, for both the Hatano-Nelson model and the NH SSH model, we find that the eigenenergy spectrum stops winding around the origin E = 0 in the complex plane precisely at the critical boundary conditions Γ c (see Fig. 4). That is, if the eigenenergy winds around the origin for PBC, it does for any Γ > Γ c . Since we have searched for the boundary conditions that feature an exact zero energy eigenmode, it is clear that a topological phase transition in the sense of Ref. [34] occurs if the base energy (here E = 0) is crossed as a function of the boundary condition parameter Γ. Even though we derived the zero energy eigenmodes for an even number of unit cells, this statement is general and valid for an odd number of cells as well. This finding provides an intuitive explanation for the discrepancy between the topological phase diagram for systems with PBC vs. systems with OBC, namely due to a topological phase transition in which the boundary condition parameter Γ plays the role of a control parameter. Thus, when extended to systems with generalized boundary conditions, the spectral winding number displayed in Fig. 4 captures the breakdown of the conventional BBC.

V. STABILITY OF EDGE STATES
We have seen that the eigenspectrum changes exponentially fast in system size when tuning its boundary conditions. In particular, we have seen for two examples that EPs at E = 0 start to move in parameter space as soon as the boundary conditions are modified. Since they coin- cide with the bulk-gap closing point in a wide parameter regime, their position in parameter space as a function of boundary conditions is a is an indicator for the transition between the qualitative PBC and OBC spectra. Hence, this transition goes with ∝ exp(−N ) [52,54,60] where N is the number of unit cells. This exponential sensitivity raises the natural question whether these apparently fragile spectral properties of systems with OBC are practically observable in realistic systems. On the other hand, an (unwanted) coupling amplitude between the first and the last site in a generic tight-binding model with linear geometry is also strongly suppressed in system size, and we encounter a problem of competing scales that we will take a closer look at in the following. The dramatic changes in the eigenspectrum of non-Hermitian Toeplitz-matrices (i.e. matrices with constant entries on all diagonals) and operators towards small perturbations have been already discussed in a mathematical context [34,63]. In particular, it was proposed to rather study so called -pseudoeigenvalues, that can become proper eigenvalues upon adding a suitable perturbation of norm . This argument was used in Ref. [34] to support the idea of investigating pseudo quasi edge states that stem from imposing half-infinite boundary conditions. As a complementary approach accounting for this spectral instability, Ref. [57] proposed to use the more stable singular value spectrum as and alternative to the directly observable eigenspectrum. Here, instead we argue that the considered perturbations should be physically motivated, and are by no means expected to be arbitrarily non-local random matrices as considered in Ref. [63]. Instead, realistic unwanted perturbations in a one-dimensional experimental setting with linear chain geometry are for instance given by couplings between sites with some larger distance that are omitted in a generic tight-binding model. Couplings that would effectively change the boundary conditions are then proportional to the overlap of the Wannier functions centered on the first and the last site. To quantify the generic scaling of this overlap, we may approximate the potential wells in the tight-binding formalism with harmonic potentials, where it is easy to see that the hopping between two sites i and j scales as ∝ exp(−|i − j| 2 ) (see the appendix VII). Thus, the matrix element coupling the last with the first site is of order exp(−N 2 ) which for large enough systems is drastically smaller than the scaling function Γ c ∝ exp(−αN ) that describes the transition between the OBC and the PBC spectrum. To include this in our analysis of the stability, we add random perturbations with the only two constraints that they should be Hermitian and physically motivated. To be precise, we change all off-diagonal zero matrix elements of H to H ij = H * ji = (a + ib) exp(−|i − j| 2 ) where a, b are random numbers ∈ [−2, 2] and thus are of the same order as the nearest neighbor hopping amplitudes in our examples in Figure 5. We hereby show that the open eigenspectrum as well as the edge states are robust against disorder that is physically motivated (see Figure  5a) and c)). To underline the importance of the locality of perturbations, we demonstrate that the edge states are not robust against perturbations that decay ∼ exp(−|i − j|), i.e. comparable to the scaling of Γ c . Specifically, adding perturbations that scale as exp(−|i−j|·ln(Γ P )/N ), where the factor ln(Γ P )/N ensures that especially the perturbation term between the first and the last site is of the same order as the sensitivity of the eigenspectrum. We find, that boundary modes disappear for such perturbations(see Figure 5b) and d)).

VI. DISCUSSION
We have studied the sensitivity of the eigen-spectrum of NH tight-binding models towards their boundary conditions. In particular, we derived analytical expressions showing how exact zero-energy modes move through parameter space as a function of the boundary condition parameter Γ interpolating between PBC and OBC. In several NH models, we find that this eigenmode at E = 0 represents an EP that can switch from embodying a bulkgap closing point to being an isolated edge mode. These two cases are found to be distinguished by the biorthogonal polarization P introduced in Ref. [54]. Corroborating earlier observations, we analytically quantify how at a critical value Γ c ∝ exp(−αN ) [52,54,60] a transition occurs between the qualitative spectral properties found for PBC and the ones for OBC.
We have further discussed the stability of surface states arising in these models towards physically motivated perturbations. In summary, the spectral instability of NH matrices loses its dramatic appearance when imposing physical locality conditions on the considered perturbations. Specifically, in a flat sample geometry, accidental couplings between opposite ends are found to scale as ∝ exp(−λN 2 ) in generic tight-binding models, and thus naturally decay even faster with system size than the critical coupling Γ c . Moreover, the topological zero-energy edge modes in NH systems are found to be robust towards generic local perturbations such as disorder potentials. Hence, despite the sensitivity of NH eigen-spectra to the choice of boundary conditions, the NH BBC as characterized by the biorthogonal polarization P and the concomitant edge modes may be seen as a topologically stable and experimentally observable phenomenon.
Finally, we note that some experiments on classical systems such as optical meta-materials, despite a certain similarity to quantum-mechanical tight-binding Hamiltonians, exhibit couplings between sites that do not follow a Gaussian (sometimes not even exponential) decay with distance. Investigating possible modifications to the definition of topological invariants and the stability of the NH BBC in such long-ranged settings is an interesting subject of future research. consisting of the Hamiltonians H atom (r−R n ) of the single atoms at positions R n and a correction potential V 1 (r). After a Bloch-wave ansatz for the eigenfunctions one arrives at with v n = d 3 rV 1 (r)|ψ n (r)| 2 (29) γ (j) n = d 3 rψ * n (r)V 1 (r)ψ n (r − R j ) (30) α (j) n = d 3 rψ * n (r)ψ n (r − R j ) where E n is the nth eigenenergy and ψ n the nth eigenfunction of the single atom Hamiltonian. γ n ) translate furthermore into the matrix elements of a tightbinding Hamiltonian in second quantization. Usually, these integrals become very small and sums in (28) are truncated for small j already. For our purposes, we are interested in the scale describing the decay of these neglected matrix elements. We approximate the single atoms with harmonic potentials since close to the equilibrium point it is a justified assumption for any potential and the eigenfunctions of the quantum oscillator are well known. Furthermore we assume V 1 to be constant so that we effectively only need to study the overlap integral α (j) n . Here, we use the eigenfunctions of the one-dimensional harmonic oscillator, noting that approximating with a 3D harmonic potential leads to the same asymptotic behavior. The overlap integral between two sites at positions R i = 0 and R j = ja, where a is the lattice spacing, in one spatial dimension is given by with eigenfunctions [73, p. 312-315] ψ n (x) = 1 √ 2 n n! mω π 1/4 e − mω x 2 2 · H n mω x (33) and H n (z) are Hermite polynomials H n (z) = (−1) n e z 2 d n dz n e −z 2 .
After the typical variable transformation z = mω π x and shifting the integrand to a symmetric form we find H n (z + ja/2)H n (z − ja/2).
(35) We want to examine the case of large distances j and only consider the first couple of eigenstates with small n. Hence, we are interested in the leading order of the Hermite polynomials H n (z) in terms of ja/2 which is 2 n z n and thus, the leading term of the integrant will be 2 2n (z 2 − j 2 a 2 /4) n . Splitting the integrant in summands with different power of z, we see that the zeroth power of z is of leading order and explicitly reads The matrix element H ij thus scales with ∝ exp(−|i−j| 2 ).