Critical non-Hermitian skin effect

Abstract

Critical systems represent physical boundaries between different phases of matter and have been intensely studied for their universality and rich physics. Yet, with the rise of non-Hermitian studies, fundamental concepts underpinning critical systems - like band gaps and locality - are increasingly called into question. This work uncovers a new class of criticality where eigenenergies and eigenstates of non-Hermitian lattice systems jump discontinuously across a critical point in the thermodynamic limit, unlike established critical scenarios with spectrum remaining continuous across a transition. Such critical behavior, dubbed the “critical non-Hermitian skin effect”, arises whenever subsystems with dissimilar non-reciprocal accumulations are coupled, however weakly. This indicates, as elaborated with the generalized Brillouin zone approach, that the thermodynamic and zero-coupling limits are not exchangeable, and that even a large system can be qualitatively different from its thermodynamic limit. Examples with anomalous scaling behavior are presented as manifestations of the critical non-Hermitian skin effect in finite-size systems. More spectacularly, topological in-gap modes can even be induced by changing the system size. We provide an explicit proposal for detecting the critical non-Hermitian skin effect in an RLC circuit setup, which also directly carries over to established setups in non-Hermitian optics and mechanics.

Introduction

Lying at the boundary between distinct phases, critical systems exhibit a wide range of interesting universal properties from divergent susceptibilities to anomalous scaling behavior. They have broad ramifications in conformal and statistical field theory^1,2,3,4, Schramm–Loewner evolution^5,6, entanglement entropy (EE)^{7,8,9,10,11,12,13,14}, and many other contexts. Recently, concepts crucial to criticalities—like band gaps and localization—have been challenged by studies of non-Hermitian systems^15,16 exhibiting exceptional points^{17,18,19,20,21,22,23,24,25,26,27} or the non-Hermitian skin effect (NHSE), which are characterized by enigmatic bulk-boundary correspondence (BBC) violations, robust-directed amplifications, discontinuous Berry curvature, and anomalous transport behavior^{28,29,30,31,32,33,34,35,36,37,38,39,40}.

We uncover here a class of criticality, dubbed the “critical non-Hermitian skin effect (CNHSE)”, where the eigenenergies and eigenstates in the thermodynamic limit “jump” between different skin solutions discontinuously across the critical point. This is distinct from previously known phase transitions (Hermitian and non-Hermitian) (Fig. 1), where the eigenenergy spectrum can be continuously interpolated across the two bordering phases. A CNHSE transition, by contrast, is characterized by a discontinuous jump between two different complex spectra along with two different sets of eigenstates. As elaborated below, this behavior appears generically whenever systems of dissimilar NHSE localization lengths are coupled, no matter how weakly. Importantly, at experimentally accessible finite system sizes, the jump smooths out into an interpolation between the two phases in a strongly size-dependent manner, such that the system may exhibit qualitatively different properties, i.e., real vs. complex spectrum or presence/absence of topological modes at different system sizes. Being strongly affected by minute perturbations around the critical point, such behavior may prove useful in sensing applications^41,42.

Fig. 1: Four different types of critical transitions.

Hermitian phase transitions (a) are marked by gap closures along the real line. In non-Hermitian cases (2nd to 4th rows, axis labels omitted), spectral phase transitions can take more sophisticated possibilities in the 2D complex energy plane. For instance, the spectral topology can change under line gap closures (b) or shrink continuously to a point and re-emerge in a different topological configuration (c), without the gap ever closing³⁸. The spectrum continuously passes through a gapless or point-like regime in the first three cases, as indicated by the black arrows. The critical non-Hermitian skin effect (d), however, is special in that OBC spectrum in the thermodynamic limit, denoted E_∞, jumps discontinuously from one configuration (left), to a different configuration (middle), and to another (right) as certain parameter changes from  −ϵ to 0 (critical border), and to ϵ, for an arbitrarily small ϵ, without ever interpolating between the configurations even though the parameter is continuously tuned.

Results

Hints of the critical non-Hermitian skin effect from the general Brillouin zone

In non-Hermitian systems with unbalanced gain and loss, the spectra under periodic boundary conditions (PBCs) and open-boundary conditions (OBCs) can be very different^{28,29,31,43,44,45}. Indeed, under OBC, eigenstates due to NHSE can exponentially localize at a boundary, in contrast to Bloch states under PBCs. This also explains the possible violation of the BBC, taken for granted in Hermitian settings.

The celebrated GBZ formalism aims to restore the BBC through a complex momentum deformation^{29,30,31,36,37,38}. Rigorously applicable for bounded but infinitely large systems, it has however been an open question whether the GBZ can still accurately describe finite-size systems. The GBZ of a momentum-space Hamiltonian H(z), z = e^ik can be derived from its characteristic Laurent polynomial (energy eigenequation)

$$f(z,E):= \det [H(z)-E]=0,$$

(1)

