Mathematical Foundations of the Non-Hermitian Skin Effect

We study the skin effect in a one-dimensional system of finitely many subwavelength resonators with a non-Hermitian imaginary gauge potential. Using Toeplitz matrix theory, we prove the condensation of bulk eigenmodes at one of the edges of the system. By introducing a generalised (complex) Brillouin zone, we can compute spectral bands of the associated infinitely periodic structure and prove that this is the limit of the spectra of the finite structures with arbitrarily large size. Finally, we contrast the non-Hermitian systems with imaginary gauge potentials considered here with systems where the non-Hermiticity arises due to complex material parameters, showing that the two systems are fundamentally distinct.


Introduction
Understanding the skin effect in non-Hermitian physical systems has been one of the hottest research topics in recent years [21,22,26,30,32,36].The skin effect is the phenomenon whereby the bulk eigenmodes of a non-Hermitian system are all localised at one edge of an open chain of resonators.This phenomenon is unique to non-Hermitian systems with non-reciprocal coupling.It has been realised experimentally in topological photonics, phononics, and other condensed matter systems [16,19,23,31].By introducing an imaginary gauge potential, all the bulk eigenmodes condensate in the same direction [20,28,32].This has deep implications for the fundamental physics of the system.For example, the non-Hermitian skin effect means that the conventional bulk-edge correspondence principle is violated.
In this paper, we study the non-Hermitian skin effect in a deep subwavelength regime using first-principles mathematical analysis.The ultimate goal of subwavelength wave physics is to manipulate waves at subwavelength scales.Recent breakthroughs, such as the emergence of the field of metamaterials, have allowed us to do this in a way that is robust and that beats traditional diffraction limits.Here, we consider one-dimensional systems of highcontrast subwavelength resonators as a demonstrative setting to develop a mathematical and numerical framework for the non-Hermitian skin effect in the subwavelength regime.
We consider finite chains of subwavelength resonators with an imaginary gauge potential supported inside the resonators.This imaginary gauge potential is realised by adding an additional first-order term to the classical Helmholtz differential problem (the γ-term in (2.4)).It is worth emphasizing that the γ-term is added only inside the subwavelength resonators, which are much smaller than the operating wavelength.Without exciting the structure's subwavelength resonances, the effect of the γ-term would be negligible.Using an asymptotic methodology that was previously developed for Hermitian systems [5,15], we can approximate the subwavelength eigenfrequencies and eigenmodes of a finite chain of resonators by the eigenvalues of a matrix, which we call the gauge capacitance matrix.When γ = 0, this capacitance matrix reduces to the capacitance matrix previously derived for Hermitian systems of subwavelength resonators [2,15].This traditional capacitance matrix was originally developed to model the analogous many-body problem in the setting of electrostatics [12,24].Given this matrix formulation, we can use the well-established theory of tridiagonal Toeplitz matrices and operators to prove the condensation of eigenmodes at one of the edges of the finite chain.
Using established spectral theory for tridiagonal Toeplitz matrices, we will show the topological origins of the non-Hermitian skin effect.That is, we will demonstrate the condensation of the eigenmodes.As a consequence of the condensation of the eigenmodes, the imaginary gauge potential leads to a mismatch between the eigenvalues of a finite chain and the band structure of the associated infinite periodic structure, when the band structure is computed over the standard Brillouin zone (as used for Hermitian problems).However, the correspondence between the spectra of finite and infinite structures can be understood if the spectral bands of the infinite periodic structure are computed over the complex plane, in the sense that the quasiperiodicities must be allowed to be complex [10].The eigenmodes behave as Bloch waves with a complex wave number, with the imaginary part describing the exponential spatial decay of the amplitude.We will show in Theorem 4.8 that the generalised Brillouin zone introduced in (4.8) gives the correct description of the limiting problem.
Although the spectrum of the infinite periodic problem computed over the standard Brillouin zone does not provide the correct approximation of the spectrum of a finite resonator chain, the problem is nonetheless insightful to investigate in its own right.We will consider exceptional points where eigenvalues and eigenvectors coincide and the quasiperiodic gauge capacitance matrix is not diagonalisable [7] and prove their existence for periodic dimer structures at the edges of the Brillouin zone.We will also show that the set of exceptional points in terms of the separating distance between the two identical resonators (which constitutes the dimer) in the unit cell separates two different phases described by different values of the winding number of the eigenvalues, which is known as the vorticity.
The imaginary gauge potential considered here is not the only way to break the Hermiticity of a system of subwavelength resonators.This is also achieved by making the material parameters complex valued, corresponding to sources of energy gain or loss in the system.An important class are PT-symmetric systems, which satisfy a combination of time-reversal and inversion symmetry.They inspired a plethora of explorations in topological photonics and phononics and their mathematical foundations [7,13,17,27].Based on the contrast with the previous work [7,9], it is natural to divide non-Hermitian systems into two classes depending on whether they can be reduced to Hermitian systems through an imaginary gauge transformation.If such reduction is possible (which we examine in Appendix D), the eigenmodes will be closely approximated by Bloch modes over the standard Brillouin zone, and the vorticity will vanish.This means, on the one hand, that vorticity can be considered as a topological indicator of the appearance of the skin effect and, secondly, that the spectral convergence of large structures with vanishing vorticity can be understood using the standard approaches for Hermitian systems, based on computations over the standard Brillouin zone [4].Conversely, in the setting considered in this work, this convergence does not hold as the system is fundamentally distinct from any Hermitian system.The current work is the first to establish the convergence theory when such reduction is not possible and the eigenmodes are described through the generalised (complex) Brillouin zone.
The paper is organised as follows.In Section 2 we present the mathematical setup of the problem and derive an asymptotic approximation of the subwavelength eigenfrequencies and their associated eigenmodes.These are accurately approximated by the spectrum of the gauge capacitance matrix.Section 3 is dedicated to the mathematical analysis of the skin effect.The condensation of the eigenmodes is proved and the non-uniform distribution of the subwavelength eigenfrequencies is illustrated.Section 4 studies the infinite periodic case.In Section 4.1 we first take the quasiperiodicities to be real.We show that using the standard Brillouin zone does not lead to the correct modeling of the limit of the set of subwavelength eigenfrequencies as the size of the system goes to infinity.Nevertheless, the band functions over the standard Brillouin zone are used to define the vorticity, which we show that it is non-trivial in the current setting.Such topological invariant can be used as an indicator of the skin effect since it is trivial for other non-Hermitian structures obtained by taking the material parameters inside the subwavelength resonators complex.In Section 4.2, we introduce the generalised (complex) Brillouin zone and prove the convergence of the subwavelength eigenfrequencies of a finite system to the bands of the corresponding infinitely periodic structure computed over the generalised Brillouin zone.In Appendix A we summarise the well-established theory of Toeplitz matrices and operators, while in Appendix B we present some additional numerical results for finite dimer systems with imaginary gauge potentials.In Appendix C we show how to obtain a generalised discrete Brillouin zone for different dimer systems.Finally, in Appendix D we show that systems with non-Hermiticity arising from complex material parameters can be described in the limit when their sizes go to infinity by bands of corresponding infinitely periodic structures computed over the standard Brillouin zone.

