Amoeba Formulation of Non-Bloch Band Theory in Arbitrary Dimensions

The non-Hermitian skin effect dramatically reshapes the energy bands of non-Hermitian systems, meaning that the usual Bloch band theory is fundamentally inadequate as their characterization. The non-Bloch band theory, in which the concept of Brillouin zone is generalized, has been widely applied to investigate non-Hermitian systems in one spatial dimension. However, its generalization to higher dimensions has been challenging. Here, we develop a formulation of the non-Hermitian skin effect and non-Bloch band theory in arbitrary spatial dimensions, which is based on a natural geometrical object known as the amoeba. Our theory provides a general framework for studying non-Hermitian bands beyond one dimension. Key quantities of non-Hermitian bands, including the energy spectrum, eigenstates profiles, and the generalized Brillouin zone, can be efficiently obtained from this approach.


I. INTRODUCTION
Non-Hermitian Hamiltonians have wide applications in many branches of physics, ranging from classical wave phenomena to open quantum systems [1,2].Among various non-Hermitian systems, those that have spatially periodic structures are especially important and extensively studied.Recently, their interplay with topological physics has stimulated the fruitful investigations of non-Hermitian topological states [2,3].One of the key phenomena uncovered in this direction is the non-Hermitian skin effect (NHSE), which refers to the counterintuitive feature that the nominal "bulk eigenstates" of a non-Hermitian Hamiltonian are exponentially localized at the boundary [4][5][6][7][8][9][10][11].Unlike the topological boundary modes whose number scales linearly with the boundary area, the number of skin modes is proportional to the system volume.Among other consequences, the NHSE implies that the energy spectra of a non-Hermitian system can be drastically different under the experimentally favored open boundary condition (OBC) and the theoretically efficient periodic boundary condition (PBC).In sharp contrast to the Anderson localization, the NHSEinduced exponential localization can occur in pristine non-Hermitian systems without any disorder.
As a pronounced deviation from the usual Blochwave picture, the NHSE implies that the standard Bloch band theory is insufficient to characterize a generic non-Hermitian band.For example, the OBC energy spectrum cannot be calculated within this framework.To address this serious issue, a non-Hermitian generalization of the standard band theory, known as the non-Bloch band theory, has been introduced and applied to various non-Hermitian systems [4,[12][13][14][15][16][17][18][19].A central concept in this theory is the generalized Brillouin zone (GBZ), which is the proper surrogate for the conventional Brillouin zone (BZ) in Hermitian bands.The GBZ allows efficient computation of the continuous OBC energy spectrum of a * wangzhongemail@tsinghua.edu.cnlarge system without the need of diagonalizing a large Hamiltonian matrix in real space.Meanwhile, the shape of the GBZ directly tells the information of real-space eigenstates, e.g., the skin localization lengths.Moreover, the topological invariants on the GBZ are able to predict the number of topological boundary modes of non-Hermitian systems, whereas the usual BZ-based topological numbers fail to do so.
So far, the GBZ has been well defined only in one dimension (1D).Finding a general definition and formula for GBZ in two and higher dimensions is challenging, because the well-known 1D approach is not amenable to a straightforward generalization [4,12,20].In certain special cases, the two-dimensional (2D) GBZ has been approximately defined and calculated [21,22].However, a general approach without resorting to uncontrolled approximations has been lacking.
In this paper, we present a general formulation of non-Hermitian band theory in arbitrary dimensions.Among other results, it tells how the GBZ and related spectral properties are quantitatively determined beyond 1D.In 1D, our formulation reduces to the well-known GBZ formulation.The formulation is universal in the sense that our main results are applicable to models in arbitrary spatial dimensions, and with arbitrary degrees of freedom in a unit cell.
Our theory is based on a natural geometrical object called the amoeba by mathematicians [23][24][25][26][27]. Inspired by the amoeba and related mathematical tools, we formulate a theory which characterizes the NHSE quantitatively in arbitrary spatial dimensions.In particular, it is possible to directly calculate the energy spectrum and density of states (DOS) in the thermodynamic (i.e., large-size) limit without the troublesome finite-size errors.We show in a theorem that the energy spectrum can be obtained from the shape of the amoeba.We also demonstrate, despite the geometry-dependent NHSE, the existence of a universal spectrum (amoebic spectrum) to which the OBC spectrum under any generic geometry converges.Furthermore, the amoeba inspires a definition and the associated algorithm of the GBZ in arbitrary spatial dimensions.Among other applications, this amoeba-based GBZ provides a starting point for calculating non-Bloch band topology beyond 1D.
The remainder of this article is arranged as follows.In Sec.II, we go through the existing method of determining the GBZ in 1D, and then try to find clues for its higher-dimensional generalization.In Sec.III we introduce the basic mathematical properties of the amoeba and the associated Ronkin function, which will be useful in calculating the DOS and the GBZ.In Secs.IV-VI, we introduce the amoeba formulation for non-Hermitian systems, and then make use of this formulation and the theory of Toeplitz matrices to establish a universal way to determine the DOS as well as the GBZ.Non-Bloch band topology based on the proposed GBZ theory is studied in Sec.VII.Finally, in Sec.VIII, several useful inequalities on the OBC and PBC spectra are proved from the amoeba approach.

II. MOTIVATION A. Review of 1D non-Bloch band theory
We start with reviewing the concept of GBZ in 1D [4,12,28,29], searching for clues of its higher-dimensional generalizations.
A general 1D tight-binding Hamiltonian with OBC can be written as |i, a⟩ (t j−i ) ab ⟨j, b| , where i, j are the position indices, and a, b are the indices for intracell degrees of freedom (band indices).Explicitly, the matrix H looks like in which each t n stands for a square matrix (t n ) ab .As the hopping matrix t j−i depends only on the spatial distance j − i, we have translational symmetry in the bulk.The Hermiticity condition t n = t † −n is not required since non-Hermitian Hamiltonians are our focus.A finite hopping range n c is assumed so that the hopping matrix t n = 0 when |n| > n c .Given this real-space Hamiltonian, the corresponding Bloch Hamiltonian is the Fourier transform of t n : h(e ik ) = n t n e ik n . (3) Note that the Bloch Hamiltonian has been written as h(e ik ) instead of the more frequently used h(k).This simplifies our notations when the real-valued wave vector k is generalized to the complex plane, which amounts to making the substitution e ik → e µ+ik = β (k and µ are real-valued).Evidently, h(β) = n t n β n is a matrixvalued Laurent polynomial of β.
To solve the real-space Schrödinger equation H |ψ⟩ = E |ψ⟩, we first note that the equations have identical form at all spatial coordinates j (except some points near the boundary): Thanks to the translational invariance in the bulk, the Schrödinger equation is a linear recurrence equation with constant coefficients.The standard approach is to take a trial solution ψ j = vβ j (v is a vector whose dimension is the number of bands), and the bulk equation gives If this trial solution can appear as a component of an eigenstate wave function, it is necessary that which is called the characteristic equation for this problem.The most general form of a wave function satisfying the bulk equations is then a linear superposition of functions v n (β n ) j , where β n is a solution of det[E − h(β)] = 0 and v n is the corresponding eigenvector.Thus, an eigenstate of energy E is expressed as where the coefficients c n are to be determined by the boundary conditions.We expand the characteristic polynomial as det with a −M (E), a N (E) nonzero, and sort its roots as It has been found that, in the thermodynamic (large-size) limit, the OBC boundary condition results in the following simple equation, known as the GBZ equation [4,12,29]: It turns out that all the β M , β M +1 solutions form a closed curve called the GBZ on the complex plane of β, which contains key information about the eigenstate profiles, including the conventional wave vector k and the spatial decay rate µ of a skin mode.Once the GBZ is obtained, one can insert β ∈ GBZ into Eq.( 7) to obtain the OBC energy spectrum.Furthermore, the topological boundary modes are dictated by the topological invariants defined in the GBZ rather than in the BZ.This phenomenon is known as non-Bloch bulk-boundary correspondence [4,12].Thus, in many senses, GBZ plays a similar role as BZ does in Hermitian systems.Band theory based on the GBZ concept is known as the non-Bloch band theory.
The condition |β M (E)| = |β M +1 (E)| lies at the heart of the 1D non-Bloch band theory.It can be intuitively justified with the simplest example, in which there are only one band and two β's, β 1 and β 2 , with The open boundary condition amounts to adding fictitious sites j = 0 and j = L + 1 at the left and right ends, with ψ 0 = 0 and ψ L+1 = 0, respectively.At the left end, we find ψ 0 = 0 =⇒ c 1 = −c 2 .Then, at the right end, it is possible for the two terms to cancel each other only if term would be much larger than the β L+1 1 term.In the general case, a rigorous treatment is to write down the system of linear equations for the boundary conditions, and it turns out that the coefficient matrix is an (M + N ) × (M + N ) square matrix.Letting its determinant vanish results in [12,29].
Unfortunately, we will see that in higher dimensions, the above treatment relying on the small rank of the coefficient matrix will become quite intractable.This has been a major obstacle in the attempt to study GBZ in higher dimensions.

B. The way to the amoeba
To demonstrate the applications of non-Bloch band theory, we consider the non-Hermitian Su-Schrieffer-Heeger (SSH) model with the Bloch Hamiltonian [4,5]: or, in terms of β, where σ x,y,z are the Pauli matrices.It is known that the eigenstates exhibit NHSE under OBC, and consequently, the OBC and PBC energy spectra are drastically different [4,5].As an illustration, the OBC and PBC spectra are plotted in Fig. 1(a), for parameter values t 1 = t 2 = 1, t 3 = 0.7, and γ = 4/3.This choice of parameters is in the topologically nontrivial regime, and therefore topological edge modes with E = 0 are found in the OBC energy spectrum.The characteristic equation Eq. ( 7 Remarkably, the GBZ equation, Eq. ( 9), enables us to find the OBC energy spectra and other physical quantities without the need of diagonalizing a large realspace Hamiltonian whose size grows with the system size.The thermodynamic-limit quantities are obtained directly from the GBZ equation.A natural question arises: What is the higher-dimensional counterpart of Eq. ( 9)?
We hope to generalize the non-Bloch band theory to 2D, such that the energy spectrum can be obtained from the GBZ instead of the real space.The characteristic equation takes the same form as in 1D (and the reasoning is also the same): For our single-band model, the determinant can be dropped, and therefore the characteristic equation is simply E−h(β) = 0. Notably, the zero locus of det[E−h(β)], namely the β-solution space of det[E−h(β)] = 0, is (real) two dimensional.In fact, the solution space can be locally parameterized by the complex-valued β x or β y , one of which is determined by the other via the characteristic equation.In contrast, the zero locus of E − h(β) for a 1D system is zero dimensional, that is, several isolated β points.
The plausible next step, following the approach in the 1D non-Bloch band theory, is to select the legitimate β's by adding the boundary conditions.For 1D systems with OBC, the number of constraint equations imposed by the boundary condition does not grow with the system size L, which greatly eases the derivation of the GBZ equation, Eq. ( 9) [4,12].For 2D OBC systems (e.g., with square or disk geometry), however, the number of constraint equations is proportional to the linear size L. It is therefore challenging to exploit all these boundary-condition equations.Consequently, it is difficult to obtain a 2D counterpart of Eq. ( 9) from the approach similar to 1D.
Although a straightforward generalization to 2D looks infeasible, we can still find some clues from the 1D GBZ construction.Equation (9) and Figs.1(b) and 1(c) suggest that the moduli of the solutions to the characteristic  27].We notice that a hole exists in the amoeba in Fig. 1(e), for which the energy E 3 does not belong to the OBC spectrum.In contrast, there is no hole in the amoeba in Fig. 1(f), with energy E 4 belonging to the OBC spectrum.Viewed from this amoeba perspective, 1D non-Hermitian energy spectra also exhibit similar behaviors.In 1D, the amoeba consists of discrete points, and the hole is simply the open interval between two adjacent points.For example, the open interval (µ 2 , µ 3 ) in Fig. 1(b) can be viewed as a hole, which closes in Fig. 1(c) with µ 2 = µ 3 .
From the above examples, we observe that the absence (presence) of a hole in the amoeba of the characteristic polynomial could be an indicator of the energy E being (not being) in the OBC energy spectrum.This is a key observation of the present work.To obtain more quantitative results from this observation, it is helpful to know some mathematical properties about the amoeba.

III. MATHEMATICAL PROPERTIES OF THE AMOEBA AND RONKIN FUNCTION
In this section, we shall introduce the basic concept of amoeba and a closely related analytic tool, the Ronkin function.As a quite recent concept in mathematics, the amoeba was introduced by Gelfand et al. in 1994 [23].Albeit elementary, the notion of amoeba has deep connections with various concepts in algebraic geometry, which has stimulated extensive studies in mathematics [24][25][26][27].
Let f be a Laurent polynomial of β j , j = 1, 2, . . ., d, where d will be identified as the spatial dimension in our study.The amoeba of f is defined as the log-moduli of the zero locus of f , The name amoeba was motivated by its appearance in 2D: It has slim "tentacles" extending to infinity, and sometimes several "vacuoles" (holes) inside its body.Importantly, a particular hole plays an important role in our formulation.It is known that the amoeba in any spatial dimensions is a closed set, and each hole is a convex set [25].
A useful analytic tool in the study of amoeba is the Ronkin function, which is defined as [24,30,31] where the domain of integration is the d-dimensional torus It is beneficial to study the gradient of the Ronkin function [25,30].To this end, we can express the integrand in R f as log |f | = Re log f .The real part can be taken at the end of the calculation.It turns out that the integral is real-valued before taking the real part, and therefore the "Re" symbol can be discarded.The derivation proceeds as In the last line, we use ∂ µj = −i∂ θj acting on f .We observe that is the winding number of the phase of f along a circle parametrized by θ j , and therefore it is always real-valued.Thus, the gradient ν j is the average of the winding number w j on the (d − 1)-dimensional torus parametrized by (θ 1 , . . ., θ j−1 , θ j+1 , . . ., θ d ): For example, in 2D, one has The next fact we see is that the Ronkin function is exactly linear on each connected component of the complement of the amoeba.We call each of these (a) (c)  11) for (a, b); it is Eq. ( 12) for (c, d).The parameter values are the same as stated in Fig. 1.In both 1D and 2D, the Ronkin function is strictly linear on each component (hole) of the complement of the amoeba, where the gradient equals the integer index.The Ronkin function is always convex.Consequently, when the central hole exists, the minimum is reached on the central hole; otherwise, the minimum is reached at a single point in the amoeba.
components (either bounded or unbounded) a hole of the amoeba.In fact, when µ is not in the amoeba, f (e µ+iθ ) ̸ = 0 is satisfied in the entire T d parametrized by (θ 1 , . . ., θ d ).It follows that w j is a constant integer as {θ 1 , . . ., θ j−1 , θ j+1 , . . ., θ d } vary, and therefore the average ν j is the same integer.Thus, we can assign to each amoeba hole an integer tuple ν = (ν 1 , . . ., ν d ) dubbed the order of the hole.The orders of two different holes cannot be the same [25].Crucially, there exists at most one hole with order ν = (0, . . ., 0), which we call the central hole.Because the order is zero, the Ronkin function is a constant in this hole.Moreover, the Ronkin function is convex in the entire µ space; i.e., ) is satisfied for any two points µ 1,2 and 0 < λ < 1 [24,30].An easy corollary of the convexity is that a Ronkin function converges everywhere: If it were −∞ at one point, convexity would imply that it is −∞ everywhere.
In 1D, the Ronkin function is closely related to Jensen's formula in complex analysis [32], which reads where g is a holomorphic function with g(0) ̸ = 0, and z k (k = 1, . . ., l) are the zeros of g enclosed by the circle |β| = R. Jensen's formula can be readily obtained from Eq. ( 16) with d = 1, in which case the gradient ν = w (the index j = 1 is redundant).In fact, the left-hand side of Eq. ( 19) is exactly the Ronkin function R g (log |R|).
To calculate it, we order the zeros of g as 16) and ( 17 which is exactly Jensen's formula Eq. (19).
We now apply the explicit formula of the Ronkin function to det so that g(β) has no pole in the complex plane and g(0) = 1.Applying Eq. ( 20) to g, we have  21) we find ∂R det(E−h) /∂µ = −M + M = 0. Therefore, Eq. ( 9) means that an energy E belongs to the OBC bulk spectrum when the central hole shrinks to zero size.Thus, the Ronkin function provides crucial information about the energy spectra.Figures 3(a In 2D and higher dimensions, it is challenging to obtain a closed form for the Ronkin function.Nevertheless, regardless of the spatial dimensions, the Ronkin function is always globally convex, and is linear in the amoeba holes, i.e., connected components of the complement of the amoeba.If there is a central hole, the Ronkin function takes minimum in the entire hole, namely, the function has a flat bottom in this hole [Fig.3 (c)], otherwise the Ronkin minimum is reached at a single point in the amoeba [Fig.3 (d)].
In view of the aforementioned relation between the GBZ and Ronkin function in 1D, it is natural to ask whether there is a deep connection between the non-Hermitian energy spectra and the Ronkin function in higher dimensions.In fact, we have already seen some numerical clues for such a connection.Our observation about 2D non-Hermitian systems at the end of Sec.II can be rephrased in terms of the Ronkin function: If and only if the minimum of the Ronkin function R det(E−h) is reached at a single point µ min (E), instead of in a hole with nonzero size, will E be in the OBC bulk spectrum.
In 1D, the single point where the Ronkin function takes the minimum is is the decay factor of an OBC eigenstate.Thus, the location of the minimum of the Ronkin function precisely determines the decay factor of an OBC eigenstate.In other words, the Ronkin function tells the shape of the GBZ.We propose that this reformulation of GBZ in terms of the Ronkin function is generalizable to non-Hermitian systems in higher dimensions, which is justified in the following sections.

IV. ENERGY SPECTRA AND DENSITY OF STATES
We now introduce the amoeba formulation of non-Hermitian energy band theory in d spatial dimensions.One of our main objectives is to calculate the DOS associated with the energy spectrum, which is defined as the number of states per area on the complex energy plane, divided by the volume of the system, in the thermodynamic limit.Because the DOS has an O(L d ) volume denominator, only the bulk states are relevant in the thermodynamic limit.All contributions from edge states, bound states, etc., vanish in the thermodynamic limit.For example, the number of possible surface states grows with the size as O(L d−1 ), which contributes O(1/L) that vanishes in the thermodynamic limit.Thus, we focus on the bulk spectrum.

A. Statement of the proposal
To describe the DOS of complex energy, it is convenient to use the language of electrostatics.Let us assign electric charge −1/N to each energy eigenvalue ϵ n , where N ∼ L d is the total number of unit cells (in this convention the total charge sums up to the number of energy bands, which is independent of the size N ).In terms of Dirac's δ function, the DOS of complex energy can be written as potential Φ(E) is given by Conversely, the DOS, or the absolute value of the charge density in the electrostatics language, can be readily obtained from the Coulomb potential by taking the Laplacian on the complex energy plane: where ∆ = , and the N → ∞ limit is taken for Φ(E).
One of our main proposals is that the Coulomb potential Φ(E), in the N → ∞ limit, can be obtained from the Ronkin function, where ϕ(E) is the minimum of the Ronkin function of det(E − h) by varying µ: or ϕ(E) = R det(E−h) (µ min ), where µ min minimizes the function.Therefore, the DOS can be directly obtained from the Ronkin function, In the special cases that the spectrum is a 1D object consisting of lines or curves in the complex energy plane, it is more preferable to define the DOS as states per length rather than states per area.In fact, the states per area will diverge in these cases, while states per length is in general finite.This is the case for 1D non-Hermitian systems, and higher-dimensional systems with real spectra.In a small neighborhood of a segment of the 1D spectrum, suppose that n is a normal vector to the segment.From the electrostatic analogy and Eq. ( 26), it is easy to see that the DOS per length denoted by ρ 1D is given by where the two derivatives are taken on opposite sides of the curve segment.

B. Numerical evidence
Before diving into a thorough analytic approach, we provide numerical evidence for the Ronkin-functionbased formula, Eq. (26).We take the model Eq.(12) [Fig.2(a)] as an example.In Fig. 4, we numerically compare the DOS derived from the Ronkin function, via Eq.( 26), and that from diagonalizing the OBC Hamiltonian.The Laplacian is implemented by discretizing the complex energy plane into grid points.The DOS obtained from the Ronkin function agrees well with that from real-space Hamiltonians, though there are some differences that can be naturally attributed to the finite-size nature of the real-space calculations.First, the spectra from real-space diagonalization in Figs.4(b) and (c) are slightly narrower in the Im E direction than the one derived from the Ronkin function in Fig. 4(a).Second, the eigenvalues from real-space diagonalization seem more likely to concentrate on the real axis.These differences can be explained by the non-Bloch PT symmetry, which features a unique size dependence of non-Hermitian spec-tra in two and higher dimensions [33].In fact, our realspace Hamiltonian [see Eq. ( 12) and Fig. 2(a)] has realvalued hoppings and therefore commutes with the complex conjugation operator: [H, K] = 0. Combined with the NHSE, it means that the model has non-Bloch PT symmetry, which implies real energy spectra when the size is small.As the size L grows to infinity, the proportion of real eigenenergies diminishes to zero [33].Since the length L adopted in Fig. 4(b) and (c) is finite, the PT symmetry breaking is incomplete, leaving a nonzero proportion of real eigenenergies.
We emphasize that the Ronkin function tells the unique universal spectrum and universal DOS of the OBC system, the precise meaning of which we explain below.It is evident that the amoeba and the Ronkin function do not use information about the shape (e.g., square or disk) or boundary details (e.g., clean or locally perturbed) of the OBC system.The Ronkin function yields the same DOS regardless of geometrical details.Therefore, this approach intrinsically assumes a detailindependent (therefore, universal) spectrum with a universal DOS, which is supported by our numerical calculation.Indeed, our numerical results indicate that the DOS of an OBC system with an arbitrary shape converges in the large-size limit to the same universal DOS, at least when certain random local perturbations are added to the boundary [e.g., in Fig. 4(b)].Note that the DOS thus obtained is independent of the specific forms of boundary randomness [34].In generic cases, even the boundary randomness is not necessary to ensure the convergence to the universal DOS.For example, the DOS of disk geometry without boundary random potential already resembles the universal DOS [Fig.4(c)].In contrast, for a polygon geometry (e.g., a square), boundary randomness significantly helps the DOS to converge to the universal DOS.Without the boundary randomness, the polygon geometry can exhibit the geometry-dependent non-Hermitian skin effect, by which different polygons may have different DOS [35].Our intuitive understanding of this phenomenon is as follows.The boundary of a polygon consists of straight line segments, which are perfectly reflective in the sense that the wave vector (momentum) parallel to the line segment is conserved during wave reflection.Thus, the boundary fails to fully mix waves with different wave vectors, and therefore the Hamiltonian can be viewed as fine-tuned rather than generic.As such, the spectrum exhibits the fingerprint of specific geometry rather than the universal spectral properties of the Hamiltonian.Boundary or bulk randomness breaks the wave vector conservation and couples waves with different wave vectors, which generates a more generic energy spectrum.An analogous phenomenon is the critical non-Hermitian skin effect in 1D [36,37].In the zero-coupling limit of two coupled 1D chains, the straightforward application of the GBZ equation, Eq. ( 9), does not yield the correct OBC energy spectrum [38].The zero-coupling limit represents a fine-tuned point, and a small interchain coupling brings the spectrum to that predicted by the GBZ theory [36,37].In 2D, our numerical results suggest that the Hamiltonian should be viewed as fine-tuned for polygon shapes, and a small randomness restores the universal spectrum characterized by the universal DOS.
To summarize, there exists a geometry-independent universal spectrum that can be calculated from the amoeba and Ronkin function.By nature, it can be called the "amoebic spectrum."The DOS of an OBC system with a generic shape always approaches the universal DOS in the large-size limit.When the DOS of an OBC system with a certain (nongeneric) shape appears to deviate from the universal DOS, this deviation can be eliminated by adding a small random local perturbation.From an experimental point of view, the universal spectrum is particularly significant because disorders are often unavoidable in realistic systems.

C. Derivation
We establish the proposal Eq. ( 26) in a few steps.We begin with Szegő's limit theorem for the determinant of a large Toeplitz matrix.Before doing so, it is appropriate here to introduce the terminology of Toeplitz matrices [39,40].A matrix A determined by A ij = a j−i is called a Toeplitz matrix; i.e., the matrix element A ij depends on the difference j −i only.It is associated with a symbol, which is a complex-valued function σ(e iθ ) = n a n e inθ , θ ∈ [0, 2π].A Toeplitz matrix is often expressed in terms of its symbol as A = T [σ].By definition, the elements of a Toeplitz matrix are the Fourier components of the symbol: The language of Toeplitz matrices is very useful in addressing tight-binding Hamiltonians.For example, a single-band real-space OBC Hamiltonian H in 1D is a Toeplitz matrix, and the Bloch Hamiltonian is its symbol; conversely, we say H is generated by the Bloch Hamiltonian.In addition, one can generalize the series {a n } from numbers to square matrices.Such a matrix A is called a block Toeplitz matrix, which corresponds to a multiband Hamiltonian.One can also generalize the indices i, j, . . .from integers to integer tuples with several components; for example, in 2D we take i = (i x , i y ).The corresponding matrix A is called a multilevel Toeplitz matrix, which can be viewed as the real-space Hamiltonian of a higherdimensional lattice model.For simplicity, we call them all Toeplitz matrices hereafter.
Szegő's limit theorem was originally established for 1D Hermitian Toeplitz matrices [41], but thereafter generalized by Widom et al. to multiband [42,43] and higherdimensional [44,45] models.To state the theorem, let us consider a subspace Ω of the d-dimensional Euclidean space.Let be a symbol, which generates the Toeplitz matrix T [σ] in Ω.For example, if we take Ω to be the d-dimensional sphere with diameter L (radius R = L/2), then each block (or element, in the single-band cases) of Szegő's limit theorem reveals the asymptotic behavior of the Toeplitz determinant in the L → ∞ limit: where N is the number of lattice points in Ω, which is proportional to L d (L stands for the linear size), and T d = [0, 2π] d is the d-dimensional torus.There are two conditions for Eq. ( 29) to hold: (i) σ(e iθ ) must be invertible for any θ on the torus T d , i.e., det σ(e iθ ) ̸ = 0. (ii) The winding number of the phase of det σ(e iθ ) along any circle in T d must be zero, so that a matrix logarithm for σ(e iθ ) is well defined.A heuristic proof of Szegő's limit theorem is available in Appendix A. Notably, when σ(e iθ ) and T [σ] are Hermitian, this theorem is consistent with the fact that the OBC bulk spectrum is asymptotically the same as the PBC spectrum.In fact, the left-hand side of Eq. ( 29) yields the logarithm of the product of all OBC eigenvalues of T [σ], while the right-hand side yields the PBC counterpart, and these two quantities should be almost equal.Szegő's limit theorem tells us that a similar relation remains valid under the aforementioned conditions even though T [σ] is non-Hermitian.Now we apply this theorem to our non-Hermitian problem.Specifically, we consider σ(e ik ) = E − h(e ik ), then the generated Toeplitz matrix is is the real-space Hamiltonian corresponding to the Bloch Hamiltonian h(e ik ).The motivation to consider E − H is the simple identity in which Φ(E) is the Coulomb potential defined in Eq. (22).We observe that this expression bears resemblance to the left-hand side of Eq. ( 29), which hints useful formulas for the Coulomb potential.Before exploiting Eq. ( 29), however, we notice that its application relies on the aforementioned two conditions.For example, it would be troublesome if det[E −h(e ik )] = 0 and therefore the symbol is not invertible at certain k points.To use Eq. ( 29), we perform a similarity transformation D −1 HD to the Hamiltonian, with D x,y = δ x,y e µ•x .It adds an e µ•(y−x) factor to H x,y , the hopping from y to x, i.e., (D −1 HD) x,y = H x,y e µ•(y−x) .As such, the corresponding Bloch Hamiltonian h(e ik ) transforms into h(e µ+ik ).
In the language of Toeplitz matrices, we have D −1 HD = T [h(e µ+ik )].In fact, given h(e ik ) = n t n e ik•n , we have h(e µ+ik ) = n t n e µ•n e ik•n , from which we can read that T [h(e µ+ik )] x,y = t y−x e µ•(y−x) = H x,y e µ•(y−x) .The transformation of E − H can be written as or It follows from Eq. ( 31) that the matrix T [E − h(β)] has the same spectrum as T [E −h(e ik )] for an arbitrary value of µ.Therefore, we have det regardless of the value of µ.We can freely choose µ such that σ(β) = E − h(β) satisfies the conditions required in applying Eq. ( 29), even if the original symbol E − h(e ik ) (with µ = 0) does not.In fact, when µ locates in the central hole of the amoeba of det(E −h), we have det[E − h(β)] ̸ = 0, and the winding number of det[E−h(β)] along any k j circle is 0 (recall the mathematical preparation in Sec.III).Thus, we can apply Eq. ( 29) to σ(β) = E −h(β) and take the real part, resulting in where β = e µ+ik with µ fixed in the central hole of the amoeba of E − h(β).As has been explained, the lefthand side of Eq. ( 32) Notably, the integral on the right-hand side is exactly the Ronkin function of det(E − h).Thus, Eq. ( 32) can be written as in which µ locates in the central hole of the amoeba.Since the Ronkin function takes its minimum in the central hole, Eq. ( 33) reduces to Eq. ( 24) in the N → ∞ limit.So far, this proof of Eq. ( 24) is incomplete because the existence of the central hole of the amoeba is assumed.When the central hole does not exist, Szegő's limit theorem Eq. ( 29) cannot be applied in its original form.Here, we propose a generalization of Eq. ( 29) to remove this limitation.The generalization makes essential use of the Ronkin function.The proposed generalization is ) in which β = e µmin+ik with µ min being the minimum location of the Ronkin function; i.e., R det σ (µ) takes its minimum at µ min .
Note that the lefthand side of Eq. ( 34) can be taken as log det T [σ] = log det T [σ(e µ+ik )] with an arbitrary µ.This is because a similarity transformation of the real-space Hamiltonian does not change its determinant [see the discussion below Eq. ( 31 vanishing of gradient ∂ µj R det σ = 0 for all j's (or simply written as ∂ µ R det σ = 0).When the amoeba of det σ(β) contains a central hole, Eq. ( 34) can be established by the same approach that leads to Eq. (33).When there is no central hole, we may intuitively view the minimum location µ min as an "infinitesimal central hole," which is consistent with the fact that the Ronkin function takes the minimum in the central hole.Although we do not find a rigorous proof of Eq. ( 34) in these general cases, our numerical results support its validity (see below).In fact, taking σ(β) = E − h(β) in Eq. ( 34) and extracting the real part lead to or which becomes Eq. ( 24) in the large-size limit L → ∞.
To confirm this, we numerically calculate the Coulomb potential Φ(E) for different system sizes, and compare it with the Ronkin minimum ϕ(E).In Fig. 5, we plot the maximal difference max E |ϕ(E) − Φ(E)| as a function of the system size L. Note that the size dependence comes solely from that of Φ(E), while ϕ(E) is independent of the size.Regardless of the shape of the OBC system, the maximal difference is in good agreement with the L −1 behavior, and converges to 0 when extrapolated to the L → ∞ limit.These behaviors are exactly what Eq. ( 36) tells.
For Hermitian bands, the DOS of OBC systems determined by our formulation using the Ronkin function is consistent with the familiar fact that the OBC bulk spectrum and the PBC spectrum are asymptotically identical.This can be proved as follows.Since the eigenenergies belong to the real axis R, one can determine the DOS by Eq. ( 27), for which knowing ϕ(E) for E / ∈ R is sufficient.For an arbitrary E / ∈ R, det[E − h(e ik )] ̸ = 0 always holds for all k ∈ [0, 2π] d because of the Hermiticity of h(e ik ).Thus, µ = 0 does not belong to the amoeba of det[E − h(β)].In other words, it belongs to one of the amoeba holes.Furthermore, the phase winding number of det[E − h(e ik )] along every k j circle is zero.Therefore, the order ν j = 0 (j = 1, . . ., d) at µ = 0, meaning that µ = 0 locates in the central hole.It follows that the Ronkin minimum is reached at µ = 0, i.e., ϕ(E) = R det(E−h) (0).On the other hand, one can readily see that R det(E−h) (0) by definition is the Coulomb potential of the PBC energy spectrum.Therefore, the OBC DOS obtained from Eq. ( 27) is the same as the PBC DOS in Hermitian cases.

V. AMOEBA HOLE CLOSING AND SPECTRAL BOUNDARY
We are now able to prove a powerful theorem about the range of the OBC spectrum in the complex energy plane.It uses only the topology of the amoeba, without having to evaluate the Ronkin function.Hence, it is more efficient when one only wants to know the range of the spectrum.
We denote by Λ the set of E where the amoeba of det[E − h(β)] does not possess a central hole.We prove the following theorem: Thus, at a certain E, if the amoeba has a central hole, the DOS is zero at this E.In other words, the bulk spectrum is restricted inside Λ.
With the results of the DOS in the previous section, the proof of this theorem is now simple.When an energy E 0 is outside Λ, one must choose µ in the central hole so that Eq. ( 33) holds.Because the shape of the central hole varies continuously as E varies, the central hole should still contain this µ for E sufficiently close to E 0 .Therefore, there exists a neighborhood of E 0 denoted by V , such that for any E ∈ V , the same µ is in the central hole of the amoeba of det(E − h).Thus, for any E ∈ V , the Ronkin minimum is ϕ(E) = R det(E−h) (µ).Recalling the definition of amoeba, we see that det[E − h(e µ+ik )] is nowhere zero in V for any k ∈ T d = [0, 2π] d .Furthermore, the Ronkin function can be expressed as (38) in which {E i (β)} are the eigenvalues of h(β) (β = e µ+ik as usual).Because det 0 and therefore its integration over T d vanishes, leading to ρ(E) = ∆ϕ(E)/2π = 0.This ends our proof of the theorem Eq. (37).Note that this proof makes use only of the original version of FIG. 6. Amoebae for several energies near the band top of the model Eq. ( 12).As we decrease the energy along the real axis, the central hole of amoeba closes at E = Et ≈ 5.959 95.Parameter values are t = 1, t ′ = 0.5, and γ = 0.2, for which the energy spectrum on a disk is shown in Fig. 1(d).
When E ∈ Λ, the DOS is generally nonzero since nothing forces it to vanish.Thus, the boundary of Λ, on which the central hole closes, coincides with the boundary of the energy spectrum.In fact, if we assume the validity of the conjecture Eq. ( 34) and therefore Eq. ( 24), we have ρ(E) = ∆ϕ(E)/2π, which is generally nonzero in Λ.In other words, the support of the DOS is exactly Λ.
In Fig. 6, we illustrate the amoeba-hole closing for the model Eq.(12).Starting from an energy well above the band top (maximum of the real part of the eigenenergies), we decrease the energy along the real axis.For E larger than a certain energy, the central hole of amoeba exists, though its size shrinks as E decreases.At E = E t ≈ 5.959 95, the central hole closes.According to our proposal, this hole-closing point is identified as the band top.Similarly, the entire spectral boundary can be delineated by the hole closing.
In Fig. 7, we compare the results obtained from the amoeba formulation and numerical calculations.We numerically diagonalize the real-space Hamiltonian with increasing sizes and obtain the values of the band top, which are then extrapolated to infinite size [Fig.7  the former [Fig.7 (a)], and ∆E b ∝ L −1 for the latter [Fig.7 (c)].It turns out that the former scaling is more accurately obeyed here.Therefore, the error of extrapolation is larger for the latter.Despite larger numerical error, the numerical results are still in good agreement with the amoeba-theoretic prediction.

VI. GENERALIZED BRILLOUIN ZONE
In this section, we establish another proposal mentioned at the end of Sec.III, namely, the location of the Ronkin minimum determines the complex momenta and hence the GBZ.According to Eq. ( 37), the DOS can be nonzero only when the central hole of the amoeba of det(E − h) is absent.Thanks to the convexity of the Ronkin function, the minimum location µ min must be unique in these cases, and therefore the proposal is unambiguous.
The GBZ essentially determines the exponential behavior of the eigenstates of a non-Hermitian lattice Hamiltonian.According to our proposal, an eigenstate with eigenenergy E is expressed asymptotically in the bulk as [46] where the sum is over all k satisfying det[E − h(e µmin+ik )] = 0, and c k are certain E-dependent coefficients.Here, µ min plays the role of the imaginary part of the wave vector if we define a complex-valued wave vector Equivalently, we can use the variable In our theory, the GBZ consists of all points β subject to For a d-dimensional non-Hermitian system, the GBZ is a d-dimensional subspace of the β space whose real dimension is 2d.In fact, there are 2d + 2 real unknowns in Eq. ( 41), namely, the real and imaginary parts of (β, E).Equation ( 41) then imposes d + 2 constraints, meaning that the solution space is d dimensional.Equation ( 41) provides a general approach to calculate the GBZ for higher-dimensional non-Hermitian systems.In practice, we often parametrize the GBZ by k, treating µ min as its function.This vectorial function µ min (k) is a complete representation of the GBZ.As an application of our theory, the GBZ thus obtained for our model Eq. ( 12) is shown in Fig. 8(a).Now we provide evidence for Eqs. ( 39) and (41).First, we show that an eigenstate |ψ⟩ subject to H|ψ⟩ = E|ψ⟩ shares the same exponential behavior with the Green's function G(E) = (E − H) −1 .To see this, we decompose the OBC Hamiltonian into its eigenstates: where the sum is over all eigenstates; |nR⟩ and ⟨nL| are the right and left eigenstates, respectively, which satisfy They are orthonormalized as ⟨mL|nR⟩ = δ mn .Suppose that x and 0 are real-space locations far from the boundary.The Green's function is related to the eigenstates in the following way: where C is an infinitesimal counterclockwise contour centered at E. We emphasize that E is taken in the OBC energy spectrum, otherwise the integral in Eq. ( 44) would vanish.The right-hand side of Eq. ( 44) is a linear combination of all eigenstates with energy E. As our objective is to obtain µ, it suffices to keep both the eigenstates and the Green's functions to the exponential level.Assuming that the Green's function does not change abruptly with respect to E, we conclude that these eigenstates share the same exponential behavior as the Green's function ⟨x| (E − H) −1 |0⟩.We would also like to remark that, while the above argument assumes that the Hamiltonian is diagonalizable, the result is general enough to hold even in the presence of exceptional points [47].
Suppose that the asymptotic behavior is To validate Eqs. ( 39) and ( 41), we need to show that µ j = (µ min ) j .We numerically confirm this relation in Figs.8(c) and 8(d).Note that in principle one may also validate Eqs.(39) and ( 41) by the exponential behavior of eigenstates themselves.However, our detour of calculating the Green's function is computationally more accurate and less expensive, enabling calculation for larger L within reasonable time.

VII. NON-BLOCH BAND TOPOLOGY
It has been found recently that the topological numbers defined on the conventional Brillouin zone fail to characterize the topological edge states and bulkboundary correspondence in non-Hermitian systems.To correctly account for the topological edge modes, the non-Bloch topological invariants defined on the GBZ have been proposed.In practice, most of their applications are restricted to 1D, because a general calculable formulation of GBZ in higher dimension has been lacking.Based on the amoeba formulation, we are now able to address the non-Bloch band topology in higher dimensions.Specifically, the non-Bloch Chern number, which was previously calculated only by continuum approximation, can now be calculated in the entire GBZ as is.
To be concrete, we consider the following non-Hermitian Chern-band model [21] h e ik = (v sin The real-space Hamiltonian reads where x are 2D integer coordinates, and e j is the unit vector in the jth direction [see Fig. 2(b)].We set t = v = 1 for simplicity.This model is known to have a Chern insulator phase as well as a trivial insulator phase in the Hermitian limit (γ = 0) [48].It has Chern number C = 1 when 0 < m < 2, C = −1 when −2 < m < 0, and C = 0 otherwise.We focus on the phase boundary between C = 0 and C = 1 .When we turn on the non-Hermitian term (γ ̸ = 0), the phase boundary extends into a curve in the m-γ plane, which should be predicted by the non-Bloch band theory.In particular, the number of chiral edge modes in each phase is given by the non-Bloch Chern number evaluated on the GBZ, which is a two-dimensional subspace of the four-dimensional β space.Alternatively, we may take k = k − iµ = log β as the coordinate system, in which the GBZ is two dimensional because µ is treated as a function of k, so that the GBZ can be parametrized by k.Thus, the non-Bloch Chern number can be viewed as the conventional Chern number of h(e i k), with k ∈ GBZ.For a band labeled by α, the non-Bloch Chern number reads where ϵ xy = −ϵ yx = 1, ∂ j refers to ∂ ∂kj , and ⟨ k, α, L| and | k, α, R⟩ are the left and right eigenstates of the α band, respectively.We shall focus on the "valence band" (Re E < 0), labeled as α = −.
More intuitively, the non-Bloch Chern number can be expressed as the change of Berry phase along circular sections of the GBZ in the y direction, i.e., C = where we have write the Bloch Hamiltonian as h = h x σ x + h y σ y + h z σ z , and h = h 2 x + h 2 y + h 2 z .Here, h x,y,z can be extracted from Eq. ( 45) by the substitution k → k: h x = v sin kx + iγ, h y = v sin ky + iγ, and h z = m − t cos kx − t cos ky .
In Fig. 9, we show the energy spectrum and the non-Bloch Chern number for the Chern-band model.As an example, we fix γ = 0.4, and take m = 2.6 and m = 1.2.From Figs. 9(b) and 9(d), we read for each case that C = 0 and C = 1, respectively.In Fig. 9(c), topological edge states are indeed seen, being consistent with the nonzero non-Bloch topological invariant C = 1.
For the present model, the precise phase boundary can even be analytically determined by the amoeba formulation.In the m-γ plane, the phase boundary between the C = 0 and C = 1 phases is a curve, which can be viewed as a function m c (γ); i.e., the phase transition point is m c (γ) for a fixed γ.On the phase boundary, the energy gap about E = 0 closes, and E = 0 belongs to the energy spectrum.Thus, we can find m c by inspecting when the central hole of the amoeba at E = 0 closes.The exact formula for m c is found to be A perturbative formula for m c was obtained in Ref. [21], and recently a result to the sixth order of γ was reported in Ref. [49].Our derivation of Eq. ( 49) is as follows.We focus on the special solution of the characteristic equation that satisfies β x = β y ≡ β, as the symmetry of the model allows.The characteristic equation det In a neighborhood of the phase transition, we numerically notice that the four roots of the above quartic equation are all real, and that the second-largest and third-largest (the middle two) roots, under the logarithm, are on the boundary of the central hole of the amoeba.Hence requiring the central hole to vanish is equivalent to requiring the middle two roots to be equal.Solving this requirement leads to the final result Eq. (49).Equation ( 49) is verified up to high precision by numerically locating where the central hole closes, as listed in Table I.From the amoeba-based approach, Fig. 10 illustrates locating the critical value m c for the case γ = 0.4.TABLE I. Topological phase boundary mc as a function of γ.The decay factor µ (equal in both x and y directions) at E = 0 and "mc amoeba" is determined using the amoeba method, with numerical error < 1 × 10 −6 ."mc numerical" is adopted from Ref. [21], which is obtained by numerical diagonalization of the real-space Hamiltonian with numerical error < 3 × 10  In practice, it is convenient to search for m c (γ) by iteration: 1. Plot the amoebae at a consecutive sequence of m's.
2. Find the interval with the smallest central-hole area (closest to hole closing).
3. Terminate if the desired accuracy is reached.
4. Otherwise, on this new interval, assign a new sequence of m's and go back to step 1.
On the other hand, m c can also be determined by numerically diagonalizing the real-space Hamiltonian and finding where the band gap closes [21].We compare the amoeba-based and numerical-diagonalization results, which agree well with each other.
In the amoeba formulation, an exact expression for the decay factor µ x = µ y ≡ µ at the band-closing point E = 0 can also be obtained: This result is very close to the µ values listed in Table I, which are numerically obtained by locating the amoebahole-closing points.We remark that a first-order approximation µ = γ + o(γ) was reported in Ref. [21].

VIII. SPECTRAL INEQUALITIES
As one of the applications of the amoeba formulation, we prove a general inequality about the spectral radius.For a matrix or operator H, the spectral radius ρ(H) is defined to be the largest absolute value of the eigenvalues.Our spectral inequality states that for a non-Hermitian lattice Hamiltonian with translational symmetry in the bulk, the OBC spectral radius is less than or equal to its PBC counterpart: Note that we are interested in the thermodynamic limit in which all edge-state contribution vanishes.It is evident that Figs.1(a) and (d) satisfy the inequality.The motivation of considering this inequality is a neat inclusion relation in 1D, which states that the OBC spectrum is in the interior of the PBC spectrum [29,50].As the non-Bloch band theory (which the 1D theorem is based on) aligns with the amoeba formulation in 1D, it is natural to ask if a parallel statement can be established in higher dimensions.We begin the proof with a simple fact about the electrostatic potential.Suppose that a Coulomb potential Φ is generated by a DOS (charge) distribution ρ in the 2D complex plane, i.e., We are interested in the average potential on a circle, say, |E| = R, which reads To calculate this average, one just needs to treat all charges inside the circle as distributed uniformly on the circle, and all charges outside as they are.Because of the additivity of the potential, one needs only to prove this property for a single point charge.The contribution of a point charge q with radial coordinate r < R is Φ = q log R, while that of a point charge with r > R is Φ = q log r.This can be readily proved using Jensen's formula Eq. (19).Taking g(Re iθ ) = Re iθ − re iα in Eq. ( 19), the left-hand side of the equation is the average potential generated by a point charge at E ′ = re iα , and Eq. ( 19) becomes Now we proceed with the spectral radius inequality.From Sec.IV, we know that the Coulomb potential generated by the OBC DOS is equal to the minimum of the Ronkin function Φ OBC (E) = min µ R det(E−h) (µ).On the other hand, we can also write the Coulomb potential generated by the PBC DOS in terms of the Ronkin function at µ = 0: Hence, It follows that the average satisfies ΦOBC ≤ ΦPBC .
Taking a circle |E| = R that surrounds the whole PBC spectrum, the average potential on this circle is ΦPBC = Q log R, where Q is the total charge, which is equal to the number of bands of the specific lattice model.Since the total charge of the OBC DOS is also Q, the aforementioned electrostatic fact implies that ΦOBC ≥ Q log R, with "=" reached when the charge density vanishes outside the circle.Combining this with Eq. ( 58), we obtain the identity ΦOBC = Q log R. Consequently, the OBC DOS is zero outside the circle, which means that the OBC spectral radius is no larger than the PBC counterpart.
An alternative proof of the spectral inequality is based on Eq. (37).We denote the PBC spectral (outer) boundary as S, which consists of one or several closed curves (note that the spectrum may also have an inner boundary, which is not included in S; for example, an annulusshaped spectrum has an outer and an inner boundary).Consider an energy E outside S. Since E is not in the PBC spectrum, µ = 0 must belong to the complement of the amoeba of det(E − h).Furthermore, one can see that µ = 0 locates in the central hole.In fact, in the |E| → ∞ limit (with µ = 0 fixed), the phase of det[E −h(β)] is determined solely by E and independent of k, and therefore the phase winding number along each k j circle is zero [cf.Eqs. ( 16) and ( 17)].It follows that, for sufficiently large |E|, the order ν = 0 at µ = 0, i.e., µ = 0 locates in the central hole.For any E outside S, one may connect E to infinity by a path that does not intersect S. Since µ = 0 is always in the same amoeba hole as E varies along this path, this amoeba hole must be the central hole.Thus, for any E outside S, µ = 0 locates in the central hole of the amoeba.It follows from Eq. (37) that the OBC DOS ρ OBC (E) = 0 outside S. Therefore, the OBC spectrum is enclosed by the PBC spectral boundary.This statement is slightly stronger than the spectral inequality Eq. ( 52).Particularly, when the energy spectrum is not convex in the complex plane, this statement implies Eq. ( 52) but not vice versa.
We can further extend the above result as follows.Taking µ to be any value other than 0, we readily find that the OBC spectrum is enclosed by any TBC(µ) spectrum, where TBC(µ) refers to a PBC Hamiltonian generated by symbol h(e µ+ik ), i.e., a Hamiltonian under "twisted" periodic boundary condition.We also mention that the development of our theorem is reminiscent of a similar theorem in 1D [50], as they both rely on the winding numbers of the symbol.
A corollary is immediately implied: The spectral range of the real (or imaginary) part of the OBC bulk spectrum is within its PBC counterpart.To be specific, one has the following general inequalities: max Re(ϵ n,OBC ) ≤ max Re(ϵ n,PBC ), (59) max Im(ϵ n,OBC ) ≤ max Im(ϵ n,PBC ), (60) min Re(ϵ n,OBC ) ≥ min Re(ϵ n,PBC ), (61) min Im(ϵ n,OBC ) ≥ min Im(ϵ n,PBC ). ( In view of the significance of the energy spectrum, the spectral inequalities have various implications, some of which have already been exploited.In one dimension, Eq. ( 60) has played an important role in the directional amplification [51][52][53].In the context of open quantum systems, the spectral inequalities are closely related to the boundary sensitivity of Liouvillian gap and relaxation time [54][55][56].Our general statement and proof of the spectral inequalities lay the groundwork for their higherdimensional applications.

IX. CONCLUDING REMARKS
In this work, we have formulated an amoeba theory of the non-Hermitian skin effect and non-Bloch band theory in arbitrary spatial dimensions.It provides a theoretical framework for studying periodic non-Hermitian systems without the serious dimensional limitation.Among other applications, our theory offers a general yet efficient approach to compute the key physical quantities of non-Hermitian systems, such as the energy spectrum, density of states, and generalized Brillouin zone.
Although the initial version of non-Bloch band theory was formulated under the OBC, the concept of GBZ in 1D is also generalizable to other boundary conditions such as the domain wall cases [57].Nevertheless, it seems that the amoeba approach, as we now understand, naturally corresponds to the standard OBC systems.Thus, we have focused on the OBC case throughout the present paper.Compared to other boundary conditions, the OBC is especially important in many senses.First, the OBC is experimentally the most relevant because realistic systems often have OBC.Second, taking the OBC enables the investigation of both bulk and boundary physics, while other boundary conditions including the PBC are blind to boundary phenomena.Third, even certain measurable physical quantities far from the boundary can be naturally expressed in terms of the GBZ from rather than the conventional BZ associated with PBC [13,58], though it is in principle free to choose the boundary condition when investigating the physics deep in the bulk.In this sense, the OBC seems to be an advantageous choice even for studying certain bulk physics.
We would also like to remark that not all aspects of this work are mathematically rigorous.Although numerical evidence is supplied whenever a mathematically strict derivation is unavailable, a fully rigorous proof of all our main results is of course desirable.
In view of the ubiquity of periodic structures in both natural and synthetic systems, it is hoped that our theory can find wide applications in the abundant non-Hermitian phenomena.For example, our formulation can be naturally applied to open quantum systems in the free-particle limit, for which the energy spectrum of the non-Hermitian Liouvillian superoperator determines the dynamics and relaxation [54,55].Our theory immediately enables calculating the relevant quantities beyond 1D.Furthermore, for many-body non-Hermitian systems, our amoeba theory may still be a good starting point for including the interaction effects, which will be left for future work.
Therefore, a condition for Eq.(A12) is that T [σ] −1 ∞ , which is the largest of its singular values, is finite in the thermodynamic limit L → ∞.This amounts to asking if T [σ] itself has a zero singular value when L → ∞.Since most of its singular values form bulk bands and have a finite gap from zero, the zero singular values should be of certain topological origin.They are intimately related to the index theorem, or the bulk-boundary correspondence in free systems.If the symbol has a nontrivial topological index, then the corresponding Toeplitz matrix has either nonempty kernel or nonempty cokernel.Either case implies that T [σ]T [σ] † is not fully ranked and has a zero eigenvalue.As to 1D Toeplitz matrices, the topological index is well known to be the negative of the winding number of det σ(e iθ ) [60].To ensure that the Toeplitz matrix does not have zero singular values, we impose a condition that the symbol σ is homotopically trivial [61].This means that there exists a continuous path in the space of invertible symbols subject to Eq. (A10), connecting the symbol σ(e iθ ) to the constant symbol σ const (e iθ ) = 1.This requirement may be too strong for now, but we see in a moment that it arises again when defining the matrix logarithm.
The asymptotic inversion formula Eq. (A12) aims at finding a good surrogate for the inverse of T [σ] in the bulk, T [σ −1 ], which is calculable by integration alone.We remark that Widom also proposed an explicit expression for the subleading term in T [σ] −1 = T [σ −1 ] + O(L d−1 ) + o(L d−1 ) in one-dimensional case [43,59], and also in higher-dimensional continuous case [44].It is not used in the current paper and therefore omitted.Szegő's limit theorem also naturally inherits an expression for the subleading term O(L d−1 ).
Next, we move on to Szegő's limit theorem Eq. ( 29).Since log det A = tr log A, it is equivalent to the following trace formula: where N is the number of lattice points in Ω.In order to properly define a logarithm for the symbol, we must rely on the homotopy path imposed earlier.Let σ λ (e iθ ) be one of these paths, which is a continuous function from λ ∈ [0, 1] to the space of invertible symbols.To ensure the path independence of the integral in Eq. (A16) below, we further require that [σ ′ λ , σ −1 λ ] = 0.For example, a possible construction of such a homotopy path in the context of taking σ(e iθ ) = E − h(e µ+iθ ) is as follows.We can find a continuous path E λ connecting E 0 = E and E 1 = ∞ that avoids the spectral range of h(e µ+iθ ) (with E λ ̸ = E − 1), then the path of symbol can be taken as A consequence of trivial homotopy is that the winding number of σ along any circle on T d is identically zero.
Moreover, it is an interesting question what role is played by point-gap topological invariants other than the winding numbers Eq. (17).Our current conjecture is that the (1D) winding numbers are all that are related to the appearance of NHSE.Mathematically, it means that the Szegő's limit theorem holds even if other topological numbers are nonzero.For example, the 3D winding number can protect topological surface states in non-Hermitian band systems, but it does not induce NHSE [62].Possibly the subleading terms in Szegő's limit theorem will be helpful in quantifying such surface states.
The logarithm in Szegő's limit theorem is now defined as log σ = To ensure that the new basis spans the entire lattice, we require all elements of U d×d to be integers, and det U = ±1.Thus, U is an element drawn from GL(d, Z).
Let the real-space Hamiltonian be H = x,y |x⟩ t y−x ⟨y| , (B2) where we suppress the intracell indices for notational simplicity, but they can be restored without effort.The Bloch Hamiltonian (also known as the symbol) then reads We may write it alternatively as a function of the complex momentum kj = −i log β j : h(e i k1 , . . ., e i kd ) = (B4) After a U transformation, the new Bloch Hamiltonian becomes h ′ (e i k′ 1 , . . ., e i k′ d ) = n t n exp i k′ We immediately see that h(e i k1 , . . ., e i kd ) = h ′ (e i k′ Renaming θ → U θ and noting that the integration measure is unchanged, we find the right-hand side to be equal to R det(E−h) (µ).
FIG. 3. The Ronkin function in (a, b) 1D and (c, d) 2D, taken at energies (a) E1, (b) E2, (c) E3, and (d) E4 stated in Fig. 1.The Bloch Hamiltonian is Eq.(11) for (a, b); it is Eq.(12) for (c, d).The parameter values are the same as stated in Fig.1.In both 1D and 2D, the Ronkin function is strictly linear on each component (hole) of the complement of the amoeba, where the gradient equals the integer index.The Ronkin function is always convex.Consequently, when the central hole exists, the minimum is reached on the central hole; otherwise, the minimum is reached at a single point in the amoeba.

( 21 )
for log |β l | ≤ µ ≤ log |β l+1 |.Thus, R det(E−h) is a piecewise linear function of µ in 1D.To relate to the concept of amoeba, we note that each log |β k | is a component of the amoeba, and each interval (|β k |, |β k+1 |) is an amoeba hole.Particularly, the open interval (|β M |, |β M +1 |) is exactly the central hole on which the Ronkin function is flat, because taking a derivative on Eq. ( ) and 3(b) are two examples of the 1D Ronkin function.The Ronkin function is flat in the central hole of Fig.3 (a)corresponding to E outside the energy spectrum, and the hole shrinks to a point in Fig.3(b) corresponding to E in the energy spectrum.

FIG. 4 .
FIG. 4. DOS from the Ronkin function and diagonalization of the real-space Hamiltonian.The Hamiltonian used here is Eq.(12) [Fig.2(a)], with parameter values t = 1, t ′ = 0.5, and γ = 0.2.(a) DOS from the Ronkin function, via Eq.(26).(b) DOS from diagonalizing the real-space Hamiltonian on a square with side length L = 130.An on-site random potential distributed uniformly in [−0.5, 0.5] is added at each boundary site.(c) DOS from diagonalizing the real-space Hamiltonian on a disk with diameter L = 140 (without random potential at the boundary).The insets show the Coulomb potential, ϕ(E) in (a), and Φ(E) in (b) and (c).To facilitate comparison with (a), the DOS in (b) and (c) is also obtained from the Coulomb potential via Eq.(23), in which Φ(E) is generated by diagonalizing the real-space Hamiltonian.

2 ×10 − 2 FIG. 5 .
FIG. 5.The maximal difference maxE |ϕ(E) − Φ(E)| between the OBC-spectrum-based and Ronkin-function-based Coulomb potentials as a function of the linear size L. The adopted Hamiltonian is Eq.(12) [Fig.2(a)].The straight line is a linear fitting from data points with L ≥ 55.(a) The real space is a square of side length L, with disorders added on the boundary [the same as in Fig. 4(b)].Each data point comes from averaging over six disorder configurations.(b) The real space is a disk of diameter L, without disorder [the same as in Fig. 4(c)].

40 FIG. 7 .
FIG. 7.Band top Et and bottom E b obtained from amoeba formulation and numerical calculations.The model is Eq.(12), with t = 1, t ′ = 0.5 fixed.(a)(c), γ is fixed to 0.2.(a) Band top obtained from diagonalizing the real-space Hamiltonian on the disk with diameter L. The extrapolation to L → ∞ agrees well with the amoeba-hole-closing point marked as the orange square.(b) Band top as a function of γ, for increasing L values.The extrapolation to L → ∞ is shown as black dots, which are in excellent agreement with the amoeba-hole-closing points marked as orange squares.The PBC result is shown as the dotted line.(c)(d) The counterparts of (a)(b) for the band bottom.The linear fitting in (c)(d) is based on data from L ∈ [80, 240].Error bars indicate 95% confidence intervals.

FIG. 8 .
FIG. 8. (a)The generalized Brillouin zone of model Eq.(12).Because of the symmetry of interchanging x and y in this model, we have (µmin)x = (µmin)y ≡ µ, which is plotted as a function of the real part of the wave vector k = (kx, ky).(b) Typical profile of a bulk eigenstate that exhibits the non-Hermitian skin effect.The eigenstate is taken at E = −0.003+ 0.056i.The system is a square of length L = 130, with certain random on-site disorders on the boundary.(c)(d) Comparison between µmin and the exponential decay rate of the Green's function.The latter is defined by linearly fitting µx, µy from log ⟨x| (E − H) −1 |0⟩ ∼ µxx + µyy, in which µx ≈ µy for our model and therefore we present the average µ = (µx + µy)/2.Here, H is a real-space Hamiltonian defined on a disk with diameter L, with random on-site disorder (distributed uniformly in [−0.5, 0.5]) on the boundary.Each circle represents the result from a disorder configuration.In the fitting, we discard x's that are within distance 20 from the boundary.In (c), we fix L = 400.In (d), we vary L at fixed E = 1, 3, 5; the corresponding µmin's are shown as horizontal lines.

9 .
FIG. 9. Energy spectrum and band topology of the non-Hermitian Chern-band model Eq.(45).(a) Energy spectrum by diagonalizing the real-space Hamiltonian on a square of side length L = 60, with random on-site disorders on boundary.Parameters are taken in the topologically trivial regime: m = 2.6, γ = 0.4.(b) The Berry phase Wy(kx) along a circle in the ky direction (kx = const).The non-Bloch Chern number C = [Wy(2π) − Wy(0)]/2π = 0. Parameter values are the same as in (a).(c) Energy spectrum by diagonalization, calculated under the same setting as in (a).Parameters are taken in the topologically nontrivial regime: m = 1.2, γ = 0.4.Topological edge states are observed.(d) The Berry phase as a function of kx, with Wy(2π) − Wy(0) = 2π and therefore C = 1.Parameter values are the same as in (c).

FIG. 10 .
FIG. 10.The amoeba-based procedure of locating the phase boundary in Table I.The figures show the amoebae, zooming in on the central hole, with γ = 0.4 fixed and m varied.The central hole closes at exactly one point m = mc ≈ 2.154 066, which is where the topological phase transition occurs.

T
[σ ′ λ ] T [σ λ ] −1 dλ, (A17)where a prime denotes a derivative with respect to λ.Using Eqs.(A12) and (A1), we deduce log T [σ] we focus on the 1-norm becomes clear now, since tr A ≤ ∥A∥ 1 for any matrix or operator A. Taking the trace on both sides, Szegő's limit theorem is proved, because the trace of the matrix on the last line is by definition tr T [log σ] =

1 ,d
. . ., e i k′ d ) as long as k1 b 1 + . ..kd b d = k′ 1 b ′ 1 + . . .k′ d b ′ d .In other words, h is an invariant scalar function.It follows that h is invariant if the wave vector transforms as amoeba (coordinates being µ j = − Im kj ) undergoes the same linear transformation U when we change the basis.Its properties of physical importance, especially the existence of the central hole, are hence unchanged under a change of basis.Furthermore, for the Ronkin function the following identity holds:R det(E−h) (µ) = R det(E−h ′ ) (U µ),(B7)where µ stands for the column vector (µ 1 , . . ., µ d ) T .This is because, by definition,R det(E−h ′ ) (U µ)log det E − h ′ (e U µ+iθ ) .(B8) −4 .