where E is the eigenenergy. While the ordinary BZ is given by the span of allowed real quasimomenta k, the GBZ is defined by the complex analytically continued momentum k → k + iκ(k), with the NHSE inverse decay length $\kappa (k)=-\mathrm{log}\,| z|$ determined by the smallest complex deformation z → e^ike^−κ(k) such that f(z, E) possesses a pair of zeros z_μ, z_ν satisfying ∣z_μ∣ = ∣z_ν∣ for the same E^29,31,38. Due to the double degeneracy of states with equal asymptotic decay rate at these E, there exist a pair of eigenstates ψ_μ, ψ_ν that can superpose to satisfy OBCs, i.e., zero net amplitude at both boundaries. As such, provided that the characteristic polynomial is not reducible, the OBC spectrum in the thermodynamic limit (denoted as E_∞) can be obtained from the PBC spectrum via E(e^ik) → E(e^ike^−κ(k)), apart from isolated topological modes. Thus it is often claimed that the BBC is “restored” in the GBZ defined by k → k + iκ(k) or, at the operator level, with the surrogate Hamiltonian H(e^ik) → H(e^ike^−κ(k))³⁸. In general, different E (energy band) solutions can admit different functional forms of κ(k), leading to band-dependent GBZs that have recently also been described with the auxiliary GBZ formalism³⁷. Since e^ike^−κ(k) is generically non-analytic, it represents effectively non-local hopping terms³⁸. As such, the GBZ description challenges the very notion of locality, which is central to critical systems, by effectively “unraveling” the real-space eigenstate accumulation through replacing local hoppings with effectively non-local ones.

Due to the robustness of the NHSE, eigenspectra predicted from the GBZ typically are approached rapidly by the exact numerically obtained OBC spectra even for small system sizes (${\mathcal{O}}(1{0}^{1})$ sites). In principle, the convergence should be exact in the thermodynamic limit, but in practical computations, floating-point errors ϵ₀ are continuously amplified as they propagate across the system. We hence expect accurate numerical spectra only when $L\,<-\mathrm{log}\,({\epsilon }_{0})/\max (\kappa )$, a condition always checked to be satisfied here to ensure that physical phenomena presented below are not due to numerical errors common in computations with non-reciprocal systems. However, the numerical agreement in eigenspectra between finite-size systems and the GBZ predictions fails spectacularly near a critical point where f(z, E) changes from being reducible to irreducible. To understand the significance of this algebraic property of reducibility, consider a set of coupled irreducible subsystems described by the characteristic polynomial

$$f(z,E)={f}_{0}+\mathop{\prod }\limits_{i}{f}_{i}(z,E),$$

(2)

where f_i(z, E) is the characteristic polynomial of the i-th subsystem, and f₀ is a constant that represents the simplest possible form for the subsystem coupling. When f₀ = 0, f(z, E) completely factorizes into irreducible polynomials, as expected from a Hamiltonian H(z) that block-diagonalizes into irreducible sectors associated with the individual f_i(z, E)’s. In particular, the OBC spectrum of this completely decoupled scenario is derived from the independent κ_i(k)’s of each subsystem, each determined by z_μ, z_ν from the same subsystem.

Yet, a nonzero coupling f₀, no matter how small, can have marked physical consequences by hybridizing different sectors of f_i significantly. Indeed, such hybridization is inevitable in the thermodynamic limit, with OBC eigenstates formed from superpositions of eigenstates ψ_μ, ψ_ν from dissimilar subsystems, each corresponding non-Bloch momenta $-i\mathrm{log}\,{z}_{\mu /\nu }$. Hence the GBZs, i.e., κ(k)’s of the coupled system, which are defined in the thermodynamic limit, are thus determined by all pairs of ∣z_μ∣ = ∣z_ν∣ not necessarily from the same subsystem. Therefore, the GBZs in the coupled case, no matter how small is f₀, can differ from the decoupled GBZs at f₀ = 0. That is, the thermodynamic limit and the f₀ → 0 limit are not exchangeable. However, since an actual finite physical system cannot possibly possess very different spectrum and band structure upon an arbitrarily small variation in its system parameter, the GBZ picture becomes inapplicable when describing finite systems (small or large) in the presence of CNHSE.

Anomalous finite-size scaling from CNHSE

For illustration, we turn to a minimal example of two coupled non-Hermitian 1D Hatano–Nelson chains⁴⁶ each containing only non-reciprocal (unbalanced) nearest neighbor (NN) hoppings (Fig. 2a). Its Hamiltonian reads as

$${H}_{2-{{chain}}}(z)=\left(\begin{array}{ll}{g}_{a}(z)&{t}_{0}\hfill\\ {t}_{0}\hfill&{g}_{b}(z)\end{array}\right),$$

(3)