Preliminaries
We consider a one-dimensional chain of N disjoint subwavelength resonators for any 1 ≤ i ≤ N .We fix the coordinates such that x L 1 = 0. We also denote by the length of the i-th resonators, and by s i = x L i+1 − x R i the spacing between the i-th and (i + 1)-th resonator.The system is illustrated in Figure 1.We will use to symbolise the set of subwavelength resonators.
As a scalar wave field u(t, x) propagates in a heterogeneous medium, it satisfies the following one-dimensional wave equation: The parameters κ(x) and ρ(x) are the material parameters of the medium.We consider piecewise constant material parameters where the constants ρ b , ρ, κ, κ b ∈ R >0 .The wave speeds inside the resonators D and inside the background medium R \D, are denoted respectively by v b and v, the wave numbers respectively by k b and k, and the contrast between the densities of the resonators and the background medium by δ: Up to using a Fourier decomposition in time, we can assume that the total wave field u(t, x) is time-harmonic: for a function u(x) which solves the one-dimensional Helmholtz equations: In this work, we will consider the following variation of (2.2): for a piecewise constant coefficient This new parameter γ extends the usual scalar wave equation to a generalised Strum-Liouville equation via the introduction of an imaginary gauge potential [33].Alternatively, one may think about (2.3) as a damped wave equation where the damping acts in the space dimension instead of the time dimension.In these circumstances of piecewise constant material parameters, the wave problem determined by (2.3) can be rewritten as the following system of coupled one-dimensional equations: (2.4) We are interested in the resonances ω ∈ C such that (2.4) has a non-trivial solution.We will perform an asymptotic analysis in a high-contrast limit, given by δ → 0. We will look for the subwavelength modes within this regime, which we characterise by satisfying ω → 0 as δ → 0. This limit will recover subwavelength resonances, while keeping the size of the resonators fixed.
We can immediately see that ω = 0 is a trivial solution to (2.4) corresponding to an eigenmode which is constant across all the resonators.In what follows, we will restrict attention to other solutions of (2.4).
The use of the Dirichlet-to-Neumann map has proven itself to be an effective tool to find solutions to (2.4) that lend themselves to asymptotic analysis in the high-contrast regime.Hereafter, we will adapt the methods of [15] to (2.4).
The first step is to solve the outer problem for some boundary data f L,R i ∈ C 2N .Its solution is simply given by where a i and b i are given by as shown in [15,Lemma 2.1].The second and harder step is to combine the expression for the outer solution (2.7) with the boundary conditions in order to find the full solution.We will handle this information through the Dirichlet-to-Neumann map.Definition 2.2 (Dirichlet-to-Neumann map).For any k ∈ C which is not of the form nπ/s i for some n ∈ Z \{0} and 1 ≤ i ≤ N − 1, the Dirichlet-to-Neumann map with wave number k is the linear operator where w is the unique solution to (2.6).
We refer to [15, Section 2] for a more extensive discussion of this operator, but recall that T k has a block-diagonal matrix representation where, for any real ∈ R, A k ( ) denotes the 2 × 2 symmetric matrix given by (2.9) Consequently, T k is holomorphic in k in a neighbourhood of the origin and admits a power series representation T k = n≥0 k n T n , where and A 0 (s) := lim k→0 A k (s).The above properties of the Dirichlet-to-Neumann map will be crucial to find subwavelength eigenfrequencies.We will use T k and T k interchangeably.

