Sensitivity of non-Hermitian systems

Understanding the extreme sensitivity of the eigenvalues of non-Hermitian Hamiltonians to the boundary conditions is of great importance when analyzing non-Hermitian systems, as it appears generically and is intimately connected to the skin effect and the breakdown of the conventional bulk boundary correspondence. Here we describe a method to find the eigenvalues of one-dimensional one-band models with arbitrary boundary conditions. We use this method on several systems to find analytical expressions for the eigenvalues, which give us conditions on the parameter values in the system for when we can expect the spectrum to be insensitive to a change in boundary conditions. By stacking one-dimensional chains, we use the derived results to find corresponding conditions for insensitivity for some two-dimensional systems with periodic boundary conditions in one direction. This would be hard by using other methods to detect skin effect, such as the winding of the determinant of the Bloch Hamiltonian. Finally, we use these results to make predictions about the (dis)appearance of the skin effect in purely two-dimensional systems with open boundary conditions in both directions.


I. INTRODUCTION
In the past few years, the interest in non-Hermitian Hamiltonians, and in particular the topological properties of such systems , has become huge. The differences between non-Hermitian systems and their Hermitian counterparts are many and give rise to interesting effects.
Mathematically, these differences arise from the fact that a non-Hermitian operator can lose several of the key properties that make Hermitian operators so nice to work with. Most importantly, the eigenvalues are no longer necessarily real, the right and left eigenvectors might not be related to each other via Hermitian conjugation or be orthogonal amongst themselves. Also, the operator might not be diagonalizable, leading to the existence of exceptional points. This makes computing probabilities and expectation values non-trivial, something which is discussed in e.g. [25], where the biorthogonal inner product is defined and used to replace the standard inner product in probability calculations.
In particular, there are several new phenomena that occur in non-Hermitian tight-binding models compared to Hermitian ones. Perhaps most notably, the bulkboundary correspondence breaks down [5,7,[13][14][15][16][17][18], and is replaced by a biorthogonal bulk-boundary correspondence as described in [5,14]. The breakdown of the bulk-boundary correspondence is related to the extreme sensitivity of eigenvalues and eigenvectors to boundary conditions that often occurs in non-Hermitian systems. In fact, in a system with open boundary conditions, the left and right eigenstates tend to be exponentially localized to one of the boundaries of the system, and the number of such localized states is extensive. This phenomenon is called the skin effect [13]. The breakdown of the bulk-boundary correspondence and the skin effect are intimately connected to each other and have lead to an active field of research [26][27][28][29][30][31][32][33][34][35][36][37]. In particular, the phenomena have been experimentally verified in mechanical systems [38][39][40][41], topoelectrical circuits [42,43], photonic systems [19,44], and in ultracold atoms [45]. This phenomenology has also been suggested to be of practical use in sensors whose sensitivity increases exponentially with the size of the system [46][47][48].
In one-dimensional systems it is clear what is meant by this description of the skin effect, but in higher dimensions a little more care is needed in defining what is meant as states can be localized to different boundaries of the system in different ways. Nevertheless, the fact that most of the eigenstates are exponentially localized to some boundary of the open system, gives an intuitive explanation of why the spectrum and the eigenstates should change radically when coupling the endsif the boundary is removed, the states cannot be localized to it anymore and the exponential localization must disappear. Because the eigenstates drastically change when coupling the edges, the eigenvalues also undergo a large shift.
Now, it turns out that not all non-Hermitian systems exhibit a skin effect, and one important question is how to determine which systems do and which not. In a onedimensional system, described by the Bloch Hamiltonian H(k), it turns out that there is a relation between the winding number around some base energy E B ∈ C and the existence of the skin effect in the corresponding system with open boundary conditions. Namely, if there exists an E B ∈ C such that w(E B ) = 0, then the system with open boundary conditions has a skin effect [2,8], otherwise not. This also implies that if we have a non-zero winding number for some E B , the spectrum of the system will be very sensitive to boundary conditions. The winding number is related to the types of gaps that can be found in a non-Hermitian system [8,12]. Namely, if w(E B ) is nonzero for some E B , we say that the system exhibits a point gap. This shows up as loops in the spectrum of the Bloch Hamiltonian. On the other hand, if w(E B ) = 0 for all E B , we say we have a line gap, which typically implies that the spectrum of the Bloch Hamiltonian consists of line segments that do not form loops.
In [14], it was numerically found that if the ends of an SSH-chain are coupled by a parameter δ, the δ required to change the spectrum by a fixed amount ∆ is proportional to e −ξ(∆)N , where N is the length of the chain and ξ(∆) > 0 also depends on system parameters. This kind of sensitivity was also studied in [7], and we will call it exponential sensitivity to boundary conditions, because an exponentially small change in the boundary condition δ leads to a finite change ∆ in the spectrum. If we do not have this kind of sensitivity, we say the spectrum has non-exponential sensitivity. Clearly, when numerically calculating the spectrum of non-Hermitian operators, one should be careful if the system has exponential sensitivity to the boundary conditions. For larger system sizes, even a machine-precision deviation can lead to substantial errors in the eigenvalues obtained. It is thus advantageous to know in advance when the system is exponentially sensitive and when it is not. We note that algorithms calculating the spectrum of non-Hermitian PT-symmetric systems to arbitrary precision we considered in [49,50].
In this paper, we aim to study the sensitivity of the eigenvalues analytically for different systems described by tight-binding Hamiltonians of the form where, c † m and c n are creation and annihilation operators respectively and the t mn are hopping parameters for which t mn is not necessarily equal to t * nm . We begin with purely one-dimensional tight-binding models where we interpolate between open and periodic boundary conditions using a parameter δ. For a class of such systems, we develop a method to analytically find how the eigenvalues depend on the parameter δ and the system size. Then we move on to two-dimensional systems, constructed by stacking one-dimensional chains. In a rectangular geometry, the Hamiltonian can be represented by a block-tridiagonal block-Toeplitz matrix, with extra blocks in the corners to account for the boundary conditions. Such matrices can in general not be diagonalized exactly, unless we have periodic boundary conditions, but using the result from the one-dimensional case, we can still draw some conclusions about these systems.
From the analytical expressions for the eigenvalues, we can see for which parameter values the spectrum should be exponentially sensitive and for which not. It turns out that in order for the spectrum to have non-exponential sensitivity there has to be some kind of balancing of the hopping parameters; somewhat loosely we need the amount of hopping to the right to be balanced by a similar amount to the left. Because of this, we will say that a system with parameter values such that the spectrum has non-exponential sensitivity is balanced. We find that the parameter values for which the system is balanced are in agreement with what the winding number predicts about the skin effect. For a model, with the hopping parameters explicitly specified, one can often easily obtain the winding number by looking at the plot of H(k) in the complex plane. However, finding out for which parameters the model has (or does not have) a non-trivial winding, is often much more complicated. In those cases, it might be much easier to determine if the system is balanced or not.
The paper is organized as follows. In Sec. II, we introduce the method to analyze the non-Hermitian systems we are interested in, and apply it to several models, including the Hatano-Nelson model and the SSH-chain. In Sec. III, we use solved one-dimensional systems, and construct several two-dimensional models, including the triangular lattice, which shows particularly interesting behavior. We discuss the results in Sec. IV. The appendices provide more details on some of the models we studied.

II. ONE-DIMENSIONAL SYSTEMS
We begin by describing a method to find the eigenvalues of a one-dimensional, one-band, non-Hermitian system. Any such system with N sites and a maximum hopping range of m < N/2 , can be described by an N × N Toeplitz matrix. Now, in order to simplify the notation, we will here only show the method for a system with next-nearest neighbour hopping, but it is straightforwardly generalizable to include longer range hopping. The Hamiltonian of a one-dimensional system with nextnearest neighbour hopping and the coupling strength between the ends determined by the parameter δ can be represented by the following matrix: For δ = 1, we have a system with periodic boundary conditions and for δ = 0, we have a system with open boundary conditions. In the case of periodic boundary conditions, the eigenvalues are found by Fourier transforming the system, and are given by where k = 0, 1, . . . , N − 1 and ω = exp(2πi/N ). Now, we make the ansatz that the eigenvalues for the case with general δ are given by λ α = t 11 + t 12 e iα + t 13 e 2iα + t 21 e −iα + t 31 e −2iα , (5) where α is a complex number, the possible values of which need to be determined.
From the Schrödinger equation, Hψ α = λ α ψ α , we get a system of equations containing N − 4 bulk equations of the form t 31 ψ α,n−2 + t 21 ψ α,n−1 + (t 11 − λ α )ψ α,n + t 12 ψ α,n+1 + t 13 ψ α,n+2 = 0 (6) with n = 3, . . . , N − 2, and four boundary equations The strategy to find the eigenvalues is to first find the general solution of Eq. (6) for arbitrary integer n, and then to determine for which values of α there exists a non-trivial solution to the boundary equations in Eq. (7). Because we solve the bulk equation in Eq. (6) for arbitrary integer n, we can use it to simplify the boundary equations, which gives Eq. (6) is a linear recurrence relation for the elements of the vector ψ α and has a characteristic polynomial of the form for t 13 = 0. (In case t 13 = 0, we get a polynomial of lower degree.) The zeros of this polynomial, x i , give us the general solution of the recurrence relation, namely where the constants c i are determined by inserting the expression for ψ α,n into the boundary equations in Eq. (8). This gives us a linear system of equations for the constants of the form where M (α) is a 4×4-matrix and C is the vector containing the constants c i . To find the values of α for which this system has a non-trivial solution, we compute the determinant of the matrix M (α) and solve the equation The form of this equation can tell us how sensitive α, and in turn the eigenvalues λ α are to changes in the boundary conditions.
In general, the bulk and boundary equations cannot be solved analytically as this involves finding roots of polynomials of high degree, but for some systems the method provides useful information. We point out that similar methods were used in, for instance, [51,52] in the special cases of the Hatano-Nelson model and the SSHchain. In these papers, an explicit ansatz was used for the eigenvectors, instead of using the recurrence relation.
To analyze the results, we are interested in the properties of the eigenvectors. For non-Hermitian systems, the left and right eigenvectors are generically not orthonormal, and one can therefore consider different ways in which to normalize them. The eigenvectors considered above, ψ α , with components ψ α,n , are the right eigenvectors, and we denote those by ψ r α and ψ r α,n , respectively, if confusion with the left eigenvectors can arise. The left eigenvectors, ψ l α with components ψ l α,n , are obtained by diagonalizing H T instead of H. In the cases we consider below, this means a simple swap of parameters of the model.
Following [5,14,25], we consider different expectation values of the site projection operator Π n = c † n c n . On the one hand we have the left and right expectation values Π n l,l α = (ψ l α ) † Π n ψ l α = |ψ l α,n | 2 and Π n r,r α = (ψ r α ) † Π n ψ r α = |ψ r α,n | 2 , which can be used to determine if the system exhibits skin effect. Namely, the right and left expectation values tell us where in the system the eigenvectors are localized, which, as explained in the introduction tells us if the skin effect is present or not. In the case of one-dimensional systems, as we consider in this section, skin effect simply means that an extensive number of right eigenvectors are localized on one side of the system, while the left eigenvectors are localized near the other.
The third kind of expectation value we consider, is the biorthogonal expectation value, which, as is shown in [25] is the one most similar to what we call an expectation value in Hermitian quantum mechanics. The biorthogonal expectation value is given by Π n l,r α = (ψ l α ) † Π n ψ r α = (ψ l α,n ) * ψ r α,n and can be used to distinguish between bulk and boundary states in the system [5,14].

A. The Hatano-Nelson model
We begin with a simple, non-Hermitian nearestneighbour hopping model, also known as the Hatano-Nelson model [53][54][55]. The Hamiltonian for a chain of length N with boundary conditions regulated by the parameter δ is given by To analyze the situation where we interpolate between open and periodic boundary conditions, it is easier to write the Hamiltonian in matrix form. Then, it is given by the N × N matrix where δ = 1 corresponds to periodic boundary conditions and δ = 0 corresponds to open boundary conditions. Below, we only give the Hamiltonians in matrix form as in Eq. (14). Often, we consider the case 0 ≤ δ ≤ 1, but unless otherwise stated, the formulas are valid for arbitrary (complex) δ. In Appendix A 2, we derive the results for more general boundary conditions, where we allow for different coupling parameters in different directions.
In accordance with Eq. (5), we get the following ansatz for the eigenvalues in the case with arbitrary boundary conditions: where α is some complex number. The recurrence relation for the elements in the eigenvector is now given by t r ψ α,n−1 + t l ψ α,n+1 − (t l e iα + t r e −iα )ψ α,n = 0, (16) which has the solution where r = t r /t l . We see here that as long as c 2 = 0 and |r| = 1, we have an exponential localization of all the eigenstates to one of the ends of the chain, depending on the absolute values of t r and t l . To find out when it is possible to have c 2 = 0 (that is, when we do not have an exponential localization of the eigenstates near the edges), we use the boundary equations, which in this case are given by By setting c 2 = 0, we can use Eq. (18) to determine the values of δ for which we get N inequivalent solutions for α, giving a complete set of eigenvalues. This procedure leads to two possible values of δ, namely δ = ±1. For other values of δ, we have c 2 = 0. Below, we give the explicit eigenstates for completeness, even though we do not need them to determine if the eigenstates are exponentially localized near the edges or not.
We continue by determining the equation for α, such that we obtain a non-trivial solution for the coefficients c 1 and c 2 . It turns out that it is convenient to write the eigenvalues in terms of The eigenvalues are then given by while the eigenvectors read In this form, one finds, provided the solutions forα are real, that the eigenvalues lie on a straight line in the complex plane. Similar variable changes, though in slightly different settings, for example in the context of the generalized Brillouin zone, were considered in [13,18]. The possible values ofα correspond to the solutions of the equation which is found by computing the determinant obtained from the boundary equations, as explaind above. We note that the form of the eigenvalues in Eq.(20) is similar, up to values ofα, to the eigenvalues of a tridiagonal Toeplitz matrix. For completeness, we give a more explicit form of the (right) eigenvectors for arbitrary δ, though still in terms of the parameterα. Solving the boundary Eq. (18) results in The left eigenvectors can be obtained by diagonalizing H T , which means swapping t l ↔ t r .
We have now obtained expressions for the eigenvalues and eigenvectors, such that we can investigate the their behavior when we change δ.

Interpolating the boundary conditions
We start by noting that as long as |t l | = |t r |, the equation forα, Eq. (22), contains a term proportional to δ that is exponential in the system size N . Therefore, the spectrum has an exponential sensitivity to the boundary conditions when |t r | = |t l |. When, on the other hand, |t l | = |t r |, the spectrum is non-exponentially sensitive to boundary conditions. We note that this happens when there seems to be an overall balancing of the hoppings in the system in the different directions, and therefore, as mentioned in the introduction, we will say that the hoppings are balanced when this happens. In fact, when this is the case, r is equal to a phase, say e iθ , which means that Eq. (22) can be written as 2δ cos(θ) sin(α) = sin([N + 1]α) − δ 2 sin([N − 1]α). (25) We show in Appendix A 1 that for δ ∈ [−1, 1] and θ real, this equation only has real solutions forα. This implies that Eq. (20) describes a line segment in the complex plane and that the system does not have a point gap at the balanced parameter values. Therefore, the system does not exhibit a skin effect when it is balanced, which is in correspondence with the fact that the spectrum has non-exponential sensitivity to boundary conditions at these parameter values.
We show this difference in behavior between balanced and unbalanced systems when going from periodic to open boundary conditions in sites, and parameters t l = 1 and t r = 2e iπ/4 such that |t l | = |t r |. In the lower left panel, we show the same results, but for the parameters t l = 1 and t r = e iπ/4 instead, such that |t l | = |t r |. In the former case, we see a 'jump' in the eigenvalues between δ = 0 and δ = 1/100, indicating a significant change in the spectrum, while in the latter case, such a jump is absent and the eigenvalues stay on the same line in the complex plane, as explained above.
In the right panel, we plot, for both sets of parameters, the left and right expectation values, Πα ,n l,l = |ψ lα ,n | 2 and Πα ,n r,r = |ψ r α,n | 2 , for a representative eigenstate, showing a pronounced skin effect in the upper panel, while in the lower panel, the skin effect is absent.
The large change in the eigenvalues upon a small change in δ also shows up in the eigenvectors. In Fig. 2, we plot both the (logarithm of the) right expectation values Πα ,n r,r = |ψ r α,n | 2 , as well as all the biorthogonal expectation values Πα ,n l,r = (ψ lα ,n ) * ψ r α,n . The right eigenvectors clearly show the presence of the skin effect, because they are exponentially localised on the right hand side of the system. For δ = 0, there is a clear oscillatory behavior, which quickly disappears upon increasing δ, as follows from Eq. (24).
In addition, the correlation length changes by several orders of magnitude upon changing δ = 0 to δ = 1/100, while it remains almost the same upon changing δ = 1/100 to δ = 2/100. The left eigenvectors show similar behavior, but are localised on the left hand side of the system.
The biorthogonal expectation value also shows the exponential sensitivity to small deviations of δ from zero. In this case, the weight of the eigenvectors is spread out over the bulk of the system (so we do not have edge states), but the way in which this occurs changes drastically with δ. For δ = 0, the inner product oscillates as a function of position (see Eq. (24)), while upon increasing δ, these oscillations disappear in the bulk, where the biorthogonal inner product becomes constant.
In conclusion, we see a perfect correspondence between the exponential localisation of eigenstates to one side of the chain, the sensitivity of the eigenvalues to the boundary conditions and the existence of point gaps in the system. In particular, we see that when |t r | = |t l |, we have a non-Hermitian system whose eigenvalues and eigenstates have properties that closely resemble what we have in a Hermitian system in terms of sensitivity to boundary conditions and localization respectively.

Application to impurity-like systems
We can also consider the case δ > 1. In this case, we do not think of δ as a parameter interpolating between open and periodic boundary conditions, but rather as a parameter that sets the strength of an impurity in the system, namely a single link with enhanced hopping.
For δ > 1, we observe the following behavior: Both for balanced |t l | = |t r | as well as unbalanced parameters, |t l | = |t r |, there is one impurity state, signified by an exponential localization of the biorthogonal expectation around the impurity [60]. The remaining states are bulk states.

B. The SSH-chain
A system with more interesting properties than the Hatano-Nelson model, is the SSH-chain. It is described by the N × N -matrix where δ ∈ [0, 1] and interpolates between open and periodic boundary conditions and N determines the length of the chain. Here N is assumed to be even since we want to be able to interpolate between open and periodic boundary conditions, but the case N odd can be dealt with in a similar way, and the result is stated in Appendix B.
It has previously been observed [13,14] that there is a pronounced skin effect for most parameter values, and that the eigenvalues are highly sensitive to changes in boundary conditions. This can be seen in the upper panel of Fig. 3, where we to the left plotted the eigenvalues for a generic system for different values of δ, and to the right an example of what the left and right eigenstates look like (for more details on this figure, see below).
In this section we study this sensitivity in more detail. Even though this system is not described by a Toeplitz matrix, and the method described in the beginning of this section at first glance thus seems to be unusable in this case, we notice that the square of the Hamiltonian is given by This is almost a Toeplitz matrix with changes only in the parts of the matrix that would affect the boundary equations, so we can use the previously described method for H 2 0 , and instead solve the problem H 2 0 ψ α = λ α ψ α , with the appropriate boundary equations. The eigenvalues of H 2 0 will be the squares of the eigenvalues of H 0 (which come in (λ, −λ) pairs when N is even). Using the eigenvalues of the periodic SSH-chain together with the bulk equations of H 2 0 , we get the following ansatz for the eigenvalues of H 2 0 and general δ: The bulk equations now give us the following recurrence relation for the elements of ψ α : This is solved, for arbitrary integer n, by where the constants c i are determined by the boundary equations As for the Hatano-Nelson model, it turns out that it is convenient to apply a shift to α Then, the eigenvalues of H are given by ± t l,1 t r,1 + t l,2 t r,2 + 2 cos(α) t l,1 t l,2 t r,1 t r,2 , whereα follows from the determinant equation, which takes the form We find that unless the spectrum has exponential sensitivity. In this special case, we have for some θ ∈ [0, 1]. As was the case for the HN-model, hermiticity is not required to get eigenvalues that are insensitive to the boundary conditions. In Fig. 3, we illustrate these results. In the upper left panel, we plot the eigenvalues of a system with N = 30 sites, and hopping parameters t l,1 = 1, t r,1 = 2, t l,2 = 3, t r,2 = 4, with δ varying from δ = 0 to δ = 1. We see that the system has a point gap, and that the spectrum is exponentially sensitive to the variation in δ. In the lower left panel of Fig. 3, we plot the case t l,1 = −i, t r,1 = 1/2, t l,2 = 4, t r,2 = 8, for which, on the contrary, the eigenvalues show a nonexponential sensitivity to the variation in δ. In this case, the system has a line gap.
In the rightmost panel of Fig. 3, we plot, for both sets of parameters, the left and right expectation values, Πα ,n l,l = |ψ lα ,n | 2 and Πα ,n r,r = |ψ r α,n | 2 , for a representative eigenstate, showing a pronounced skin effect in the upper panel, while in the lower panel, the skin effect is absent.
As was the case for the HN-model, we observe a strong skin effect in the unbalanced case, which also exhibits a point gap, and a strong δ dependence of the eigenvalues, in accordance with results in previous works.
For both sets of parameters in Fig. 3, we observe a double zero-mode E ≈ 0 for δ = 0, which gradually merges with the rest of the spectrum upon varying δ from δ = 0 to δ = 1. Even in the case t l,1 = 1, t r,1 = 2, t l,2 = 3, t r,2 = 4, where the rest of the spectrum is exponentially sensitive to varying δ away from zero, the energy of the δ = 0 zero mode changes gradually with δ. This corresponds to what was seen in e.g. [46], where it was shown that the zero mode can be insensitive to boundary conditions even though the rest of the spectrum is not. In appendix B, we obtain the condition for having a zero mode, which happens for t l,2 tr,2 t l,1 tr,1 > 1 in the large N limit. This agrees with what was found in [56] and generalizes the condition found in e.g. [14]. We note that in the balancing condition, for which the spectrum is insensitive to changes in δ, the left and right hopping parameters are 'paired up', while in the condition for having a zero mode, the hopping parameters with index one and index two are 'paired up'.
To make the system more useful for later, we add alternating on-site potentials according to To find the eigenvalues, we note that where and Since we can directly use the result from Eq. (28) and add v 2 . The final result, then, is whereα is again given by Eq. (34), and we see that alternating on-site potentials do not change the sensitivity of the eigenvalues to boundary conditions.

C. Longer range hopping
Both the Hatano-Nelson model and the SSH-chain are examples of models with nearest neighbor hopping. In this section we study some examples of models with longer range hopping and see if they can be balanced. In general, longer range hopping makes the recurrence relation in Eq. (45) more complicated, but when one restricts to the case of two parameters, the equations simplify.

Unidirectional hoppings
First, we consider a system with nearest and nextnearest neighbour hopping to the left, but no hopping to the right. This is described by the matrix We note that for δ = 0, the matrix only has one eigenvector and all its eigenvalues are zero. For δ = 0, we make the ansatz The recurrence relation for the elements of ψ α becomes and is solved by The constants c 1 and c 2 are again found using the boundary equations Dropping an unimportant overall factor, this gives us the following determinant equation for α: One can use either factor to determine α, both give rise to the same eigenvalues. By using the first factor, one finds that the solutions for α only depend on δ, not on the other parameters of the model. Parametrising δ = |δ|e iφ , we obtain the following eigenvalues, with j = 0, 1, . . . , N − 1 which means that if δ is real and δ ∈ (0, 1), the solutions will approach the solutions of the periodic case for large N . In the special case that δ = 0, all eigenvalues will be 0 regardless of N . This implies an extreme sensitivity to the boundary conditions for all values of t l and u l , which makes sense, since there is no way to balance the hopping to the left with hopping in the opposite direction. In Fig. 4, we plot the eigenvalues of a chain of length N = 50 with t l = 1 and u l = 2 for values of δ between 0 and 1, and see that even for very small values of δ, the eigenvalues deviate a lot from 0 and approach the eigenvalues of the periodic case.

Hoppings in different directions
Next, we study the Hamiltonian given by the matrix which describes a system with hoppings both to the right and to the left, but with different range. We make the following ansatz for the eigenvalues: The recurrence relation for the elements of ψ α becomes which has the characteristic equation This equation has roots such that the elements of ψ α are given by where We note that The constants c 1 , c 2 , c 3 are determined using the boundary equations, which in this case read The values of α for which this system of equations has solutions correspond to when the determinant of the system equals zero. The expression for the determinant is long, but can be simplified a lot using Eq. (57). Introducing the polynomial p(n), recursively defined by with p(0) = 0 and p(1) = −1, the equation for α can be written in the following way: where we introduced y = e iα and we dropped a factor (x − − x + ). Though not manifestly so, the equation is proportional to (r − 2y 3 ) (we note that r = t r /(u l y 3 )). After dividing out by this factor, the remaining polynomial equation has degree 3N in terms of y, and degree N in terms of t r /u l . Thus, we find 3N solutions for y (and We note that the change in the eigenvalues when we go from δ = 0 to δ = 0.01 is much larger when tr = u l . There is, however, still a significant change in the spectrum also for tr = u l . hence α), but they come in triples which give rise to the same eigenvalue, so we obtain N eigenvalues as wanted.
For δ = 0, the equation reduces to (61) and for δ = 1, we have (up to unimportant factors) We note that Eq. (60) contains terms with different powers of r. Since r depends on t r /u l , one could expect the behavior of the eigenvalues to be insensitive to perturbations when |u l | = |t r |, but this is not the case for this model, in contrast to previous models. We can, however, see that the winding number of the Bloch Hamiltonian is non-zero when |u l | = |t r |, so the result is somewhat expected. Technically, the reason for this behavior is that the eigenvectors Eq. (55) generically have an exponential behavior, even when the parameters appear to be 'balanced' as |u l | = |t r |.
We illustrate this in Fig. 5, where we plot the eigenvalues for a chain of length N = 30, with parameters (u l , t r ) = (2, 1); (1, 1); (1, 2), and δ varying from δ = 0 to δ = 1. We observe that in all three cases, the eigenvalues are sensitive to the boundary conditions, and there is a non-zero winding number for δ = 1, even in the case where the hopping to the left and right seem to be balanced. Also the eigenstates have a pronounced skin effect (not shown). 3. Example of a system not solvable by the method Since the described method to find eigenvalues relies on finding zeros of a polynomial equation of a degree proportional to the hopping range, it will typically fail for more complicated systems since we cannot solve such equations exactly. Nevertheless, we will briefly look at one such system, which will be of interest for later use.
We show this model in Fig. 6. The periodic version of this chain is described by the Bloch Hamiltonian We first note that an interesting special case of the model described by Eq. (63) is the case u l = t r and u r = t l . For these parameters, one ends up with a chain of triangles and a Bloch Hamiltonian given by H(k) = t l e ik + t r e −ik + t r e 2ik + t l e −2ik = cos 3k 2 t r e ik/2 + t l e −ik/2 .
Unless t r /t l is a phase, the eigenvalues in general form a loop in the complex plane, and the system will thus be exponentially sensitive to boundary conditions. We note that there are several choices of the parameters, such that H(k) in Eq. (63) does not wind.
One way of achieving this is, is by setting t l = te i(φ+φ1/2) , t r = te i(φ−φ1/2) , u l = ue i(φ+φ2/2) , u r = ue i(φ−φ2/2) . For these parameters, the Hamiltonian becomes Assuming that t, u and the angles φ are real, this forms a straight line segment in the complex plane, and thus there is no winding of the spectrum in this case, which means that for these parameter values the spectrum has non-exponential sensitivity. We note that the condition we get on the hopping parameters for this to be the case, is slightly different from previous conditions in systems with nearest neighbor hopping. Namely, we get a conditions also on the relative phases of the parameters, and not only on the absolute values.
A second way of achieving this, is the case t l = te iφ1 u l = te iφ1+iφ t r = te iφ2 u R = te iφ2−iφ . For these parameters, the Hamiltonian becomes (66) Again, for t and the phases φ real, this forms a straight line segment in the complex plane, implying that the spectrum has non-exponential sensitivity.
(67) In this case, the values of H(k) (with t, u and the phases φ real) do not generically form a straight line in the complex plane. However, because the argument of the second cosine is twice the argument of the first cosine, one obtains a curve where the part with π ≤ k ≤ 2π traces back the part of the curve with 0 ≤ k ≤ π. Again, we conclude that there is no winding for these parameters, and hence the spectrum has non-exponential sensitivity.
We believe that the three cases above exhaust the ways in which one can obtain an H(k) that does not wind.

III. TOWARDS TWO-DIMENSIONAL SYSTEMS
In this section, we consider two-dimensional systems that we construct by stacking one-dimensional chains studied in Sec. II. Contrary to the one-dimensional case, we will here allow for more complicated boundary conditions, so in the direction of stacking, we will have two parameters, δ 2 and δ 2 determining the boundary conditions. The one-dimensional chains, however, have the regular coupling by a single parameter δ 1 , which takes values in the interval [0, 1]. Suppose we stack N 2 chains, each containing N 1 sites. The Hamiltonian for such a system can be described by the where A(δ 1 ), B(δ 1 ) and C(δ 1 ) are N 1 × N 1 -matrices describing the one-dimensional chains and the coupling between them respectively (we drop the dependence of A, B, and C on the other parameters of the model) [61]. To find analytic expressions for this kind of matrix for arbitrary δ 1 and δ 2 , δ 2 is in general a very hard problem.
In particular it is difficult to find the eigenvalues for open boundary conditions in both directions, i.e. for all δ's being zero. If one were to try to implement the method used for one-dimensional systems, by considering a unit cell of size N 2 , one would end up with several bulk equations that depend on the size of the system, and this would not be solvable. Alternatively, one can simplify the bulk equations, by making use of the two-dimensional nature of the problem, and write the elements of the eigenvectors as ψ i,j . This indeed leads to a simple bulk equation, but also results in the complication that the boundary equations in one direction should be satisfied for all possible locations in the other direction. It turns out that this approach does not work either.
Below, we consider two special types of boundary conditions, for which we can obtain an analytic solution, but we start even simpler, by briefly noting that one can construct solvable models on the square lattice (for details, see App. C). To do this, we simply take a one-dimensional model that we solved, i.e., for which we know the functional form of the eigenvalues λ α , as well as the equation for α, denoted by eq α . The square lattice model now has the hopping parameters of this model in the horizontal direction, and the hopping parameters of a second solved model in the vertical direction, with eigenvalues λ β , where β satisfies eq β . The eigenvalues of the square lattice model are then simply given by λ α,β = λ α + λ β , since the two directions are completely independent of each other.
We will, however, be interested in studying more complicated models that cannot be solved exactly if we have open boundary conditions. Thus, in this section, we will study two special cases of boundary conditions that can be solved analytically and then use them to try and make predictions of the system with open boundary conditions. BC 1. We will begin by assuming that δ 2 = δ 2 = 1, making the system periodic in one direction, while we can interpolate between open and periodic boundary conditions in the other direction. Taking the Fourier transform of H in the periodic direction, we get a Bloch Hamilto-nian of the form where now ω j = exp(2πij/N 2 ), which has the property that if λ jk is an eigenvalue ofH j with corresponding eigenvector v jk , then λ jk is an eigenvalue of H with corresponding eigenvector BC 2. Next, we consider the case δ 2 = δ −1 2 . Define the matrix Then (71) This describes a system with periodic boundary conditions in one direction, and its eigenvalues are given by the eigenvalues of the matrices where ω j = exp(2πij/N 2 ). For δ 2 = 1, this case reduces to BC 1 as it should.

A. Stacking Hatano-Nelson chains
Now we study an explicit model, which is insensitive to boundary conditions, for correctly chosen parameters. Consider the lattice shown in Fig. 7 with boundary conditions according to BC 1. In this system, we have and The eigenvalues of the matrix H can be found by diagonalizing the matrixH j defined in Eq. (69). In this case, we haveH where The matrixH j is precisely of the form discussed in Eq. (14) and its eigenvalues are thus given by whereα is determined by the equation We note that for the spectrum has non-exponential sensitivity. This condition is independent of t d , u d and u u , which is reasonable since these parameters are either on-site potentials or hoppings in the direction of periodic boundary conditions. In case Eq. (80) is fulfilled, as was previously shown,α is real and thus cos(α) ∈ [−1, 1]. This means that the eigenvalues given in Eq. (78) are contained in the set given by the values of the function where and with c ∈ [−1, 1] and t ∈ [0, 2π]. We obtained z 1 (t) and z 2 (t) from h d and 2 √ h l √ h r respectively by replacing ω j with e it . We note that for each t, the values of f (t, c) form a line segment in the complex plane between the points z ± (t) = z 1 (t) ± z 2 (t). That is, where t ∈ [0, 2π]. This means that irrespective of the value of δ 1 , the eigenvalues of H will be contained in the area formed by connecting the points z ± (t) on the two curves by straight lines. We also note that as we increase N 1 , the number of different solutions to the equation for α will increase, and as we increase N 2 the number of ω j increases, and thus, as we make the system larger, more and more of the area will be covered by eigenvalues. Therefore, as we will see is useful later, we argue that in the thermodynamic limit, the spectrum in this area will be dense. We are interested in finding out when the parameter values fulfill Eq. (80). This happens when h l = e 2πir/N2 h r , which means that where r is a real number. Here we will study some particularly interesting cases. From now on we will assume that the parameters of H are real.
Case 1. For h r = h * l , Eq. (80) is fulfilled and we have r = − arg(h r )N 2 /π. Furthermore, we get t r = t l , v ur = v dl and v dr = v ul . We note that the Hermitian case is a special case of this. In the upper panel of Fig. 8, we plot the curves z ± (t) together with the eigenvalues of a system with t r = t l = 2, v ur = v dl = 4, v dr = v ul = 3, t d = 1, u u = −3 and u d = 2 for N 1 = N 2 = 30 when we change δ 1 in steps of 0.01 from δ 1 = 0 (blue dots) to δ 1 = 1 (red dots). We see that the eigenvalues lie along straight line segments connecting the two curves and that the spectrum has non-exponential sensitivity. In the lower panel, we plot the right expectation value Πα ,n r,r = |ψ r α,n | 2 for a representative eigenvector, which shows that we have no skin effect at these parameter values.
Case 2. The case r = 0 implies h r = h l , i.e. t l = t r , v ul = v ur , v dl = v dr , which means that the three matrices A, B and C are Hermitian, while the full Hamiltonian H might not be. The eigenvalues are in this case given by whereα is given by the equation In this case, z ± (t) reduces to where t ∈ [0, 2π]. That is, two ellipses centered at t d ±2t r .
In the upper panel of Fig. 8, we plot the curves z ± (t) together with the eigenvalues of a system with t r = t l = 2, v ur = v ul = 3, v dr = v dl = 4, t d = 1, u u = −3 and u d = 2 for N 1 = N 2 = 30 when we change δ 1 in steps of 0.01 from δ 1 = 0 (blue dots) to δ 1 = 1 (red dots). In the lower panel, we plot the right expectation value Πα ,n r,r = |ψ r α,n | 2 for a representative eigenvector, which shows that we have no skin effect at these parameter values.
Case 3. The case r = j (i.e., e 2πir/N2 = ω j ) implies that t l = v ur , v dl = t r , v ul = v dr = 0. This is a non-Hermitian matrix with eigenvalues given by whereα is given by the equation Rewriting the eigenvalues as shows that z ± (t) in this case reduce to where t ∈ [0, 2π]. These two curves are actually two segments that together form a closed loop s(t ) in the complex plane, described by where t ∈ [0, 2π]. The eigenvalues of H will thus lie on line segments joining two points on this curve.
In the upper panel of Fig. 8, we plot the curves z ± (t) together with the eigenvalues of a system with t r = v dl = 1, t l = v ur = 2, v dr = v ul = 0, t d = 1, u u = −3 and u d = 2 for N 1 = N 2 = 30 when we change δ 1 in steps of 0.01 from δ 1 = 0 (blue dots) to δ 1 = 1 (red dots). In the lower panel, we plot the right expectation value Πα ,n r,r = |ψ r α,n | 2 for a representative eigenvector, which shows that we have no skin effect at these parameter values.
Case 4. The case r = −j (i.e., e 2πir/N2 = ω −1 j ) implies that t l = v dr , v dl = v ur = 0, v ul = t r . This case is similar to the previous one, and therefore not shown in Fig. 8.
Unbalanced. Finally, we consider an unbalanced case, with parameters t d = 1, t r = 2, t l = 3, u u = 4, v ur = 5, v ul = 6, u d = 7, v dr = 8, v dl = 9. From the plot in the upper panel of Fig. 8, where we have plotted the eigenvalues of this system for N 1 = N 2 = 30 when changing δ 1 in steps of 0.01 from δ 1 = 0 (blue dots) to δ 1 = 1 (red dots), it is clear that the eigenvalues do not lie on straight lines between two curves. In the lower panel, we plot the right expectation value Πα ,n r,r = |ψ r α,n | 2 for a representative eigenvector, which shows that in this case, we do have a skin effect.
Comparing the plots in Fig. 8, we see a clear difference in behavior of the eigenvalues between the balanced and unbalanced cases. Similar to what we saw for the Hatano-Nelson model and the SSH-chain, we see that the eigenvalues of the unbalanced system change significantly even for a very small change in δ 1 , while there is a much less drastic change for the balanced system.

Special case: A triangular lattice
Up until now, we have effectively been studying a onedimensional system with a large unit cell. We can, however, attempt to make some predictions of the behavior of the eigenvalues of the system with non-periodic boundary conditions in both directions. Namely, let us implement the balancing conditions from Case 3 in both directions of the lattice. Then, we end up with and which describes a triangular lattice, schematically shown in Fig. 9. For this subclass of Case 3, it turns out that if one implements the boundary conditions BC 1, the eigenvalues are actually contained within the curve Eq. (93). One would expect that having implemented the balancing condition in both directions, the eigenvalues would show non-exponential sensitivity to changes in δ 1 and δ 2 . Surprisingly enough, however, we find this to be the case only when N 1 = N 2 . We numerically study the triangular lattice with open boundary conditions in one direction when going from open to periodic boundary conditions in the other. As expected, the eigenvalues still lie within the region bounded by the curve in Eq. (93) for all values of N 1 and N 2 . Moreover, the spectrum seems to have non-exponential sensitivity when N 1 = N 2 . However, when N 1 = N 2 , and in particular when N 2 N 1 or N 1 N 2 , we start to see a higher sensitivity to boundary conditions and that the eigenstates start to localize at the end of the lattice. This can be seen in Fig. 10, where we plot a part of the eigenvalue spectrum for a triangular lattice with t l = 1 and t r = 5 for different values of N 2 when keeping δ 2 = δ 2 = 0 and changing δ 1 from 0 to 1 in steps of 0.01. In this figure, we also plot the right expectation value Πα ,n r,r = |ψ r α,n | 2 for a representative eigenvector, for the different values of N 2 . We see how the states show a significant localization to one of the sides of the system when N 2 N 1 , but not when N 2 = N 1 . We already encountered this behavior in Sec. II C 3 for the Hamiltonian in Eq.(64), which describes the case N 2 = 2, where we saw that we do have a skin effect. Here, we thus see a gradual disappearance of the skin effect when we let N 2 approach N 1 . We checked that the localization of the eigenstates is always 'in the long direction of the system', so the localization due to the skin effect we observe here, is similar to the localization due to skin effect in the one-dimensional systems we studied earlier. To explain the presence of the skin effect in the triangular lattice model, we interpret the model as a one-dimensional model with a large unit cell of size N 2 and length N 1 , such that N 2 N 1 . In this way, we can use the phase winding formula (1), which implicitly assumes the thermodynamic limit N 1 → ∞. To apply Eq. (1), we assume, for simplicity, that the system is open in the '2' direction. That means that the Bloch hamiltonian is a N 2 × N 2 tridiagonal matrix where a = t l e −ik +t r e ik , b = t r +t l e ik , c = t l +t r e −ik . We denote the determinant by det(N 2 ). It is straightforward to show that det(N 2 ) = a det(N 2 − 1) − bc det(N 2 − 2), while det(1) = a and det(2) = a 2 − bc. Generically, det(N 2 ) will wind around some point in the complex plane, implying that the system shows a skin effect, consistent with our results that we for generic parameters observe a skin effect when N 2 = N 1 . The argument above does not imply that we should have a skin effect for all values N 1 , because we need the assumption that N 1 N 2 . Indeed, when N 1 is of the order of N 2 , the system cannot be considered one-dimensional. Nevertheless, it is interesting to see that we do not observe a skin effect when N 1 ∼ N 2 and both large, because in principle, there could have been.
Using the above approach, we can also try to find parameters for which the skin effect is absent, regardless of the system sizes. It turns out that by setting t r /t l = e iφ , with φ real, the determinant takes the form det(N 2 ) = e iN2φ/2 f (k, φ), where f (k, φ) is real and periodic in k. So in this case, det(N 2 ) does not wind in the complex plane, and we do not have a skin effect. This is consistent with what we saw in Sec. II C 3.
The dependence of the skin effect on the ratio N 1 /N 2 we observe here seems to contradict the results in [57]. Namely, in this reference, it is claimed that we have a The dots show how the eigenvalues change when we change δ1 in steps of 0.01 from 0 (blue dots) to 1 (red dots). Lower panel: plot of the right expectation value of a representative state for each of the respective system sizes. We see the gradual disappearance of the skin effect as N2 approaches N1 when the right eigenstates become delocalized and the eigenvalues become less sensitive to boundary conditions. skin effect in a two-dimensional system if and only if the spectral area is finite. As we have argued before, the spectral area of the triangular lattice is finite when we let the system size go to infinity, also when keeping the ratio N 1 /N 2 = 1 constant. This does also not seem to be taken into account by the exception of the geometry dependent skin effect, described in the same reference, where it is said that if the boundary of the system coincides with a mirror symmetry line of the lattice, the skin effect will disappear, even if the spectral area is finite. In our case, the boundaries of the system do coincide with mirror symmetry lines, but we still have, for N 1 = N 2 , a localization of the eigenstates to one of the sides of the lattice, which implies that the mirror symmetry is not enough to prevent the skin effect from appearing.
We do not have analytical expressions for the eigenvalues in the case that δ 2 = δ 2 = 0, and δ 1 interpolates between 0 and 1, so we cannot see analytically how the difference between N 1 and N 2 would affect the behavior of the eigenvalues in this case. However, we can get a hint of this if we instead consider the system with the different boundary conditions from BC 2. Then we get where We know the eigenvalues of this matrix. They are given by whereα is determined by the equation and we see that unless N 1 = N 2 there will be an exponential dependence on system size inα. To derive Eq. (101), we assumed that N 1 is even, so that we do not have to worry about square roots. We note that this exponential dependence on the system size is somewhat different form the exponential dependence we encountered so far.
Here, it is the perturbing parameter δ 2 that is raised to the power N 1 /N 2 , while earlier, it was some combination of the hopping parameters that was raised to the system size.
As we stated above, the boundary conditions we considered here are not the same as in Fig. 10. However, with the boundary conditions BC 2, we do get similar behavior of the spectrum as in Fig. 10. Namely, depending on the ratio N 1 /N 2 , there is a skin effect (that is pronounced when N 1 N 2 ). In the presence of the skin effect, the eigenvalues are also sensitive to changes in δ 2 . Because of the similarity in behavior, we think that the mechanism at play in the case of Fig. 10, is similar to the mechanism we just obtained analytically for the boundary conditions BC 2.
We close the discussion on the triangular lattice system, by making a few remarks about the sensitivity of the eigenvalues when we do have a skin effect. Because the exponential dependence on the system size in Eq. (101) has exponent N 1 /N 2 , one has to specify in which way one takes the large system size limit, when discussing the sensitivity of the eigenvalues as a function of system size. If one fixes the ratio N 1 /N 2 , and take both N 1 and N 2 large, the sensitivity of the eigenvalues does not change in the large system size limit. If, on the other hand, one fixes the size of one of the directions, the exponential sensitivity increases when making the other direction large. In any case, if the eigenvalues are sensitive, this is due to the presence of the skin effect. As stated above, the interesting aspects of the triangular lattice system are the presence of the skin effect, despite the fact that we have 'balanced' parameters in both directions, and that the appearance of the skin effect depends on the ratio N 1 /N 2 .

The Kagome lattice
Perhaps a more interesting system, which, as is shown in [5], supports boundary states, is the Kagome lattice. This system has similarities with the triangular lattice, and the results from the triangular case seem to be generalizable to the Kagome bulk states. Here we look at the Kagome lattice with hoppings according to Fig. 11.
Also in this case, we observed numerically that for N 1 = N 2 , the bulk spectrum has non-exponential sensitivity to the boundary conditions, while for N 2 N 1 , we seem to have a higher sensitivity among the bulk states. By studying the right eigenvectors of these insensitive eigenvalues, we notice that, just as in the case of the triangular lattice, for each eigenvector ψ, Π n rr becomes more localized the bigger the difference between N 1 and N 2 is, which implies that we get a skin effect when N 1 and N 2 deviate from each other. This phenomenon with the emerging skin effect shows that the explanation for the absence of the skin effect given in [5] is incomplete. There, the closed loops formed in the system were used as an explanation for not getting this build-up of the eigenstates in the system, but since the loops exist also in the case of N 1 = N 2 there must be something more behind this. We do currently not have an explanation for this behavior, but as we observed for the triangular lattice in the case of BC2, the sensitivity depends on N 1 /N 2 , which is at least an indication that it is not unreasonable to expect a similar dependence on system size also in this case.

B. Stacking SSH-chains
Finally, we study the lattice in Fig. 12, which is obtained by stacking SSH-chains. In this case, we have (for N 1 even) which means that we get where and ω j = e 2πij/N2 . From Eqs. (42) and (34) in section II B, we know how to find the eigenvalues of this matrix. Namely, they are given by whereα is determined by the equation As before, we note that if the exponential sensitivity of the spectrum disappears. Similar to the results in the previous section, we see that this gives rise to several different cases, and we will now list some of them together with plots of the eigenvalues of some representative systems. Assuming that all parameters are real, we first list four cases where individual hoppings are related to each other: The case h r,1 = h l,1 and h r,2 = h l,2 . This implies We notice that all these cases imply an overall balancing of the hoppings in the unit cell. In general, however, Eq. (109) implies where r ∈ R, which means that we also can get conditions consisting of relations between products of hopping parameters. We list a few such cases.
While the cases 1-4 give us equalities between individual hopping parameters, we see that cases 5-7 give us relations between products of adjacent hopping parameters. We also note that in the cases 5-7 we still need some kind of balance in the left-right direction, even though it is not as straightforward as in previous cases. It is interesting to see that by stacking SSH chains, there are indeed more complicated ways in which the system can be balanced, because of the fact that one has more parameters to play with. We note, just as in the case with stacked HN-chains, that the hopping parameters in the direction of the periodic boundary conditions do not affect the sensitivity of the spectrum at all.
From the plots in Fig. 13, we see a significant difference in the behavior of the spectrum between the balanced and a representative unbalanced case when we change the boundary conditions.
We would like to comment on Case 4 in Fig. 13. This case does not look 'as neat' at the Cases 1-3. One could wonder if this is due to numerical inaccuracy, despite the fact that this is a balanced case. We therefore checked the numerical results against our analytical results, and found perfect agreement between the two (that is, agreement up to machine precision).

IV. DISCUSSION
In this paper, we obtain (almost) analytical expressions for the eigenvalues of non-Hermitian one-dimensional one-band models with arbitrary boundary conditions. We use the analytical expressions to analyze the sensitivity of the eigenvalues to the boundary conditions in several one-dimensional systems. Using the one-dimensional models, we construct two-dimensional models with arbitrary boundary conditions in one direction, and either periodic or a particular deformation of periodic boundary conditions in the other direction.
We find that for most parameter values, we have eigenvalues that are exponentially sensitive to boundary conditions, but that there are parameter values for which this exponential sensitivity disappears. In this case, the behavior of the eigenvalues, apart from the fact that the eigenvalues can still be complex, is similar to what we expect in a Hermitian system.
In one-dimensional chains with nearest neighbor hoppings, we find the rather intuitive result that in order for the spectrum to behave similar to the Hermitian case, the hoppings to the left should balance those to the right. This seems to be a necessary, but not sufficient condition for 'Hermitian behavior', which is a bit surprising. Even more surprising is the fact that in the case of longer range hoppings it might not be possible at all to balance the system. This indicates that the intuitive understanding of why the system shows 'Hermitian behavior' in some cases is lacking, and it would be interesting to understand why this is.
In two dimensions, the situation is more complicated. Ideally, one would analytically like to study arbitrary boundary conditions in both directions, but this turns out to be out of reach. One can, however, study the case with open boundary conditions in one direction, and a deformation of periodic boundary conditions in the other. One can then use these results to explain the numerically obtained behavior of two-dimensional systems with open boundaries in both directions.
In particular, we studied a rather general two-dimensional model obtained by stacking HN-chains. This model contains a hopping model on the triangular lattice as a special case. By studying this model, with periodic boundary conditions in one direction, and arbitrary boundary conditions in the other, we obtain a condition on the hopping parameters, such that the system is not exponentially sensitive. This condition amounts to a local balancing of the hopping in the direction with arbitrary boundary conditions, as in the one-dimensional case. If one implements this balancing condition in both directions, one could expect that the system is not exponentially sensitive. Our numerical results show that this is however not the case. Only if the system size is the same in both directions, the spectrum has nonexponential sensitivity. We obtained similar results for the Kagome lattice. Clearly, obtaining a full understanding of the behavior and sensitivity of the eigenvalues of two-dimensional systems with open or arbitrary boundary conditions in both directions is a question of great interest for future investigation. So far, in contrast to the one-dimensional case where the winding number can be used, there is no general condition for the existence of the skin effect in two-dimensional systems, even though there have been attempts to come up with such conditions [57], but these conditions fail to capture the dependence on the aspect ratio of the system. In addition, our results also indicate that the picture put forward in [5], namely that "loops" in the system would prevent the skin effect, is incomplete.
Apart from stacking HN-chains to obtain twodimensional models, we also studied two-dimensional models obtained by stacking non-Hermitian SSH-chains. In this case, the balancing conditions become more complicated, in comparison to the case where we stacked HN-chains. Namely, the balancing condition for the SSH case, becomes an equation for sums of product of hopping parameters. Despite the more complicated condition, one can still find non-trivial solutions, and indeed cases for which the system is balanced, thus without the exponential sensitivity to the boundary conditions.
In the current paper, we considered the eigenvalue stability of non-Hermitian hopping models. It would certainly be interesting to extend the current analysis to models including pairing terms. To deal with pairing terms, one could use the method of Lieb, Schultz and Mattis [58], and apply it to the non-Hermitian case. One can speculate that it should be possible to find balanced systems, even in the presence of pairing terms. Council (VR).
Appendix A: Some details for the NH model In this appendix, we give some more details of the non-Hermitian hopping model we consider in Sec. II A. Namely, we show thatα is real under certain conditions, and we provide the solution for a slightly more general model.

Proof thatα is real for the HN model under certain conditions
In this section, we show that Eq. (25), repeated here for convenience, only has real solutions forα, provided that −1 ≤ δ ≤ 1 and θ ∈ R, that is, for |t l | = |t r |. Under these conditions, the eigenvalues of the Hatano-Nelson model lie on a linesegment in the complex plane. We recall that the trivial solutionsα = 0, π should be discarded (unless there are multiple solutions at these values), because they do not lead to valid eigenstates. From the argument below, it will become clear that we can assume that cos(θ) = ±1, because | cos(θ)| < 1 follows easily. For concreteness, we assume that cos(θ) = 1.
We find that the independent solutions forα are indeed real, and given byα = pπ/N , with p even for δ = 1 and p odd for δ = −1. Discarding the solutions coming from sin(α) = 0, we find that all solutions are double solutions. To pick an independent set of solutions, we take the double solutions with 0 < p < N , but only one solution for p = 0 and p = N .
The proof works similarly when −1 ≤ cos(θ) < 1, as well for −1 ≤ δ < 0, which means we covered all cases, and we are done.

Solution for a slightly more general hopping model
In this section, we provide the solution for a slightly more general non-Hermitian hopping model. In particular, we consider the following hamiltonian (here given for N = 5 sites) H =      1 t l 0 0 δ r t r t r 0 t l 0 0 0 t r 0 t l 0 0 0 t r 0 t l δ l t l 0 0 t r N      . (A2) The method to obtain the eigenvalues is identical to the one used in the main text in Sec. II A, so we simply state the result here. The eigenvalues of this modified non-Hermitian hopping model are given by whereα satisfies the transcendental equation For δ l = δ r = 0, this case was studied in [59].
Appendix B: Some details on the SSH chain In this appendix, we give some more details on the non-Hermitian SSH model.

Solution for even length chains
For completeness, we state the solution for the eigenvalues of the non-Hermitian SSH chain, Eq. (26), but in the presence of the perturbations δ r and δ l in the upper right and lower left corner respectively (i.e., with matrix elements δ r t r,2 and δ l t l,2 respectively).

Solution for odd length chains
In the main text, we give the solution of the system for even system sizes. Here we provide the solution for chains of odd length. In this case, the most natural perturbing elements are δ l √ t l,1 √ t l,2 and δ r √ t r,1 √ t r,2 . That is, the hamiltonian takes the form (B3) The method to obtain the eigenvalues is the same as for the chains of even length, that is, one determines λ 2 α by considering H 2 . The functional form of λα is given in Eq. (B1). The equation to determineα is now given by λ 2 α sin(α(N + 1)/2) sin(α) − δ l δ r sin(α(N − 1)/2) sin(α) 2 = δ 2 r t r,1 t r,2 t r,1 t r,2 t l,1 t l,2 (N −1)/2 + δ 2 l t l,1 t l,2 t l,1 t l,2 t r,1 t r,2 (N −1)/2 + 2δ l δ r t l,1 t r,1 t l,2 t r,2 .
In terms of x = e iα , this equation has degree 2L (after multiplying with x L ), with solutions coming in (x, 1/x) pairs, both members of a pair giving rise to the same value of λ 2 α . In contrast with the case N even, the spectrum is not symmetric, so for each independent solutioñ α, one needs to determine the actual sign of the eigenvalue.
We note that the evaluation of σα can be somewhat unstable numerically, for given values ofα, in particular when δ l , δ r are large in comparison to the hopping parameters in the model. zero eigenvalue, and determine for which parameters the equation forα is satisfied, up to corrections that are exponentially small in system size.