with ${g}_{a}(z)={t}_{a}^{+}z+{t}_{a}^{-}/z+V$ and ${g}_{b}(z)={t}_{b}^{+}z+{t}_{b}^{-}/z-V$, ${t}_{a/b}^{\pm }={t}_{1}\pm {\delta }_{a/b}$ being the forward/backward hopping of chains a and b. This model can be also realized with a reciprocal system with the NHSE in a certain parameter regime (Supplementary Note 1). When t₀ = 0, the two chains are decoupled, and the characteristic polynomial is reducible as f(z, E) = [g_a(z) − E][g_b(z) − E]. Each factor f_a/b(z, E) = g_a/b(z) − E determines the skin eigensolutions of its respective chain. However, even an infinitesimal coupling t₀ ≠ 0 generically makes f(z, E) irreducible, providing that the two chains correspond to different GBZs (see “Methods” section). Specifically, consider the simple case of ${t}_{a}^{+}={t}_{b}^{-}=1$ and ${t}_{a}^{-}={t}_{b}^{+}=0$. Without couplings (t₀ = 0), the two chains under OBC respectively yields a Jordan-block Hamiltonian matrix in real space, with the spectrum given by E = ±V. Because the eigenstates of the decoupled chains are exclusively localized at the first or the last site, their GBZs collapse⁴⁵. By contrast, for any t₀ ≠ 0, $f(z,E)={E}^{2}-E(z+{z}^{-1})+(z+V)({z}^{-1}-V)-{t}_{0}^{2}$ is irreducible (here $-{t}_{0}^{2}={f}_{0}$ from Eq. (2)), insofar as the eigenenergy roots $E=\cos k\pm \sqrt{{t}_{0}^{2}+{(V+i\sin k)}^{2}}$ are no longer Laurent polynomials in z = e^ik that can be separately interpreted as de facto subsystems with local hoppings. In fact, in higher degree polynomials, an algebraic expression for z may not even exist as implied by the Abel–Ruffini theorem. Importantly, the corresponding OBC E_∞ spectrum and the GBZ for t₀ ≠ 0 are now qualitatively different. As derived in the Methods section, setting ∣z_a∣ = ∣z_b∣ gives OBC spectrum (in the thermodynamic limit): ${E}_{\infty }^{2}=\frac{1-{\eta }^{2}}{1+{\eta }^{2}}+{V}^{2}+{t}_{0}^{2}\pm 2\sqrt{{t}_{0}^{2}-{\eta }^{2}+{\eta }^{2}{t}_{0}^{2}}/(1+{\eta }^{2})$, with $\eta \in {\mathbb{R}}$. Clearly, even one now takes the t₀ → 0 limit, ${E}_{\infty }^{2}$ only simplifies to ${E}_{\infty }^{2}\to {V}^{2}+\frac{1\pm i\eta }{1\mp i\eta }$, which is not the above-mentioned OBC spectrum of the two decoupled chains. Likewise, the t₀ → 0 limit of the coupled GBZ, which can be shown to be the locus of $z=\pm \sqrt{{V}^{2}+{e}^{i\theta }}-V$ and $z=1/[V\pm \sqrt{{V}^{2}+{e}^{i\theta }}]$, θ ∈ [0, 2π], has nothing in common with the collapsed GBZs of the decoupled case.

Fig. 2: Critical non-Hermitian skin effect and finite-size scaling in the model of two coupled Hatano–Nelson chains.

a The two-chain model [Eq. (3)] with hopping asymmetry in chains a, b denoted by δ_a/b, and on-site energy offset  ±V. A small inter-chain t₀ can cause significant coupling when δ_a ≠ δ_b. b Open-boundary spectra (black dots) and eigenstate profiles (insets) with N = 10, 20, and 80 unit cells and the coupling parameter t₀ = 0.01, showing very different spectral behavior at different system sizes N. At small N ≈ 10, coupling effects are negligible, with the spectrum coinciding with the real-value open-boundary E_∞ spectrum (green) in the decoupled thermodynamic limit. As N increases, the spectrum gradually approaches the open-boundary E_∞ spectrum (red) for the coupled thermodynamic limit, with hybridization becoming sharper. c The generalized Brillouin zones of the two chains in the decoupled limit (yellow and green), and that of the two bands with t₀ = 0.01 (red circles, numerically obtained at N = 80). The Brillouin zone with ∣z∣ = 1 is given by the gray dash circle. Other parameters are t₁ = 0.75, δ_a = − δ_b = 0.25, and V = 0.5.