Characterisation of the subwavelength eigenfrequencies
The Dirichlet-to-Neumann map allows us to reformulate (2.4) as follows: (2.11) We can then further rewrite (2.11) in weak form by multiplying it by a test function w ∈ H 1 (D) and integrating on the intervals.Explicitly, we obtain that u is solution to (2.11) if and only if a(u, w) = 0 (2.12) for any w ∈ H 1 (D), where We also introduce a slightly modified bilinear form This new form parametrised by ω and δ is an analytic perturbation of the a 0,0 form, which is continuous coercive on H 1 (D) and so a ω,δ inherits this property.We exploit this, by defining the h j (ω, δ) functions as the Lax-Milgram solutions to the variational problem for every w ∈ H 1 (D) and 1 ≤ j ≤ N .In the following lemma we show that the functions h j allow us to reduce (2.12) to a finite dimensional N × N linear system by acting as basis functions.
Lemma 2.3.Let ω ∈ C and δ ∈ R belong to a neighbourhood of zero such that a ω,δ is coercive.The variational problem (2.12) admits a non-trivial solution u ≡ u(ω, δ) if and only if the N × N non-linear eigenvalue problem has a solution ω and x := (x i (ω, δ)) 1≤i≤N , where C(ω, δ) is the matrix given by (2.16) Proof.This can be shown using the arguments presented in [15,Lemma 3.4].

Discrete approximation of the subwavelength eigenfrequencies and the gauge capacitance matrix
In the high-contrast, low-frequency scattering problems, recent work has shown that capacitance matrices can be used to provide asymptotic characterisations of the resonant modes [5].Capacitance matrices were first introduced to describe many-body electrostatic systems by Maxwell [24], and have recently enjoyed a resurgence in subwavelength physics.
In particular, for one-dimensional models like the one considered here, a capacitance matrix formulation has proven to be very effective in both efficiently solving the problem as well as providing valuable insights in the solution [2,15].Therefore, we seek a similar formulation for the case of subwavelength resonators with an imaginary gauge potential.
In the case of systems of Hermitian subwavelength resonators (i.e., when γ = 0), the entries of the capacitance matrix are defined by where ν is the outward-pointing normal, D i the i-th resonator and V i : R → R the solutions of the problems where δ ij denotes the Kronecker symbol; see [5].
Here, the formulation for the non-Hermitian system (2.4) is different from the Hermitian case.The following definition is in agreement with the three-dimensional case for the skin effect [1].We let R * := R \{0}.
Definition 2.4 (Gauge capacitance matrix).For γ ∈ R * , we define the gauge capacitance matrix where V j is defined by (2.18).
It is easy to see that the gauge capacitance matrix is tridiagonal, non-symmetric, and is given by while all the other entries are zero.
The following proposition is a central result.It shows that to leading order in δ the gauge capacitance matrix encodes all the information of subwavelength eigenfrequencies.Proposition 2.5.Let ω ∈ C and δ ∈ R belong to a small enough neighbourhood of zero.Then, the matrix C(ω, δ) from (2.16) has the following asymptotic expansion:

