Quantized classical response from spectral winding topology

Topologically quantized response is one of the focal points of contemporary condensed matter physics. While it directly results in quantized response coefficients in quantum systems, there has been no notion of quantized response in classical systems thus far. This is because quantized response has always been connected to topology via linear response theory that assumes a quantum mechanical ground state. Yet, classical systems can carry arbitrarily amounts of energy in each mode, even while possessing the same number of measurable edge states as their topological winding. In this work, we discover the totally new paradigm of quantized classical response, which is based on the spectral winding number in the complex spectral plane, rather than the winding of eigenstates in momentum space. Such quantized response is classical insofar as it applies to phenomenological non-Hermitian setting, arises from fundamental mathematical properties of the Green’s function, and shows up in steady-state response, without invoking a conventional linear response theory. Specifically, the ratio of the change in one quantity depicting signal amplification to the variation in one imaginary flux-like parameter is found to display fascinating plateaus, with their quantized values given by the spectral winding numbers as the topological invariants.

T opological quantization has captivated a generation of physicists since the discovery of the quantum-Hall effect 1 . In the more recent years, its standing as a novel quantum phenomenon was further strengthened by the observation of discretely varying Hall conductances in the quantum anomalous Hall 2,3 and quantum spin Hall effects 4 . Rigorously established through linear response theory, the link between topological winding numbers and quantized conductivity has indeed earned a place as a classic result of quantum condensed matter physics. Indeed quantized Hall response and nontrivial Chern topology are now widely regarded to be almost synonymous.
Yet, the concept of quantized response has so far eluded classical systems. Without a quantum mechanical ground state, classical systems are not amendable to conventional linear response theory, which expresses quantized responses as perturbations upon a ground state. While classical metamaterials like photonic crystals, acoustic structures, and electrical circuits can possess an integer number of topologically protected boundary states [5][6][7][8] , their response behaviors are not based on the number of accessible channels, but on analog solutions of differential equations that are by no means quantized.
In this work, we introduce the paradigm of quantized response based on the winding of the spectrum in the complex energy plane, rather than the winding of eigenstates in momentum space. For a long time, the complex spectral winding number has been exploited in predicting the non-Hermitian pumping under the open boundary conditions (OBCs) [known as the non-Hermitian skin effect (NHSE)] 9-12 , leading to the breaking of bulk-boundary correspondence and various anomalous topological phenomena  . In a recent experiment, arbitrary spectral winding is observed by visualizing the frequency band structure of optical frequency modes 38 . However, no directly measurable quantity has been associated with the spectral winding. Notably, spectral winding as a topological feature is not directly related to quantum physics and is hence a classical concept. Furthermore, the notion of classical response can be also a property of the Green's function alone, which has been recently shown to exhibit exponential growth in the presence of nontrivial spectral winding numbers [39][40][41] . While nontrivial complex energy winding is in principle well defined for quantum systems, it is most physically relevant in classical settings like mechanical, photonic, plasmonic, and electrical systems where non-Hermiticity does not present significant measurement difficulties 30,31,42 . Indeed, electrical circuits are governed by circuit Laplacians whose complex eigenvalues merely indicate phase shifts or steady-state impedances, rather than ephemeral excitations.
Specifically, what we find quantized is the response of the logarithm of the Green's function components with respect to an imaginary flux-like local parameter that continuously adjust the system boundary conditions, forming some quantum-Hall-like plateaus quantized according to the spectral winding number. This discovery was inspired by the observation that deforming a non-Hermitian system from periodic to open boundaries (PBCs to OBCs) 12,13,24,28,32 , which is closely related to complex flux insertion, always reduces the spectral winding number one by one till it reaches zero. In a classical setup subject to a steady-state drive e.g., a circuit lattice with an input current, the quantized quantity can be the logarithmic impedance experienced by the response field, e.g., the voltage. This intriguing result is rooted in the way non-Bloch eigenstates explore the interior of spectral loops, and has been investigated in the new context of scale-free non-Hermitian pumping 43 .