This paradoxical singular behavior of GBZs leads to anomalous scaling behavior in finite-size systems that are more relevant to experimental setups. The discontinuous critical transition illustrated above becomes a smooth crossover between the different OBC E_∞ solutions. As the size N of a coupled system is varied, its physical OBC spectrum interpolates between the decoupled and coupled OBC E_∞ solutions. As illustrated in Fig. 2b for the two-chain model Eq. (3) at small coupling t₀ = 0.01 (with t₁ = 0.75 and δ_a = −δ_b = 0.25 for well-defined skin modes), the OBC spectrum (black dots) changes markedly from N = 10 to 80 unit cells. For small N = 10, the spectrum approximates the OBC E_∞ (green) for t₀ = 0 (which lies on the real line), with the associated GBZs given by two perfect circles in the complex plane (Fig. 2c). At large N = 80, the spectrum converges toward the true OBC E_∞ (red curve) with nonzero coupling, where the associated respective GBZs of the two bands (also shown in Fig. 2c) are much different from the two circles as decoupled GBZs. Indeed, the eigenstates for N = 10 are almost entirely decoupled across the two chains, while those for N = 80 are maximally coupled/decoupled depending on whether they approach the red/green E_∞ curves. In the intermediate N = 20 case, the OBC spectrum lies far between the two E_∞’s, and cannot be characterized by their GBZs. The size-dependent behavior of the OBC spectrum is further elaborated through a spectral-flow study³¹ in Supplementary Note 2.

Let us now explain the above-observed marked size-dependent spectra via the competition between dissimilarly accumulated skin modes and the couplings across them. The general conditions for such are unveiled in the “Methods” section. In our model (Eq. 3), the inverse decay lengths in chains a, b are given by ${\kappa }_{a/b}=\frac{1}{2}\mathrm{log}\,({t}_{a/b}^{+}/{t}_{a/b}^{-})$, which will be dissimilar as long as δ_a ≠ δ_b. After performing a similarity transformation that rescales each site j by a factor of ${e}^{j{\kappa }_{b}}$, chain b becomes reciprocal with ${\kappa }_{b}^{\prime}=0$ while chain a has a rescaled inverse decay length ${\kappa }_{a}^{\prime}={\kappa }_{a}-{\kappa }_{b}$. If ${\kappa }_{a}^{\prime}\,\ne\,0$, chain a always possesses exponentially growing skin modes scaling like $\,{e}^{{\kappa }_{a}^{\prime}N}$ at one end. As such, the coupling t₀, even if being extremely small, still affects the spectrum and eigenstates markedly as the system size N increases, as further elaborated in the “Methods” section.

Scale-free exponential wavefunctions

A hallmark of conventional critical systems is scale-free power-law behavior, particularly in the wavefunctions. Interestingly, such scale-free behavior can also be found in the exponentially decaying wavefunctions, i.e., skin modes. Shown in Fig. 3a are the profiles of the slowest decaying eigenstates ψ(x) of H_2-chain at different system sizes N = 20, 40, 60, and 80, with the horizontal axis normalized by N. These featured eigenstates belong to the top of the central black ring in Fig. 2b, with their distance from the coupled OBC E_∞ ring (red) decreasing as  ~N⁻¹. Unlike usual exponentially decaying wavefunctions with fixed spatial decay length, here ∣ψ(x)∣ ~ e^−κx with κ ~ N⁻¹ (Fig. 3b), such that the overall profile ψ(x) has no fixed length scale. Such unique scale-free eigenmodes result from the slow critical migration of the eigenstates between E_∞ solutions (Fig. 2a, inset).

Fig. 3: Anomalous scaling behaviors of the critical non-Hermitian skin effect.

a Scale-free open-boundary skin eigenstate of the largest ${\rm{Im}}[E]$ eigenenergy of H_2-chain at system sizes N = 20, 40, 60, and 80 (red, purple, blue, green, respectively). Its rescaled profile, despite decaying exponentially rather than power-law, remains invariant across different N. This scale invariance persists in the N > 20 regime, and is due to the N⁻¹ decay (dashed line) of the inverse skin depth (red dots), as plotted in the inset. Parameters follow Fig. 2’s, except with t₀ = 10⁻³. b Entanglement entropy S (blue) of a half-filled open-boundary H_2-chain at odd system sizes N, with real-space cut at $\lfloor \frac{N}{2}\rfloor$ (the floor function of N/2) and parameters t₁ = 0.58, V = 1, t₀ = 0.4, and δ_a = − δ_b = 0.25. It saturates near zero in the gapless decoupled small N regime, but scales like $\sim \frac{1}{3}\mathrm{log}\,N$ (yellow) in the gapless coupled large N regime.

Anomalous correlations and entanglement entropy

The OBC spectra can be gapped for certain system sizes where EE obeys an area law scaling, and then become gapless at other sizes where the EE scaling is replaced by a logarithmic dependence on system size⁴⁷. This indicates that the CNHSE can lead to an unusual scaling behavior of the EE. Consider for instance the OBC H_2-chain (Eq. 3) with parameters chosen to gap out the OBC spectrum at small system sizes N. With all ${\rm{Re}}[E]\,<\,0$ states occupied by spinless-free Fermions, the real-space entanglement entropy S (blue curve in Fig. 3b) exhibits a crossover from the decoupled gapped regime at N ≤ 5 where it remains a constant due to the size-independent area of boundaries (two ends), to the gapless regime N > 20 where it approaches the $\frac{1}{3}\mathrm{log}\,N$ behavior of a gapless system (yellow line). In generic CNHSE scenarios with multiple competing OBC E_∞ loci, S can scale differently at different system size regimes, choices of fillings, and entanglement cuts, challenging the notion of single well-defined scaling behavior. As further shown in Supplementary Note 3, the two-Fermion correlator 〈ψ(1)ψ(x)〉 characterizing the EE also crossovers from rapid exponential decay at small N to 1/x power-law decay at large N. Remarkably, the probability of finding another Fermion nearby generally increases drastically when the system is enlarged (with filling fraction maintained).