.21)
Here, is the volume matrix.
Proof.We need to show that unique solution h j (ω, δ) with 1 ≤ j ≤ N to the variational problem (2.12) has the given integral asymptotic behaviour as ω, δ → 0. From the definition of a ω,δ , the function h j ≡ h α j (ω, δ) satisfies the following differential equation written in strong form: (2.23) Since T ω v is analytic in ω 2 , it follows that h j (ω, δ) is analytic in ω 2 and δ: there exist functions (h j,2p,k ) p≥0,k≥0 such that h j (ω, δ) can be written as the following convergent series in H 1 (D): (2.24) By using the power series expansion of the Dirichlet-to-Neumann map and identifying powers of ω and δ, we obtain the following equations characterising the functions (h j,2p,k ) p≥0,k≥0 : with the convention that h j,2p,k = 0 for negative indices p and k.It can then be easily obtained by induction that for any p ≥ 0, 1 ≤ j ≤ N.
Then, for p = 0 and k = 1, we find that h j,0,1 satisfies (2.26) (2.27) The general solution for the ODE −y − ay + b = 0 is given by for two constants κ 1 and κ 2 .So the solution to (2.26) on 1 Imposing the transmission conditions on the derivatives, we obtain a value for the integral which combined with (2.27) yields the result after some careful algebraic manipulation.
The following corollary, whose proof is similar to the analogous results in [2,14], describes how the gauge capacitance matrix characterises the non-trivial solutions to (2.4).
Corollary 2.6 (Discrete approximations of the eigenfrequencies and eigenmodes).The N subwavelength eigenfrequencies ω i satisfy, as δ → 0, where (λ i ) 1≤i≤N are the eigenvalues of the eigenvalue problem (2.28) Furthermore, let u i be a subwavelength eigenmode corresponding to ω i and let a i be the corresponding eigenvector of C γ .Then where V j are defined by (2.18) and a (j) denotes the j-th entry of the eigenvector.
Lemma 2.7 (Properties of the gauge capacitance matrix).
(i) Recall the Hermitian capacitance matrix C, defined in (2.17).Then, for any (ii) For equally spaced identical resonators: s i = s, i = for all 1 ≤ i ≤ N , the gauge capacitance matrix is almost Toeplitz taking the form (2.29)

Skin effect and localised interface modes
We now turn to the question of eigenmode condensation.We will begin by proving that the skin effect occurs in a finite system of subwavelength resonators.Then, once this fundamental result has been established, we will consider more exotic structures.For example, we will consider a structure with an interface formed by adjoining two half-structures with opposite signs of γ, leading to wave localisation along the interface, as shown Section 3.2.

Skin effect
Using the extensive theory of Toeplitz matrices and operators, we will now show that the system (2.4) is a prototypical example of the non-Hermitian skin effect.We will consider a system of equally spaced identical resonators as described in Lemma 2.7 (iii).That is, a chain of N resonators with s i = s and i = for all 1 ≤ i ≤ N .For this case, the particular structure of the capacitance matrix allows us to apply existing results concerning the spectra of tridiagonal Toeplitz matrices, which we briefly recall in Appendix A. We will be able to derive a variety of explicit results that would not be possible for more complex situations.Nevertheless, Corollary 2.6 applies for more general cases and the numerical results (which we present for dimers in Appendix B) show that similar results hold for systems with multiple resonators in the unit cell.
Applying the explicit formula for the eigenpairs of perturbed tridiagonal Toeplitz matrices from Lemma A.6 to C γ , we obtain the following theorem.
Theorem 3.1.The eigenvalues of C γ are given by Furthermore, the associated eigenvectors a k satisfy the following inequality, for for some Here, a k denotes the i-th component of the k-th eigenvector.It is easy to show that the first eigenvector a 1 is a constant vector (i.e. a (i) 1 is the same for all i).From Theorem 3.1, we can see that the eigenmodes display exponential decay both with respect to the site index i and the factor γ. We show this decay in a simulation of a large system of resonators in Figure 2. Furthermore, (3.2) shows that changing the sign of γ swaps the edge at which the localisation occurs.
Of particular interest is the application of Lemmas A.2, A.3 and A.4 to the gauge capacitance matrix.The combination of these lemmas shows that the localisation of the eigenvectors of both the finite and semi-infinite capacitance matrix depends on the Fredholm index of the symbol of the associated Toeplitz operator at the corresponding eigenvalue or, equivalently, of its winding number.Specifically, every eigenvector having the corresponding eigenvalue laying the region of the complex plane where the Fredholm index (or equivalently, the winding number of the symbol) is negative, has to be localised.This shows the intrinsic topological nature of the skin effect.
In Figure 3 we display the convergence of the pseudospectrum and the topologically protected region of negative winding number.We remark that the trivial eigenvalue 0 is outside of the region as the winding number is not defined there.As a consequence, the corresponding eigenvector is not localised.
The explicit formula for eigenvalues from Lemma A.6 also gives some insight on the distribution of the density of states.The property of the cosine to have minimal slope when 0.16 0.17  it is close to its maxima and minima, causes eigenvalues to cluster at the edges of the range of the spectrum.This is demonstrated by the non-uniform distribution of the black dots in Figure 3.
We conclude this section with a qualitative analysis of the spectral decomposition of C γ .In Figure 4a we show the eigenvector localisation as a function of the site index for finite arrays of various sizes.The localisation of a vector v is measured using the quantity v ∞ / v 2 .After rescaling the site index, we expect there to be some invariance to the array size N , based on the formulas of Lemma A.6. Figure 4b shows, on the other hand, the singular values of the eigenvector matrix, again with the indices normalised.The number of non-zero singular values is a proxy for the dimension of the range of the matrix.This has an exponential dependence on N .As a result, Figure 4 shows that as N increases, despite the number of eigenvalues growing like N , the rank of the matrix of eigenvectors is fixed.As a result, we have an increasing number of collinear eigenvectors.This can be interpreted as there existing an exceptional point of "infinite" order in the arbitrarily large system.As N → ∞, the finite structure exhibits an exceptional point of "infinite" order: arbitrarily many eigenmodes become close to parallel.Here shown in the case s = = 1 and γ = 0.5.