Results
Motivation for quantized Green's function response. In conventional literature, quantized response usually refers to quantized linear response in a quantum setting, where an occupied quantum state is driven by a time-dependent perturbation (see Supplemental Note 1 for a brief introduction). However, in the classical settings that we shall focus on i.e., photonic waveguides, acoustic lattices, and electronic circuit, the system does not settle into a ground state, and the Kubo formula describing the linear response, which concerns the response due to occupied quantum states, is inapplicable. While an integer number of topological modes have been observed in various photonic 7,[44][45][46] and acoustic [47][48][49][50][51] settings, each classical topological mode can be excited, however, weakly or strongly depending on the energy put into it, and is thus not quantized in this sense. The classical response corresponds to the flow of arbitrary amounts of optical, phononic or electric current, and not the modified (quantized) occupancy of quantum mechanical eigenstates. Consider an external coherent drive ϵ = (ϵ 1 , ϵ 2 , . . . ϵ L ) with different amplitudes applied to each of the L sites of a lattice 40 . For a steady-state drive with a fixed frequency ω, ϵðtÞ ¼ ϵðωÞ expðÀiωtÞ, and the resultant classical response field at the same frequency can be written as ϕðtÞ ¼ ϕðωÞ expðÀiωtÞ with the response field amplitude given by analogous to the Kubo formula, which is exclusively for quantum settings. Here G is the Green's function matrix and γ represents an overall gain/loss in the system. For a signal entering the system from a single site x, ϵ only possesses one nonzero component ϵ x , and the induced field at another site y is ϕ y = G yx ϵ x . In particular, the directional signal amplification of a signal entering one end of a 1D chain and measured at the other end is described by the two matrix elements G 1N and G N1 39,40 .
In the spectral representation, the Green's function Eq. (1) takes the form where Ψ L=R n are the left/right eigenstates corresponding to the nth eigenenergy E n . Its matrix elements G xy can be computed by evaluating Ψ L=R n at x and y. Equation (2) is valid in both classical and quantum settings, since the Hamiltonian is just the operator that describes time evolution, and is well defined regardless of whether position and momentum commute.
Ordinarily, we do not expect the matrix elements of G (or functions of them) to respond to an external influence β in a quantized manner, since there is no reason why the derivatives of ðω À E n Þ À1 and the eigenstates should conspire to add up to an integer. However, when translation symmetry is broken, the eigenstates can potentially become exponentially localized likẽ e κx , such that the matrix elements of G are dominated by the largest spatial decay rate κ ¼ κ max , with G xy $ e κ max ðyÀxÞ . In such special scenarios, the response of ln G xy for a fixed interval x − y is wholly dependent on how κ max varies with the external influence. In particular, if κ max were to vary with an external parameter in a quantized manner, so will the response quantity ln G xy .
In this work, we discover that ln G xy indeed possess such a quantized response if the external influence β were to be an impurity parameter tuning the boundary conditions, which coincides with tuning an imaginary flux when the latter is sufficiently weak 12 . This quantized quantity is furthermore equal to the winding number of the energy spectrum in the complex energy plane. In the following sections, we shall elaborate on these rather surprising findings, and show the classical quantized response can be measured. While this quantization applies to both classical and quantum systems, we shall call it the quantized classical response to distinguish it from the topological Hall response that exclusively exists in quantum systems.
Point-gap topology and PBC-OBC spectral evolution. To elucidate the role of spectral winding and motivate a suitable notion of classical response, we consider a generic one-band non-Hermitian system with N lattice sites, described by the following tight-binding Hamiltonian with t j the hopping amplitude across |j| lattice sites,ĉ x the annihilation operator of a (quasi-)particle at the xth lattice site, c xþN ¼ĉ x representing the PBC, and r (l) the maximal range of the hopping toward right (left) direction. The associated momentum-space Hamiltonian is given by with k the quasi-momentum, z ≔ e ik , and P r+l (z) a (r + l)th-order polynomial of z. For any t j ≠ t Ã Àj (assuming t j = 0 if j > l or j < − r), the Hamiltonian becomes non-Hermitian and possesses a point-gap topology, characterized by a nonzero spectral winding number w.r.t. a reference energy E r enclosed by the PBC spectrum 10,11,52 , with C being the Brillouin zone (BZ), i.e., k varying from 0 to 2π. Simply put, ν(E r ) gives the number of times the PBC spectrum winds around E r , as illustrated by a representative example in Fig. 1a, corresponding to the Hamiltonian of Eq. (3) with r = l = 2. As a side note, ν(E r ) also reflects the degeneracy of eigenstates at E r when the system is placed under semi-infinite boundary conditions (SIBC) 10 . Unlike the loop-like PBC spectrum depicting a nontrivial point-gap topology, the OBC spectrum must not cover any finite area in the complex plane, and so generically must take the form of curves within the PBC spectral loops 12,13 , e.g., the Y-shape lines formed by the gray dots in Fig. 1a. That is, any reference energy E r inside a PBC spectral loop is enclosed ν(E r ) times by the PBC spectrum as k varies from 0 to 2π, but the same E r cannot be enclosed by the OBC spectrum. The important qualitative insight is hence the following: if we continuously deform the system from PBC to OBC, the evolving spectrum gradually collapses from the PBC loop spectrum to the OBC line spectrum, and is therefore expected to pass E r for ν(E r ) times.
To further appreciate this understanding, we consider the realspace Hamiltonian of Eq. (3) with the following substitution only at the system's boundary, i.e., t j → e −β t j , if x + j > N or x + j < 1. This introduces an additional scaling factor, or an impurity, to the boundary couplings of a 1D chain with N unit cells. If one now continuously varies β from 0 to infinity, a PBC-OBC spectral evolution can be examined in detail.
As shown in Fig. 1a, by considering only t ±1 and t ±2 , we already have a representative and intriguing model whose spectral winding number can be either 1 or 2 in the topologically nontrivial regime. Next we consider the tuning of the boundary coupling via t ±1 → e −β t ±1 for x = N and t ±2 → e −β t ±2 for x = N, N − 1. Let the nth right eigenvalue of this model system with such PBC-OBC interpolations be E n (β). Then the spectral evolution is all captured by E n (β) vs β. To motivate the connection with the Green's function, and to capture incidences when E n (β) comes close to E r , we also define a quantity i.e., the absolute sum of the inverse energy spacings between the evolving spectrum E n (β) and a reference energy E r . For an actual system always of finite size, E n (β) is discretized, but it can still be made to be very close to the PBC reference eigenvalue E r for a sufficiently large system. As such, the quantity I β (E r ) can be a diagnosis tool to examine how many times E n (β) visits (the proximity of) E r as β is tuned. Moreover, the OBC limit can be essentially reached once β is beyond a critical value β = β OBCÑ η with η the effective localizing length of the eigenstates 43,53 . With these understandings, one infers that as β varies from 0 to β OBC , I β (E r ) is expected to display high peaks whenever the complex spectral evolution passes through E r . As explained above, the total number of such local peaks then reflects the spectral winding number ν(E r ). In Fig. 1b we illustrate I β (E r ) as a function of β for several E r denoted by the stars of different colors in Fig. 1a, corresponding to different spectral winding numbers ν(E r ). It is indeed observed that the number of peaks of I β (E r ) directly reflects the spectral winding number ν(E r ). In the Methods section, we offer more insights based on the so-called generalized Brillouin zone (GBZ), to better understand why ν(E r ) The quantity I β (E r ) = ∑ n |1/(E n (β) − E r )|, where E r is the reference energy marked in (a) with the same color (yellow, red, and blue), and E n (β) is the nth eigenenergy of the system with the hoppings across the boundary multiplied by e −β . I β (E r ) diverges exactly ν times when E r sits in a region of spectral winding ν. Parameters are set at t 1 = 1, t −1 = 0.5, t 2 = 2, and t −2 = 0, with N = 100.
can be captured by the number of singularities encountered throughout the complex spectral evolution.
Quantized response in signal amplification. While the previously defined quantity I β (E r ) is useful in diagnosing the spectral winding, it is not directly measurable. Below, we show how it inspires another analogous quantity that is directly associated with signal amplification. In particular, the quantity introduced below displays quantized plateaus that precisely match spectral winding numbers, making it possible to distinguish between one nontrivial point-gap topology from another. This is a true advance as compared with earlier interesting attempts where signal amplification was only used to probe NHSE under OBC 39,40 .
For our system with the PBC-OBC interpolation parameter β (0 < β < β OBC ), what enters into the expression of the Green's function G in Eq. (2) is 1/[E r − E n (β)] again (as in I β (E r )) and hence the Green's function should be able to capture the complex spectral evolution. More importantly, it is found that the associated eigenstates under PBC-OBC interpolation pile up at its boundary in a manner qualitatively different from that of NHSE 43 . With further derivation in the Method section, we find the amplification ratio from site N to site 1 for a one-dimensional system with only the nearest-neighbor couplings, provided t 1 > t −1 and E r falls within the loop spectrum of the system at a given β (see Method). For t −1 > t 1 , one can analogously obtain d ln jG N1 j dβ ¼ 1, corresponding to a winding number ν(E r ) = − 1. We emphasize that the variation in Eq. (7) can be directly measured in a steady-state response experiment. In an electrical circuit setting, this task can be done via impedance measurements, as will be elaborated later.
For systems with solely mth-nearest-neighbor coupling, our analysis detailed in Methods section shows that the system can be viewed as m independent subchains and the overall amplification can be captured by the m × m diagonal block at the top-right (bottom-left) corners of the overall Green's function G, denoted as G ←,m×m (G →,m×m ), corresponding to measuring the output at the first (last) m sites of a signal entering from the last (first) m sites. Since each subchain yields an amplification factor proportional to e β , one directly obtains that jG ;m m j det½G ;m m / e mβ or jG !;m m j det½G !;m m / e mβ , depending on the amplification direction. Specifically, or analogously ν !;m d ln jG !;m m j=dβ is quantized at m. Clearly then, ν ←,m (ν →,m ) counts the number of independent amplified modes whose amplification factor has the e β dependence. On the other hand, the PBC spectrum of the system winds m times around the origin of the complex plane, a fact obviously true since the associated momentum-space Hamiltonian is dominated by the terms e ±imk . For more general cases with coexisting couplings across r (l) lattice sites to the right (left), the system can still be effectively understood as m ¼ Max½r; l different subchains, with t ±m viewed as the nearest-neighbor coupling on each subchain, and the rest understood as intersubchain or intrachain longer-range couplings. For example, lattice sites j, j + m, j + 2m, ⋯ , with each pair of neighbors coupled by t ±m , form the jth subchain, for j = 0, 1, 2, ⋯ , m − 1.
Although the subchains are now coupled, the above-defined det½G ;m m (det½G !;m m ) remains physically relevant: given by the product of all its eigenvalues, it represents the product of the amplification factors of all the independent eigenmodes. Therefore, ν ←,m (ν ←,m ) still counts the number of effectively independent amplified modes with the amplification factor e β . As such, there is a fascinating correspondence between quantized response as captured by ν ←,m or ν →,m and the spectral winding number. Figure 2 presents computational results for ln G ;2 2 and its derivative with respect to β as functions of β, denoted ν ←,2 , again for the same system as that in Fig. 1. (A more complicated example with hopping range up to 3 is demonstrated in Supplemental Note 2.) We also compare these results with I(E r ) defined previously for a chosen reference energy point E r = ω + iγ, with ν(E r ) = 2. It is observed that ν ←,2 as a measurable physical response shows three clear plateaus quantized at ν ←,2 = 2, 1, 0. Echoing with the jumps between these plateaus during the spectral evolution, I β (E r ) shows local peaks whenever the spectral evolution passes through E r . As shown in Fig. 2c, these transitions during the spectral evolution as a result of increasing β match precisely with the β values for which the complex spectrum touches the reference energy point and hence the spectral winding number is about to jump. That is, the transitions between these quantized plateaus have a clear topological origin and are hence identified as topological transitions. As a side remark, observing all the plateaus in our example here happens to require a broad regime of β and hence cases with very weak boundary coupling. This is, however, not a concern because the first plateau at small β is the one to reflect the topology under PBC. Being a topologically quantized response, these plateaus are also robust against disorder, as shown in Supplemental Note 3.
The results presented in Fig. 2 are particularly stimulating. Indeed, therein neither the next-nearest-neighbor coupling nor the nearest-neighbor coupling is dominating. Yet quantized plateaus at m = 1 and m = 2 are still obtained. Returning to our decoupled subchain picture above, this indicates that for different reference energy points, the behavior of this system is topologically equivalent to that of one single chain or that of two weakly coupled subchains. With this perspective, we propose to examine ν ←,m vs. different choices of m in order to fully map out the phase boundaries, without any prior knowledge of spectral winding. We thus proceed to find the maximal ν ←,m by scanning m, for sufficiently small β. The obtained value is then expected to yield ν(E r ) of the studied system under PBC. To this end we define ν ¼ Max½ν ;1 ; ν ;2 , given that our model system at most has effectively two subchains (one can similarly define ν → ). We present in Fig. 3a our results of ν ← at β = 0 for different ω and γ, with the reference energy E r = ω + iγ. The phase boundaries identified there are in excellent agreement with the actual PBC spectrum shown in Fig. 1a. The only subtlety is that we also obtain a negative value ν ← = −1 in the topological trivial regime. In fact, the quantity ν ← does not contain the topological information in this regime, as ν ←,0 is not defined in our formalism. In retrospect, the central quantity ν ←,m with arbitrarily chosen m as a positive integer is expected to yield a nonpositive value in this regime, because of the absence of a directional amplification for E r with zero spectral winding (more details in the Method section). That is, in these topologically trivial regimes both ν ← and ν → are found to be negative, the amplification factor is far less than unity, and hence there is actually no signal amplification after all. One may just exploit this additional feature to locate regimes with zero spectral winding. By contrast, a topologically nontrivial regime with negative winding ν(E r ) < 0 corresponds to the directional amplification toward the opposite direction, and can be measured by ν →,m .
Our qualitative analysis and quantitative results are so far based on single-band systems. To establish a more general and intriguing connection between spectral winding topology and signal amplification, one must verify whether quantized response still emerges in multiband systems. The answer is affirmative. This is a remarkable finding because in multiband systems, the spectral winding topology can be jointly induced by different bands on the complex plane connected all together, leaving each individual band contribute some fractional winding only. To treat multiband systems, we consider the analogous m × m Green's function block by identifying one fixed sublattice of each unit cell, a natural route because the steady-state profile on the chosen sublattice (sublattice chosen to have nonzero support of the steady state) from one unit cell to another does contain information of amplification. As an example for demonstration, we have investigated an extended Su-Schrieffer-Heeger 54 model with a high-spectral winding number in a variety of parameter regimes, with plateaus corresponding to different spectral winding numbers. [See Supplemental Note 4].
Measurement of quantized response in electrical circuits. We next elaborate on how the quantized response can be directly extracted by measuring the impedance in an electrical circuit setting. Instead of external perturbations, a circuit is most The topological phase diagram and simulation of impedance measurement in a circuit realization. a The topological phase diagram, in excellent agreement with that shown in Fig. 1, is mapped out from examining physical response functions ν ¼ Max½ν ;1 ; ν ;2 for the model depicted in the main text, with parameters set at t 1 = 1, t −1 = 0.5, t 2 = 2, t −2 = 0, and N = 100. The reference energy for defining ν ← is given by E r = ω + iγ. At the phase boundary, E r is an eigenvalue of the system under the periodic boundary condition, leading to the divergence of the Green's function matrix. b Quantization of d ln jZ m j dβ , which is directly obtainable from simulated circuit measurement data of the impedance between the first m and last m nodes of the circuit lattice, here computed with 50 nodes. Z m (Eq. (10)) is the m × m matrix of impedances between the measured nodes, and β is the parameter controlling the strengths of the couplings connecting the two ends. The complex admittance parameter Ω indicated on panel b is connected with parameters ω and γ plotted in panel a with ω ¼ ReðΩÞ and γ ¼ ImðΩÞ. The first plateau encountered when β is increased from zero gives the nonzero topological winding ν (green, yellow for ν = 1, 2). When ν = 0, the logarithmic gradients of both Z 1,2 are close to or smaller than 0 (blue), corresponding to the absence of signal amplification toward either direction. naturally driven by a steady-state AC or DC current, with voltage response given by Kirchhoff's law I = LV, where L is the circuit Laplacian matrix and the components of I and V are, respectively, the input currents and voltages at each node. Suppose that the circuit is then grounded by identical circuit components with complex admittances −Ω. In this case, the full (grounded) Laplacian becomes J ¼ L À Ω I, and the voltage distribution due to the input current are given by 55,56 where the eigenvalues E n and L/R eigenstates Ψ L=R n are that of the Laplacian L. Notably, the quantity in the square parentheses agree exactly with our definition of the Green's function (Eq. (2)), with Ω taking the role of ω + iγ, analogous to a classical version of a tunable "Fermi energy" for probing the physics around a reference energy. We extract the physical response through the impedance Z ij between the ith and jth nodes, which is related to the Green's function via Z ij = G ij + G ji − G ii − G jj . By varying the identical grounding admittances Ω via a combination of RLC components with ±π/2 relative phase shifts, we will be able to effectively access the response from different regions of the complex spectral plane.
Quantized classical response can be obtained from a circuit whose Laplacian exhibits nontrivial spectral winding. This requires effectively asymmetric couplings between nodes, which has already been demonstrated in existing topolectrical experiments through combinations of capacitors, inductors and INICs (negative impedance converter with current inversion) comprising operation amplifiers 30,57,58 , as further elaborated in the Supplemental Note 5. The boundary couplings can also be adjusted to effect the variation with β through tunable inductors connected in series with the asymmetric couplings 43 . Arbitrarily large spectral winding numbers can always be achieved by coupling sufficiently distant nodes, which can be much more feasibly done in electrical circuits compared to other platforms.
In analogy to ν ←,m from Eq. (8), we can define Z m = det Z ij | i≤m,N−j<m , the determinant of the m × m matrix of impedances between the first m nodes at one end with the last m at the other end of the circuit chain. Although it is not exactly equivalent to ν ←,m , it is expected to vary with β in a similar manner, since the impedance Z ij is dominated by the component of the Green's function that produces the directional amplification. Therefore, by keeping track of the effective β(ω) and Ω(ω), the logarithmic gradient of the impedance determinant d ln jZ m j dβ will be expected to exhibit quantized jumps as Ω(ω) crosses the boundary between regions of different topological winding ν. For m ⩽ 2, which includes the model we had considered (Figs. 1, 2, and 3), we explicitly have whose gradients are dominated by those of terms like G 1,N−1 G 2,N and G 1,N G 2,N−1 in the presence of directional amplification (toward the first lattice site). As demonstrated via the simulated measurements in Fig. 3b, the gradient d ln jZ m¼1;2 j dβ indeed exhibits plateaus quantized at the winding number ν (blue, green, yellow for ν = 0, 1, 2, respectively) where Ω is tuned to. For higher topological winding, we take the plateau that first occurs when β is increased from 0, i.e., the plateau closest to periodic boundary conditions.