Size-dependent topological modes

Topological modes are usually associated with bulk invariants in the thermodynamic limit, with finite-size effects having a diminishing role in the face of topological robustness. This intuition is not necessarily true in non-Hermitian systems, as hinted from ref. ³⁴, where an infinitesimal instability can cause a Z₂ topological transition in the thermodynamic limit³⁴. Remarkably, the CNHSE here can cause topological edge modes to appear only at certain system size regimes. Consider replacing the non-reciprocal intra-chain couplings of our H_2-chain model with inter-chain couplings with non-reciprocity  ±δ_ab between adjacent unit cells (Fig. 4a), as described by the following CNHSE Su–Schrieffer–Heeger (SSH) model⁴⁸:

$${H}_{{\rm{CNHSE}}-{\rm{SSH}}}(z)={h}_{y}(z){\sigma }_{y}+{h}_{z}(z){\sigma }_{z}+{h}_{0}(z){\mathbb{I}},$$

(4)

where h_y(z) = iδ_ab(z + 1/z), h_z(z) = V + δ₋(z − 1/z), and h₀(z) = t₁(z + 1/z) + δ₊(z − 1/z), with δ_± = (δ_a ± δ_b)/2. H_CNHSE-SSH is so named because interestingly, at δ₋ = δ_ab, it can be transformed via a basis rotation σ_z → σ_x into an extended SSH model with non-reciprocal inter-cell couplings given by  ±2δ₋ and a uniform non-reciprocal next-nearest neighbor hopping given by t₁ ± δ₊ (Supplementary Note 4), which is known to possess a Z-type topologically nontrivial phase.

Fig. 4: Topological transition induced by changing the system’s size.

a Sketch of the H_CNHSE-SSH model with cross inter-chain non-reciprocal couplings  ±δ_ab. b Open-boundary spectra (black dots) at N = 20, 30, and 40 unit cells and coupling δ_ab = 0.5 × 10⁻³. The main part of the spectrum in each case behaves similarly as the model in Fig. 2b, except for a pair of topological edge states blue emerging within the point gap at zero energy. The open-boundary E_∞ spectrum is given by green and red colors in the decoupled and coupled thermodynamic limits respectively. Other parameters are δ_a = −δ_b = 0.5, t₁ = 0.75, and V = 1.2. c κ solutions (red, blue, green, and yellow surfaces) of f(z, E) = 0 as a function of the complex energy, with the same parameters in b. Intersecting regions (green and red dotted lines) give the open-boundary skin solutions of the system in the thermodynamic limit. Among them, green lines correspond to the skin solutions of two decoupled chains at δ_ab = 0. The solutions of red curves emerge at a small but nonzero δ_ab, and the skin solutions of the weakly coupled system are given by the intersecting regions with the smallest ∣κ∣, i.e., the red loop in the center and green lines at the two ends with large and small ${\rm{Re}}[E]$. d Emergence of in-gap degenerate modes as a function of δ_ab/δ_a and N with δ_a = − δ_b = 0.5, t₁ = 0.75, and V = 1.2, with the plotted boundary scaling logarithmically with N.

When δ_ab = 0, the system is decoupled into two Hatano–Nelson chains, which must be topologically trivial. The OBC spectrum E_∞ in the decoupled case and the associated inverse decay length κ are shown in Fig. 4b, c (green curves), with positive/negative κ corresponding to skin modes accumulating population at opposite boundaries. Also shown in Fig. 4b, c (red curves) are E_∞ in the coupled case and the corresponding κ for the hybridized skin modes. With small N = 20 unit cells in Fig. 4b, the finite-size OBC spectrum (gray dots) qualitatively agrees with the decoupled E_∞ (green), with a real-valued gap at E = 0 along the ${\rm{Im}}[E]=0$ axis (inset). Upon the size increase to N = 30 and then to N = 40, such a gap first closes on the complex plane and then develops into a point gap with two zero-energy degenerate modes lying in its center. The Z-type topological origin of such in-gap modes is also carefully verified in Supplementary Note 5. The gap closure and then the emergence of in-gap topological modes resemble the typical behavior of a topological phase transition. Yet, here it is an intriguing size-induced effect. Further, the emergence of in-gap modes only requires exponentially weaker inter-chain coupling (i.e., smaller δ_ab/δ_a) for larger N, as shown in the “phase” diagram shown in Fig. 4d.

Proposal for circuit demonstration