Non-Hermitian interface modes between opposing signs of γ
We now consider interface modes between two structures where the sign of γ is switched from negative (on the left part) to positive (on the right part).Most commonly, localised interface modes are formed by creating a defect in the system's geometric periodicity (see, for example, [8]).In non-Hermitian systems based on complex material parameters, similar localised interface modes have been shown to exist in the presence of a defect in the periodicity of the material parameters [2,9].Given the existence of the skin effect, demonstrated in the previous section, it is reasonable to expect that we might be able to produce a similar localisation effect using systems of resonators with imaginary gauge potential.With this in mind, we consider the following system of N = 2n + 1 resonators: It is not difficult to see that also with this system we can recover a capacitance matrix for which a similar result as Corollary 2.6 holds.In particular, generalising (2.19), we get The decay properties of the eigenvectors (3.2) and the symmetry property with respect to γ show that this symmetric system (3.3) has all but two modes localised at the interface.These interface modes are shown in Figure 5, superimposed on one another to portray the general trend.

Infinite periodic case
For the infinite periodic case we consider a one-dimensional system constituted of N periodically repeated disjoint subwavelength resonators.We assume that for some > 0, i = for all i ∈ Z.We use the same notation as in [2].Specifically, L is the size of the unit cell and α denotes the quasiperiodicity parameter of the Floquet-Bloch transform.