Discussion
In this work, we have introduced the new paradigm of quantized classical response, where a quantized response coefficient is established from a subblock of the Green's function matrix that varies with an imaginary flux-like parameter. Being based on the topological winding properties of the Green's function in the complex energy plane, this quantization does not assume the existence of any quantum mechanical ground state, and applies to all systems, classical, and quantum. Specifically, in a variety of situations including multiband cases where the spectral winding topology is rich, we show that the spectral winding number is directly detectable as a steady-state response coefficient to changes in the boundary condition. Indeed, through the signal amplification setting, the number of independent amplifiable modes that share a common exponential dependence on the imaginary flux-like parameter can be experimentally determined, which reveals the spectral winding number. Such correspondence between spectral winding numbers and quantized response is arguably broader in scope than in the case of momentum-space topology, because spectral winding does not even require translational invariance. Our results are relevant to a number of current experimental platforms of non-Hermitian systems 30,31,33,39,[59][60][61] . In the context of classical electrical circuits, we have shown that a quantized response can be easily extracted from extremely experimentally accessible impedance measurements.

Methods
Insights based on generalized Brillouin zone. Here we use the so-called generalized Brillouin zone (GBZ) method to elaborate why ν(E r ) can be captured by the complex spectral evolution. According to the non-Bloch band theory, the OBC spectrum can be described by the PBC one in a GBZ, using a complex deformation of the quasi-momentum k → k + iκ OBC (k) 9,12,32,62,63 . The PBC-OBC spectral evolution can then be effectively described by k → k + iκ(k) with κ(k) varying from 0 to κ OBC (k), with κ OBC (k) having the minimal magnitude to yield the OBC spectrum 12,63 . The PBC-OBC spectral evolution can hence be understood as arising from tuning κ(k) and hence deforming the BZ to the GBZ, as shown in Fig. 4. Moreover, an winding number can be defined as, analogous to that of Eq. (5), with integration in the GBZ instead of the BZ. This winding number must be zero when the OBC spectrum is reached because, again, the OBC spectrum cannot enclose any finite area 10,11 . Following the Cauchy principle, the spectral winding number is found to be ν(E r ) = N zero − N pole , where N zero and N pole are the counting of zeros and poles enclosed by the integration path (BZ or GBZ) weighted by their respective orders. The conclusion is hence as simple as follows. If we continuously tune κ(k), the PBC-OBC spectral evolution must pass through different zeros of P rþl ðzÞ z r À E r h i [colored dots in Fig. 4] for a total of ν(E r ) times, such that the spectral winding number reduces from ν(E r ) to 0 eventually when the integration path approaches the GBZ. Thus, during the complex spectral evolution, the spectrum under k → k + iκ(k) must pass the reference energy E r for a is the PBC Hamiltonian of the system and E r is the chosen reference energy for calculating the winding number. The system is chosen as the same as that in Fig. 1 in the main text, i.e., H(z) = 2z 2 + z + 1/2z. total of ν(E r ) times, constituting a rather formal argument to justify our treatment in the main text.