The CNHSE is most simply realized when the two subsystems have equal and opposite κ values, since the system is then net reciprocal. Consider the RLC circuit as illustrated in Fig. 5. It is governed by Kirchhoff’s law I = JV, where I, V are the input currents and potentials at nodes 1A, 1B, 2A, 2B, ... and J is the circuit Laplacian given, at AC frequency ω = (LC)^−1/2, by

$$J(k)= \, (2i\omega C\sin k){\sigma }_{z}+{r}^{-1}({\mathbb{I}}-{\sigma }_{x})+2{R}^{-1}(1-\cos k){\mathbb{I}}\\ \to \, \left(\begin{array}{cc}2\omega C{e}^{ik}+\Delta (k)&-{r}^{-1}\hfill\\ -{r}^{-1}\hfill&2\omega C{e}^{-ik}+\Delta (k)\end{array}\right),$$

(5)

where $\Delta (k)={r}^{-1}+2{R}^{-1}-2({R}^{-1}+\omega C)\cos k$. The second line was obtained via a unitary basis transformation ${\sigma }_{y}\to {\tilde{\sigma }}_{y}=U{\sigma }_{y}{U}^{-1}={\sigma }_{z}$ that transforms the circuit Laplacian into a form similar to Eq. (3) (with V = 0), which is susceptible to the CNHSE. In this rotated basis, we evidently have two effective chains coupled by  −r⁻¹, each with unbalanced gain/loss couplings that give rise to equal and opposite NHSE. Note that RLC components are all reciprocal and cannot realize the non-reciprocal effective chains individually. However, with the basis transformation given above, the two effective chains become entangled in a way such that they are net reciprocal and hence easy to realize with RLC components, as illustrated in Fig. 5.

Fig. 5: A circuit lattice realizing the critical non-Hermitian skin effect.

The Laplacian of it is given by Eq. (5) with N unit cells with open-boundary conditions. The nodes A and B of unit cells 1, 2, 3, …  are labeled accordingly. Grounded R resistors are added to make the boundary nodes experience the same outdegree.

One can experimentally demonstrate the CNHSE by building copies of the circuit with different numbers of unit cells N (or alternatively by adjusting its length with appropriately placed switches), and mapping their Laplacian (admittance) spectra via established approaches^49,50,51,52. For instance, one can systematically connect a current source I to each node α, one node at a time (the current exits through the ground), and measure the resultant electrical potentials V_β,α at each node β. The spectrum of J is given by the inverse of the eigenvalues of the matrix V_β,α/I. In the presence of the CNHSE, the spectral plots for different N should qualitatively resemble that in Fig. 2b, since Eq. (5) is of the form of Eq. (3). Due to the robustness of the skin effect, component uncertainties in an actual experiment should minimally affect the resultant spectrum, as verified by simulation results presented in Supplementary Note 6. In particular, the circuit Laplacian spectra, which manifest the CNHSE, are almost undisturbed by uncertainty tolerances of up to 20%.

Discussion

In mathematical terms, the CNHSE arises when the energy eigenequation exhibits an algebraic singularity that leads to inequivalent auxiliary GBZs across the transition. The CNHSE heralds a whole class of discontinuous critical phase transitions with rich anomalous scaling behavior, challenging traditional associations of criticality with scale-free behavior. Even a vanishingly small coupling between dissimilar skin modes can be consequential as the system size increases. This insight is much relevant to sensing and switching applications. Beyond our two-chain models, there are other scenarios that can engineer coupling between subsystems of dissimilar NHSE length scales and hence yield CNHSE (e.g., see “Methods” section for a discussion of general two-band models). In particular, we anticipate fruitful investigations in various experimentally feasible settings such as electric circuits^53,54,55,56, cold atom systems^57,58, photonic quantum walks⁵⁹, and metamaterials^41,60, all of which are investigated with finite-size systems and hence highly relevant to the CNHSE.

Methods

Discontinuous transition of GBZ in two-chain models

The discontinuous transition induced by an infinitesimal transverse coupling in the thermodynamic limit, and also the crossover in a finite system, exist only when the two decoupled chains have different κ of their OBC skin solutions. To see this, we consider a general two-chain model described by Hamiltonian

$$h(z)=\left(\begin{array}{cc}{g}_{a}(z)+{V}_{a}&{t}_{0}\hfill\\ {t}_{0}\hfill&{g}_{b}(z)+{V}_{b}\\ \end{array}\right),$$

(6)

where g_a,b(z) only contains terms with nonzero order of z. When decoupled, the two chains correspond to the polynomials g_a,b(z) + V_a,b, respectively, and possess the same κ solutions when and only when g_b(z) = cg_a(z), with c a nonzero coefficient. When a nonzero transverse coupling t₀ is introduced, the characteristic polynomial of the two-chain system can always be written in the form of

$$\begin{array}{rcl}{P}_{c}(z)&=&({g}_{a}(z)+{V}_{a}-E)({g}_{b}(z)+{V}_{b}-E)-{t}_{0}^{2}\\ &=&(c{g}_{a}(z)-A)({g}_{a}(z)-B),\hfill\end{array}$$