Standard Brillouin zone
Similar to the finite case, we are first interested in resonant frequencies ω α (δ) for any α in the first Brillouin zone Y * := (− π L , π L ].In the periodic case, we call α → ω α (δ) band functions.
In order to get the quasiperiodic capacitance matrix one needs to apply similar modification as in the finite case to the procedure of [2], where the reader can find the details of the derivation.Definition 4.1 (Quasiperiodic gauge capacitance matrix).For γ ∈ R * , we define the quasiperiodic gauge capacitance matrix The following results hold.
Proposition 4.2.Assume that the eigenvalues of C γ,α are simple.Then the N subwavelength band functions (α → ω α i ) 1≤i≤N satisfy, as δ → 0, where V is defined by (2.22).We select the N values of ± √ δλ i having positive real parts.
The band theory of C γ,α is very rich due to the non-Hermitian structure of the matrix.Indeed one might rapidly see that for some general 1 ≤ i, j ≤ N C γ,α i,j = C γ,α j,i .From now we will consider systems of periodically repeated dimers, that is N = 2, and i = = 1.The band functions ω α present very interesting symmetries as the following lemma shows.This behaviour is observed in Figure 6, where we see the symmetry of the real parts and antisymmetry of the imaginary parts.
Therefore, the real part of the eigenvalues is symmetric and the imaginary part is antisymmetric with respect to α → −α.
(a) Real and imaginary part of ω α i for i = 1, 2. We observe the symmetries of Lemma 4.4: the real parts are symmetric and the imaginary parts antisymmetric with respect to α → −α.The following lemma establishes the exceptional points of the system, where the eigenvalues and the eigenvectors of the quasiperiodic gauge capacitance matrix simultaneously coalesce.
Lemma 4.5.Let N = 2 and 0 < γ.Then the two eigenvalues of C γ,α coalesce if and only if α = ± π L and cosh(γ) The geometrical multiplicity of the eigenvalue is then 1.
In order to establish the geometrical multiplicity of the eigenvalue, which reads tr(C γ, π L )/2, we need to find the dimension of the kernel of C γ, π L − tr(C γ, π L )/2.After some algebraic manipulation, we obtain , so that its eigenvectors read which obviously agree for γ as above.
We denote by the unique critical γ satisfying (4.4).
In Figure 7 we consider a periodic chain of dimers having s 1 = 1 and s 2 = 2.We plot the two eigenvalues in the complex plane as α varies across the Brillouin zone, from −π/L to π/L.Lemma 4.5 tells us that, for these parameter values, the eigenvalues will coalesce when γ c = 0.73899.The crossing of this critical γ c corresponds to a fundamental behavioural change in the band functions.This observation leads to the following results.
Proposition 4.6.Consider a system of periodic dimers with spacings s 1 and s 2 and band functions ω α i for i = 1, 2. There exists complex points E i ∈ C with non-zero winding of ω α i for i = 1, 2 if and only if 0 < γ < γ c (s 1 , s 2 ).
Proof.The main observation of the proof is that for 0 < γ < γ c (s 1 , s 2 ) we have, for σ = (12) the permutation of two symbols, (ω This is easily seen from the fact that tr(C γ, π L ) 2 − 4 det(C γ, π L ) changes sign at γ = γ c (s 1 , s 2 ) as we can see from the formula in the proof of Lemma 4.5.Combining this with Lemma 4.4 we can conclude that ω α i draws a closed curve on the complex plane if and only if 0 < γ < γ c (s 1 , s 2 ).The non-triviality of the winding then follows by the fact that the real and the imaginary parts of ω α i are periodic with period 2π L .
Corollary 4.7.The vorticity of a system of periodically repeated dimers is non-zero if and only if 0 < γ < γ c (s 1 , s 2 ).In Figure 8, we plot just one of the eigenvalues in the complex plane, to fully illuminate the behaviour we observed in Figure 7.We can clearly see the fundamental change that occurs at the transition when γ crosses the value γ c (s 1 , s 2 ).Whereas when γ < γ c the eigenvalue forms a single closed loop, when γ > γ c this loop is broken open to form a Cshaped curve that connects to the other eigenvalue to form a single loop.At the transition between these two states, when γ = γ c , the two eigenvalues form two touching closed loops.In Figure 8(C), we show how the critical value γ c varies as a function of the ratio of the spacings between the two resonators.Of course, when s 1 = s 2 the system is just a single repeating resonator so this critical value vanishes.Otherwise, it is strictly positive.

Complex band theory and generalised Brillouin zone
We saw in Section 3 that the eigenmodes of a finite system of subwavelength resonators with imaginary gauge potential have exponentially decaying amplitudes.Further, we saw that the system displays signs of an "infinite" order exceptional point in the limit as its size becomes large.Our final results of this paper concern the precise nature of the limiting spectrum, as the size becomes arbitrarily large.
A crucial first question is how to understand the spectrum of the limiting operator.For periodic (Hermitian) systems, the natural tool to apply is the Floquet-Bloch transform.However, as shown in Section 4.1, this does not lead to the correct limit of the set of the eigenvalues associated with the finite structure as its size goes to infinity.This is because the eigenmodes of the infinite periodic system are not Bloch modes, in the traditional sense.Instead, we will consider the Floquet-Bloch transformation with complex quasi-periodicity, so that we are able to account for the potential growth or decay of the amplitude of the eigenvectors [10].
Recently, the truncated Floquet-Bloch transformation has been applied to recover information about the band structure from finite approximations in Hermitian systems [3,6].Here, since we are considering the quasi-periodicity to be complex, we apply a slightly different method.We consider a system of N resonators.For a resonance ω j , we compute the corresponding quasi-periodicity by where ω α is the subwavelength eigenfrequency of a system with one resonator repeated periodically.In Figure 9, we display this procedure for an array of N = 60 resonators.In Figure 9a, we show the minimising region of (4.7) and see clear local minima of the function |ω j −ω α | for values of α that are symmetric about the imaginary axis.In Figure 9b we show the resulting complex quasi-frequencies α j , plotted against their indices, and in Figure 9c, we plot the complex band functions α j → ω j .In Appendix C, we present the analogous numerical results for a system of periodically repeated dimers both in the Hermitian case and non-Hermitian case with complex material parameters.
The above numerical results suggest to introduce the following set, called the generalised Brillouin zone, Here, λ α+iβ(α) i for i = 1, 2, are defined in Lemma 4.3.The value of this set is immediately clear from the following convergence result.(ii) Consider a system of dimers, i.e. with i = = 1.Then, the same result as in (i) holds as the number of dimers goes to infinity.
Proof.One first notices that because of the boundedness of the spectrum and after a possible renaming, λ k have to converge in R as N → +∞.The key idea of the proof is to use the explicit formulations from Lemma 4.3 and Picard's little theorem.Indeed, this is enough for item (i) as C α → λ α ∈ C is entire and so in particular Y * α → λ α is surjective to R up to possibly one point.For item (ii) the application of Picard's little theorem to the argument of the square root shows that C α → λ α i ∈ C for i = 1, 2 fails to be surjective because of the branch cut, but the opposite signs in λ α 1 and λ α 2 compensate, concluding the proof.
In Figure 9c, we plot both the discrete pairs (α j , ω j ) and the continuous curve, computed over Y * .We see perfect agreement between the discrete and continuous data, to the extent that it is difficult to distinguish the continuous curve from the discrete crosses.This agreement is remarkable given that the finite structure here is still relatively small, having N = 60 resonators.Further numerical results, including the extension to the case of resonator dimers, can be found in Appendix C.

Concluding remarks
We have derived from first principles the mathematical theory of the non-Hermitian skin effect arising in subwavelength physics in one dimension.Through a gauge capacitance matrix formulation, we obtained explicit asymptotic expressions for the subwavelength eigenfrequencies and eigenmodes of the structure.This allowed us to characterise the system's  fundamental behaviours and reveal the mechanisms behind them.In particular, the exponential decay of eigenmodes (the skin effect) for a system of finitely many resonators was shown to be induced by the Fredholm index of an associated Toeplitz operator.We also showed how the system behaves as its size becomes large and developed appropriate tools to understand this spectral convergence.
A natural question to ask is how the non-Hermitian system considered here compares to other non-Hermitian systems.For example, while the non-Hermiticity here arises due to the imaginary gauge potential, it is possible to break Hermiticity by simply making the material parameters complex valued.This has been studied previously in the subwavelength setting [2,7].It has been shown that, under the assumption of PT-symmetry, these non-Hermitian systems support exceptional points [7] and that the analysis of infinite periodic structures can be restricted to the standard (real) Brillouin zone [2].As a result, the spectrum of the finite structure converges to the band structure of the infinite one.In Appendix D, we show that this is due to the fact that, away from exceptional points, these non-Hermitian systems with complex material parameters can be mapped to a Hermitian system under an appropriate transformation.That is, they are equivalent to a Hermitian one away from exceptional points.Conversely, the system considered in this work, which is non-Hermitian due to the introduction of an imaginary gauge potential, is fundamentally distinct from the Hermitian system and a generalised (or complex) Brillouin zone must be considered to understand the limiting spectrum.
The explicit theory we have developed is only possible because of the simple structure of the gauge capacitance matrix and the rich literature on (tridiagonal) Toeplitz matrices and perturbations thereof.The theory of systems with periodically repeated cells of K resonators remains incomplete as no similar result to Lemma A.6 is known for block-Toeplitz matrices.However, we have shown numerically that the phenomena generalise to systems of many resonators (in Appendices B and C).The other important generalisation of this work is to higher-dimensional systems, in which it is well known that the skin effect can be realised [34].This will be the subject of a forthcoming publication [1].
Lemma A.2. Let T be a Toeplitz operator with symbol f T ∈ L ∞ .Then where Ind(A) is the Fredholm index of A. In particular, for f T continuous we have where w(f, λ) is the winding number of f around λ.
From Lemma A.2 we can see that the spectrum of a Toeplitz operator potentially has non-zero measure in the complex plane.It is therefore impossible that the spectrum of a sequence of Toeplitz matrices converges to the one of a Toeplitz operator in general.The following result [29,Theorem 7.3] shows, however, that such a convergence result exists if we consider pseudospectra.
Lemma A.3.Let T be a Toeplitz operator with continuous symbol f T and let T N be the truncation of T N to the upper-left N × N submatrix.Then for any ε > 0 lim While Lemma A.3 describes the convergence of the pseudospectra of T N , we also know how the corresponding eigenvectors behave.For example, we have the following lemma, from [29, Section 7].
Lemma A.4.Let T be a Toeplitz operator with continuous symbol f T and let T N be the truncation of T to the upper-left N × N submatrix.Let λ ∈ C be such that w(f T , λ) < 0.Then, if f T is sufficiently smooth, the corresponding eigenvector decays exponentially.Furthermore, there exists some M > 1 which is such that λ is a M −N pseudoeigenvalue of T N with corresponding eigenvector v N satisfying N is the j-th component of v.

Tridiagonal Toeplitz matrices
As the capacitance matrices we handle are tridiagonal, we focus now on the well understood theory of tridiagonal Toeplitz matrices.We will denote N × N tridiagonal Toeplitz matrices by their diagonal entries Proof.This is now trivial using the explicit formulations from Lemmas A.5 and A.6.
Despite not having a surjective mapping, there is only one eigenvalue of the perturbed matrix that cannot be approximated via the Toeplitz matrix and this eigenpair is known explicitly (µ 1 = 0, v 1 = 1).

Appendix B. Dimer systems of subwavelength resonators with imaginary gauge potentials
In this appendix, we briefly present some numerical illustrations of the skin effect in a system of finitely many dimers.Figure 10a shows the eigenvectors in a system of subwavelength resonator dimers with imaginary gauge potential, which are clearly strongly localised at one edge.This is a skin effect analogue to the simpler case studied in Section 3, with the corresponding modes plotted in Figure 2.Then, Figure 10b is the analogue of Figure 3, which shows that the spectrum of repeated dimers presents a bipartition and remains purely real.to, reducible to Hermitian systems away from exceptional points.This implies that the vorticity ν defined by the same formula as in (4.6) is trivial, see [9].
We consider the following ordinary differential equation (ODE): where V is a diagonal matrix encoding the (complex) material parameters and C the (Hermitian) capacitance matrix.We assume that V C is diagonalisable and invertible, i.e., we are away from exceptional points.Equation

Figure 2 .
Figure 2. Eigenvector localisation for a system of N = 36 resonators with s = = 1 and γ = 0.5.Following (3.2) the eigenmodes present an exponential decay and are thus localised on the left.The top left eigenmode is associated to the trivial eigenvalue 0.

Figure 3 .
Figure3.ε-pseudospectra of C γ for ε = 10 −k for k = 1, . . ., 5. Shaded in grey is the region where the symbol f T has negative winding, for T the Toeplitz operator corresponding to the semi-infinite structure.The pseudospectra are computed using[18].

Figure 4 .
Figure 4.As N → ∞, the finite structure exhibits an exceptional point of "infinite" order: arbitrarily many eigenmodes become close to parallel.Here shown in the case s = = 1 and γ = 0.5.

Figure 5 .
Figure 5. Plot of all the eigenmodes localised at the interface associated to a system described by (3.3).The x-axis encodes the site index of the resonators.Simulation performed with a structure of N = 50 resonators, = s = 1 and γ = 1.Two trivial eigenmodes are not shown.

2 Figure 7 .
Figure 7. Eigenvalue behaviour in the complex plan as α varies over the Brillouin zone for different values of γ for a system of periodically repeated dimers with s 1 = 1 and s 2 = 2.The value γ c (1, 2) = 0.73899 corresponding to a system with exceptional point signs the change from open to close band structure.

Figure 8 .
Figure 8. Band function winding and its relation to γ c .The crossing of γ c induces a change in the winding of the band functions and thus of the vorticity ν γ .Computed for a system of N = 2 periodically repeated resonators with spacing s 1 = 1 and s 2 = 2.

Figure 9 .
Figure 9. Complex band structure of a system of N = 60 resonators with = s = 1 and γ = 1.The non-trivial imaginary part encodes the exponential decay of the eigenmodes presented in Section 3.

T. 1 ≤.for 2 ≤
N (a, b, c)The eigendecomposition of tridiagonal Toeplitz matrices is well understood, as the following result[25, Section 2]  shows.Lemma A.5 (Eigenvalues and eigenvectors of tridiagonal Toepliz matrices).The eigenvalues of T N (a, b, c) are given byλ k = b + 2 √ ac cos k N + 1 π , 1 ≤ k ≤ N, k ≤ n and 1 ≤ i ≤ N,where u(i)k is the i-th coefficient of the k-th eigenvector.The spectral theory of tridiagonal matrices that are almost Toeplitz has been object of many studies.In this work we are interested in N × N matrices with the modified Toeplitz structureT N (a, b, c)For matrices of the type T N (a, b, c) the spectral theory is fully understood as well, as the following lemma[35, Theorem 3.1]  shows.Lemma A.6.Suppose that ac = 0, a + b + c = 0 and let µ be an eigenvalue of T N (a, b, c) then either µ = µ 1 := 0 and the corresponding eigenvector isv 1 = 1 or µ = µ k := b + 2 √ ac cos π N (k − 1) , 2 ≤ k ≤ N, k ≤ N and 1 ≤ j ≤ N.It should be noted that the condition a + b + c = 0 guarantees that µ 1 = 0. Proposition A.7.Let N ∈ N and (λ i,N , u i,N ) for 1 ≤ i ≤ N and (µ i,N , v i,N ) for 2 ≤ i ≤ N be the eigenpairs of T N (a, b, c) and T N (a, b, c) respectively so that the conditions of Lemma A.6 are satisfied.Then there exists an injective map σ N : {2, . . ., N } → {1, . . ., N } such that for any i = 2, . . ., N, |µ σ N (i),N − λ i,N | → 0 as N → ∞ and v σ N (i),N − u i,N → 0 as N → ∞.