Directional signal amplification versus the PBC-OBC spectral evolution
Cases with only nearest-neighbor coupling. As mentioned in the main text, under an external drive ϵðtÞ ¼ ϵðωÞ expðÀiωtÞ and an overall on-site gain/loss parameter γ, the resultant response field ϕ(t) can be written as ϕðtÞ ¼ ϕðωÞ expðÀiωtÞ, with ϕðωÞ ¼ Gðω; γÞϵðωÞ; Gðω; γÞ ¼ 1 where E r = ω + iγ, G is the Green's function matrix. The amplification factors for a signal toward the left and the right are described by the matrix elements G 1N and G N1 , respectively. For a non-Hermitian H, the Green's function matrix can be expressed in the spectral representation 40,41 Gðω; γÞ ¼ 1 with Ψ R n the nth right eigenstate of H with eigenenergy E R n , and Ψ L n the corresponding left eigenstate.
To be more explicit, consider the Hatano-Nelson model under the PBC-OBC interpolation, described by the following Hamiltonian also with the assumption t 1 > t −1 without loss of generality. We shall also assume N ≫ 1 as topological properties are usually more significant in large systems. Let the nth right eigenstate be Ψ R for x ∈ [2, N − 1], and Taking the following ansatz eigensolutions: ð18Þ with C n the normalization constant, we obtain e Àβ t 1 t À1 þ e ÀM R n ðNÀ2Þ ¼ e Àβ e ÀM R n ð2NÀ2Þ þ t 1 t À1 e ÀM R n N : ð19Þ Now let us discuss the possible range of M R n depending on β, t 1 , and t −1 . We first note that if Re½M R n < 0, then the two sides of Eq. (19) exponentially explode to infinity with different rates when N ≫ 1 (hence this equality cannot hold). This tells us that Re½M R n ≥ 0. Next in the case that Re½M R n ¼ 0, i.e., je ÀM R n j ¼ 1, the eigensolution in Eq. (18) is extended, which may be satisfied only when either β = 0 (PBC) or t 1 = ±t −1 (Hermitian or anti-Hermitian). To check this, we denote e ÀM R n ðNÀ2Þ ¼ e iA and e ÀM R n N ¼ e iB , with real numbers A and B. Thus Eq. (19) can be expressed as Taking modulus square on both sides of Eq. (20), we find that Thus, we arrive at leading to the above mentioned conclusion. At last it is only possible with Re½M R n > 0, i.e. je ÀM R n j < 1, such that the eigensolutions decay from x = 1 to x = N. Hence, e ÀM R n N is vanishing and we may drop the higher-order infinitesimal in Eq. (19) to arrive at e Àβ e ln t 1 where δ t = t 1 − t −1 > 0 and N − 2 ≈ N is taken as N ≫ 1. Then the decaying exponent M R n is given by Note that the 1/N factor is crucial for the quantized response to be discussed later.
The corresponding eigenenergy is given by ¼ t 1 e ÀM R n þ e Àβ t À1 e Àβþln t 1 δt À Á e M R n : Note again that Eq. (24) does not hold if ln t 1 δ t >β because we require Re½M R n ≥ 0. Therefore, the condition must be satisfied in our consideration and this can be achieved without much difficulty. Indeed, δ t ≈ 0 and β ≈ 0 are close to the Hermitian and PBC limit, respectively, and the above exponentially decaying eigensolutions no longer hold. Nevertheless, our numerical results show a clear quantized response in a wide range of parameters, indicating an intrinsic topological property of the system. To further illustrate how the eigensolutions lead to the quantized response as discussed in the main text, we assume t 1 ≫ t −1 , i.e. δ t ≈ t 1 , so that ln ðt 1 =δ t Þ % 0 and we have This assumption is mainly for conceptual simplicity. The second term in the eigenenergy expression of Eq. (25) can be neglected as we are working in a parameter regime away from the PBC limit, i.e., β ≫ 0, so that we can use the above exponentially decaying eigensolutions. The following discussion is equally valid for β ! β β À ln ðt 1 =δ t Þ, as long as the relation in Eq. (26) is satisfied. Note that here jE R n =t 1 j is determined by the ratio β/N. That is, the eigenvalues E R n will be distributed on a circle on the complex plane, whose radius depends only on t 1 as well as the ratio β/N.
Likewise, the left eigenstates under the same assumptions satisfy H y β Ψ L n ¼ E L n Ψ L n and E L n ¼ E R n À Á Ã , and they are found to be ψ L x;n ¼ C Ã n e ÀM L n ðNÀxÞ ; t 1 e ÀM L n ¼ E L n ; From the biothorgonal condition hΨ L n jΨ R n i ¼ 1, we then obtain the normalization constant With preparations above one important matrix element of the Green's function G can then be found as follows: We may now attempt to rewrite the discrete sum in Eq. (31) in terms of a loop integral of a complex variable z, by defining k n ≔ 2nπ/N and z :¼ e Àβ=N e ik n . To that end one must implicitly assume that N under consideration is sufficiently large so as to use an integral to replace the discrete sum. With this in mind, for a given β, β/N is assumed to be vanishingly small, and hence essentially we are working in the regime of |z| → 1. Under these conditions, the sum in Eq. (31) can then be evaluated by the following integral which is found to be if z 0 ≡ E r /t 1 satisfies |z 0 | < e −β/N , i.e., the pole of the integrand falls within the the integral loop. This condition leads to meaning that the reference energy E r falls within the loop spectrum of H β . The above detailed theoretical considerations indicate that, so long as E r is enclosed by the loop spectrum of H β , we have which is just the claim in the main text regarding how to use a quantized physical response to detect the spectral winding number ν(E r ), as computationally verified in Fig. 5b. The above integral also leads to G 1N = 0 for E r outside the loop spectrum of H β , namely a signal enters from site N shall vanish at site 1 when the system is far away from the PBC limit (β → 0). Approaching the PBC limit, the system becomes translational invariant and sites 1 and N are actually two neighboring sites, G 1N then must be nonzero. Overall, |G 1N | decreases when gradually turning off the boundary hopping by increasing β, leading to a negative value of d ln jG 1N j=dβ [as shown in Fig. 3 in the main text].
Next we investigate what happens if we let β exceed β c , the critical β value for which a reference energy point E r falls exactly on the loop spectrum of H β (as shown in Fig. 5a). Let us first recall the result from Eq. (27), which indicates that the radius of the loop spectrum of H β scales with β/N. Specifically, for a given lattice size N and a given reference point E r under investigation, one immediately obtains that which is clearly proportional to N. As such, to probe the regime of β ≥ β c , β should at least linearly increase with N as well. This being the case, we can no longer approximate the discrete sum in Eq. (31) as a loop integral with N → ∞ in mind, simply because it has the factor e β in its numerator, which diverges with N because in the regime of interest β diverges with an increasing N. Based on the discussions above, in the regime of β ≥ β c , what is under investigation becomes the evaluation of the same discrete sum, but with β scaling proportionally with N, and hence both becoming sufficiently large in the event of using a loop integral to replace this discrete sum. Although we still use the same complex variable z ¼ e Àβ=N e ik n to invoke a possible loop integral, we see that |z| is now far from unity. To reflect this, we now use the alternative expression of the discrete sum in Eq. (63) and then rewriting it as the following integral, with the integrand having one pole of order N at z = 0, and another 1st-order pole at z = E r /t 1 , Assuming that |z 0 | = |E r /t 1 | > e −β/N , namely, the reference energy E r is outside the loop, so z = 0 is the only pole of the integrand enclosed by the integration path, we have a value independent of β, which is again consistent with the computational results in Fig. 5b2. Further, using the previous expression of β c from Eq. (36), we find that in this Interestingly, although this magnitude e β c =t 1 exponentially larger than e β /t 1 obtained earlier for β < β c , this amplification factor is saturated and no longer depends on β in the regime of β > β c . As a side note, one might wonder why we cannot also use the loop integral in Eq. (37) to treat the first case, namely, a fixed β in the regime of β < β c but with N approaching sufficiently large values. As said earlier, in this case we essentially perform the summation under the condition of |z| = 1. Under this condition we always have z N = 1 and hence the expression in Eq. (37) is no longer useful.
In the same fashion, we can now proceed to examine G N1 , which depicts how the signal is amplified or suppressed in the other direction. The matrix element G N1 is found to be the following, Rewriting Eq. (40) in terms of a loop integral, we have e Àβ z 2 ðE r À t 1 zÞ dz: The integrand has a 1st-order pole at z 0 = E r /t 1 , and a second-order pole at z 1 = 0. Therefore we have if |E r | > t 1 e −β/N (reference energy falls outside the loop of integral), and if |E r | < t 1 e −β/N (reference energy E r falls inside the loop of integral). Here because of the factor e −β in Eqs. (40) and Eq. (41), the replacement of the discrete sum by the loop integral is always valid by assuming a sufficiently large N, i.e., regardless of whether β is assumed to be fixed or assumed to scale linearly with N. Thus, results obtained above for both β > β c and β < β c are valid, which are indeed consistent with our numerical results. Note also that for fixed β, our numerical results for finite systems give a small but nonzero G N1 when β < β c (e.g. |G N1 | ≈ e −5 in Fig. 5c2), which vanishes when further increasing N (not shown).
Overall, we obtain that the gradient of ln jG N1 j with respect to β is again quantized, with if E r is NOT enclosed by the loop spectrum of H β , corresponding to β > β c [ Fig. 5c].
If we further increase β the system shall approach the OBC limit when β = β OBC ≈ αN with α ¼ ln ð ffiffiffiffiffiffiffiffiffiffiffiffi ffi t 1 =t À1 p Þ, where the spectrum falls on the same lines as the OBC spectrum 16,53 . Nevertheless, this limit is still not exactly like the OBCs, as the two boundaries are still weakly connected. For example, due to the (weak) boundary couplings, a flux threading cannot be gauged away and can lead to fluctuation of eigenergies, unlike in real OBC cases. Indeed, from our numerical results, we indeed see that G N1 keeps decreasing after β exceeds β OBC , and becomes a constant when β ≳ 2β OBC − β c [Fig. 5c3].
Cases with only mth-nearest-neighbor coupling. Consider now a 1D non-Hermitian chain with only the mth-nearest-neighbor couplings: Here we first assume N/m is an integer, thus the system is decoupled into m identical 1D subchains, and the eigenstates satisfy for x ∈ [2, N − m], and e Àβ t m ψ R s;n þ t Àm ψ R NÀ2mþs;n ¼ E R n ψ R NÀmþs;n ; ð47Þ t m ψ R mþs;n þ e Àβ t Àm ψ R NÀmþs;n ¼ E R n ψ R s;n ; with s = 1, 2, . . . , m labelling different subchains. As in the previous discussion for the case of m = 1, we assume t m ≫ t −m to obtain some simple analytical results. Here we replace the labels x and n with x s and n s for each subchain. We take the ansatz ψ R x s ;n s ¼ C n s e ÀM R ns ðx s À1Þ ; ð49Þ with n s only takeing values from 1 to N m = N/m, given that each subchain contains only N m lattice sites, and x s = (x − s + m)/m ranging from 1 to N m , with x being s, m + s, 2m + s, . . . , N − m + s, and more importantly, Similarly, the left eigenstates are given by Again, from the biothorgonal condition hΨ L n jΨ R n i ¼ 1, we have the normalization constants Note that in the Green's function matrix, the element G 1N shall always be zero as each subchain is decoupled from the others. In this case, the directional amplification of each subchain corresponds to the element G s(N−m+s) , with As n s takes value from 1 to N m = N/m, here we need to define k s = i2mn s π/N, so that the summation can be replaced by an integral with k s varying from 0 to 2π. Similar to the case with only nearest-neighbor coupling, we then have for each subchain, when E r is enclosed by the loop-like spectrum of each subchain. This result indicates that each subchain has its own spectral winding number ν s (E r ) = 1, but for the original 1D chain with mth-nearest-neighbor couplings, the element G 1N and G N1 are zero as sites 1 and N belong to different decoupled subchains.
On the other hand, the m-subchain picture here also indicates an effective unit-cell structure with m sublattices, even though the sublattices are physically equivalent on the lattice. Thus the directional amplification of the overall system comprised by these subchains/sublattices shall be described by the combination of that of each subchain, corresponding to the corner blocks of the overall Green function matrix, for signals moving toward the left and the right side, respectively. In the above schematic scenario where the subchains are fully decoupled from each other, only the diagonal elements of the above two matrices are nonzero, i.e., ðG ;m m Þ ab % Àδ ab e β =t 1 (see Eq. (33)) when E r is enclosed by the spectrum on the complex plane, and ðG !;m m Þ ab % δ ab e Àβ t 1 =E 2 r (see Eq. (42)) when E r is NOT enclosed by the spectrum on the complex plane. This being the case, we arrive at det½G ;m m / e mβ ; det½G !;m m / e Àmβ ; in the above two cases, repectively. Thus we have ν ;m :¼ d ln jG ;m m j dβ ¼ m; corresponding to the spectral winding number ν(E r ) = m for the case with only the mth-nearest-neighbor couplings. Furthermore, this conclusion shall also be valid when N/m is not an integer. In such cases, the system still possesses the m-subchain picture, only that the two ends of one subchain are connected to those of other subchains now. This conclusion is also verified by our numerical calculations. The independent subchain picture developed here offers an important method to understand more realistic situations where couplings with different hopping ranges coexist.
Cases with couplings across multiple ranges. In the main text, using the concept of counting the number of effectively independent modes whose amplification factor has the e β dependence, we have physically explained why the determinant of a subblock of the Green's function can still be used to capture the spectral winding number. Here we use another method to illuminate on the origin of classical quantized response. Let us assume that the largest hopping distance is m lattice sites. We first assume that the system's hopping is still dominated by mth-nearestneighbor hopping so that the previous m-subchain picture applies. We write the equation of the Green's function in the following form with G m×(N−m) , H r ðNÀmÞ m , and H r m m being different blocks of the G and E r − H matrices, whose sizes are indicated by their subscripts. Note that for the two blocks of H r = E r − H, the coefficient e −β is only contained in the block H r m m . On the other hand, H r ðNÀmÞ m has nonzero and β-independent elements only in the first 2m rows. Consider then the first 2m columns of G m×(N−m) in Eq. (62). Using the m-subchain picture as the starting point of consideration, the explicit expressions of G m×(N−m) are obtained as the following, This integral expression does not depend on β. With these preparation and now inspecting Eq. (62) again, it is seen that G ←,m×m should have a e β coefficient, so as to cancel the e −β coefficient in the matrix H r m m . Therefore the determinant of G ←,m×m yields a coefficient of e mβ , and can hence reflect the spectral winding number ν(E r ) = m.
Let us now gradually strengthen all other hoppings with different length scales. Before the spectral loop hits the reference energy E r , G is well defined, the above arguments are expected to hold and hence the response plateau is quantized at m, thus indicating the topological robustness of the quantized response. In other words, this current picture breaks down when the spectral loop passes the reference energy (where (E r − H) has zero eigenvalues), during which the qualitative nature of the Green's function changes drastically. In particular, after the phase transition the system has a smaller spectral winding number, say ν(E r ) = m − 1 for example, and hence the system is now topologically equivalent to a simpler case with m − 1 subchains. One should then use (m − 1)-subchain scenario to reapply the above insight self-consistently, by considering the block G ←,(m−1)×(m−1) instead to measure the topological winding number m − 1. This is indeed confirmed in our extensive numerical results for several systems with different spectral winding numbers.

Data availability
Raw numerical data from the plots presented are available from the authors upon request.

Code availability
Although not essential to the central conclusions of this work, computer codes for generating our figures are available from L.L. and C.H.L. upon reasonable request.