(7)

where A, B are two coefficients determined by other parameters. Therefore for two chains with the same κ solutions, a transverse coupling t₀ only modifies the energy offset between them, without inducing a transition of skin solutions.

Nevertheless, the above factorization does not hold when the coupling term t₀ is z-dependent, corresponding to inter-chain couplings between different unit cells. Under this condition, P_c(z) cannot be factorized into two sub-polynomials of g_a(z) and g_b(z) = cg_a(z), meaning that the skin solution is changed for the system.

GBZ solutions E_∞ for the two-chain model

For analytic tractability, we consider the case of Eq. 3 of the main text with ${t}_{a}^{+}={t}_{b}^{-}=1$ and ${t}_{a}^{-}={t}_{b}^{+}=0$ (i.e., t₁ = δ_a = −δ_b = 0.5), but nonzero b and V. We obtain

$${H}_{2-{\rm{chain}}}(z)=\left(\begin{array}{cc}z+V&{t}_{0}\hfill\\ {t}_{0}\hfill&1/z-V\end{array}\right),$$

(8)

with the characteristic polynomial given by

$$\begin{array}{rcl}f(z,E)&=&{E}^{2}-E({z}^{-1}+z)+[(z+V)({z}^{-1}-V)-{t}_{0}^{2}]\\ &=&\hskip-1.8pc\frac{V-E}{z}-z(V+E)+[{E}^{2}-{V}^{2}-{t}_{0}^{2}+1].\end{array}$$

(9)

To find the GBZ solutions E_∞ for comparison with the actual OBC solutions, we solve for roots ∣z₊∣ = ∣z₋∣ of f(z, E) = 0 (with $\Sigma ={E}^{2}-{V}^{2}-{t}_{0}^{2}+1$):

$$\begin{array}{rcl}{z}_{\pm }&=&\frac{\left(\Sigma \pm \sqrt{{\Sigma }^{2}+4({V}^{2}-{E}^{2})}\right)}{2(V+E)}\\ &=&\frac{\Sigma \pm \sqrt{{(\Sigma -2)}^{2}-4{t}_{0}^{2}}}{2(V+E)}.\end{array}$$

(10)

For ∣z₊∣ = ∣z₋∣ to hold, the square root quantity must differ from Σ by a complex argument of π/2³⁸ i.e.,

$$\sqrt{{(\Sigma -2)}^{2}-4{t}_{0}^{2}}=i\eta \Sigma,$$

(11)

where $\eta \in {\mathbb{R}}$. Simplifying, we obtain $\Sigma =\frac{2}{1+{\eta }^{2}}\left(1\pm \sqrt{{t}_{0}^{2}+{\eta }^{2}({t}_{0}^{2}-1)}\right)$ or, in terms of ${E}^{2}\to {E}_{\infty }^{2}$,

$${E}_{\infty }^{2}=\frac{1-{\eta }^{2}\pm 2\sqrt{{t}_{0}^{2}-{\eta }^{2}+{\eta }^{2}{t}_{0}^{2}}}{1+{\eta }^{2}}+{V}^{2}+{t}_{0}^{2}.$$

(12)

as in the main text, with η tracing out a one-parameter continuous spectrum. The GBZ can be numerically obtained by substituting Eq. (12) into the expression for z_± in Eq. (10) with E = E_∞. From that, we obtain two momentum values ${k}_{\pm }={\rm{Re}}[-i\mathrm{log}\,{z}_{\pm }]$ with $\kappa ({k}_{+})=\kappa ({k}_{-})=-\mathrm{log}\,| {z}_{+}| =-\mathrm{log}\,| {z}_{-}|$ as the inverse length scales. Note however that because of the proximity to the t₀ = 0 critical point, this value of κ(k_±) is significantly different from the actual inverse OBC skin depth for a large range of finite system sizes.

Dissimilar skin modes in general two-band models

In a more general picture, the CNHSE and the size-dependent variation may exist when different parts of the system have dissimilar skin accumulation of eigenmodes. In the two-chain model, we mainly consider regimes with small inter-chain couplings, thus the two energy bands (overlapped or connected in most cases) with dissimilar skin modes are mostly given by one of the two chains respectively. To unveil the condition of having dissimilar skin modes in a general two-band system, we consider an arbitrary two-band system described by a non-Bloch Hamiltonian $H(z)={h}_{0}(z){\mathbb{I}}+{\sum }_{n = 1,2,3}{h}_{n}(z){\sigma }_{n}$, with z = e^ike^−κ(k), and κ(k) a complex deformation of momentum k describing the NHSE. Its characteristic polynomial is given by

$$f(z,E)={[E-{h}_{0}(z)]}^{2}-P(z)=0,$$

(13)

with $P(z)={\sum }_{n = 1,2,3}{h}_{n}^{2}(z)$. NHSE can be described by a GBZ where the solutions of f(z, E) = 0 satisfy E_α(z_μ) = E_α(z_ν) with ∣z_μ∣ = ∣z_ν∣ and α = ±  the band index, and $\kappa (k)=-\mathrm{log}\,| z|$ gives the inverse decay length. Conventionally, NHSE is studied mostly for a system with only nonzero h₀(z) (i.e., a one-band model) or P(z) (e.g., the non-reciprocal SSH model), where the zeros of f(z, E) lead to E_± = h₀(z) and ${E}_{\pm }^{2}=P(z)$, respectively. In either case, we can see that the two bands of E_± must have the same inverse skin localization depth κ(k), as E_α(z_μ) = E_α(z_ν) must be satisfied for α = ±  with the same z_μ,ν. To have dissimilar skin modes for the two bands, h₀(z) and P(z) must both be non-vanishing, and possess different skin solutions. That is, although h₀(z_μ) = h₀(z_ν) and $P({z}_{\mu ^{\prime} })={h}_{0}({z}_{\nu ^{\prime} })$ can still be satisfied with ∣z_μ∣ = ∣z_ν∣ and $| {z}_{\mu ^{\prime} }| =| {z}_{\nu ^{\prime} }|$, we cannot have ${z}_{\mu }=z^{\prime}$ and ${z}_{\nu }=z^{\prime}$ at the same time, otherwise the same κ(k) can be obtained for the two bands.

Competition between skin localization and inter-chain coupling

As mentioned in the main text, if two coupling chains have inverse NHSE decay lengths (non-Hermitian localization length scales) κ_a, κ_b, a change of basis will bring their coupling to be effective between a chain with no NHSE, and another with an effective skin depth κ_a − κ_b. Since that entails exponentially growing skin modes scaling like ${e}^{({\kappa }_{a}-{\kappa }_{b})N}$ at one end, we expect the effect of even an infinitesimally small inter-chain coupling t₀ to scale exponentially with N, and eventually change the OBC spectrum substantially.

Consider increasing the inter-chain coupling t₀ in our two-chain model (Eq. 3 of main text) from zero. At sufficiently small t₀, we have two practically independent OBC Hatano–Nelson chains with real spectra. Their infinitesimal coupling only shifts their eigenenergies slightly along the real line. But at a critical t₀ = t_c, the OBC spectrum is rendered complex as one or more pairs of eigenenergies coalesce and repel along in the imaginary direction. Shown in Fig. 6a is the inverse exponential scaling of the critical t₀ = t_c with N. We observe that ${t}_{\mathrm{c}}^{2}{e}^{({\kappa }_{a}-{\kappa }_{b})N} \sim {\mathcal{O}}(1)$, in agreement with the intuitive expectation that t_c should scale inverse exponentially with N because the effect of t₀ scales exponentially with N. Yet, the fact that ${t}_{\mathrm{c}}^{2} \sim {e}^{-({\kappa }_{a}-{\kappa }_{b})N}$ signifies that the CNHSE is fundamentally a non-perturbative effect since it differs from ${t}_{\mathrm{c}} \sim {e}^{-({\kappa }_{a}-{\kappa }_{b})N}$ as expected from first-order perturbation theory with left and right eigenstates that are oppositely exponentially localized spatially.

Fig. 6: Inverse exponential scaling of the critical bare coupling t₀ = t_c required for the open-boundary spectrum of H_2-chain to make the transition from real to complex.

a, b t_c Versus the system’s size N and effective skin depth κ_a − κ_b, respectively. The numerical data (blue) fits very well with the predicted scaling law ${t}_{\mathrm{c}} \sim {e}^{-({\kappa }_{a}-{\kappa }_{b})N/2}$ (dashed lines) with ${\kappa }_{a}-{\kappa }_{b}=\mathrm{log}\,2$ in a and N = 40 in b. Unless specified in the figure, the parameters are t₁ = 0.75, δ_a = − δ_b = 0.25 as in Fig. 2 of the main text. In b, κ_a − κ_b is obtained from Eq. (14) with δ_a = −δ_b varying from 0.1 to 0.4.

The scaling behavior of ${e}^{({\kappa }_{a}-{\kappa }_{b})N}$ also suggests that increasing N has similar consequences as increasing the non-reciprocity in the system, the strength of which is reflected by the absolute value of (κ_a − κ_b). Therefore it is also expected that the CNHSE shall emerge when we enhance the non-reciprocity but fix N. In Fig. 6b, we show the inverse exponential scaling of the critical t₀ = t_c with κ_a − κ_b, where the inverse NHSE decay lengths are given by

$${e}^{{\kappa }_{a,b}}=\sqrt{\frac{{t}_{1}+{\delta }_{a,b}}{{t}_{1}-{\delta }_{a,b}}}$$

(14)

for the two decoupled chains. The scaling behavior versus κ_a − κ_b further confirms that ${t}_{\mathrm{c}}^{2} \sim {e}^{-({\kappa }_{a}-{\kappa }_{b})N}$.