Quantum-squeezing-induced algebraic non-Hermitian skin effects and ultra spectral sensitivity

Abstract

The well-established non-Bloch band theory predicts exponential localization of skin-mode eigenstates in one-dimensional (1D) non-Hermitian systems. Recent studies, however, have uncovered anomalous algebraic localization in higher dimensions. Here, we extend these ideas to Hermitian bosonic quadratic Hamiltonians incorporating quantum squeezing, offering a genuine quantum framework to explore non-Hermitian phenomena without external reservoirs. We construct a two-dimensional (2D) bosonic lattice model with two-mode squeezing and study the spectral properties of its bosonic excitation within the Bogoliubov–de Gennes (BdG) formalism. We demonstrate an algebraic non-Hermitian skin effect (NHSE), characterized by quasi-long-range power-law localization of complex eigenstates. The system shows ultra spectral sensitivity to perturbations involving either two infinitesimal on-site impurities or a weak long-range hopping impurity, while remaining insensitive to a single weak on-site impurity. Analytical treatment via the Green’s function reveals that this sensitivity originates from the divergence of the nonlocal Green’s function associated with the formation of nonlocal bound states. Our study establishes a framework for realizing novel higher-dimensional non-Hermitian physics in Hermitian bosonic platforms such as superconducting circuits, photonic lattices, and optomechanical arrays, with the demonstrated ultra spectral sensitivity suggesting potential applications in quantum sensing and amplification.

Graphical abstract

Keywords

non-Hermitian skin effects / ultra spectral sensitivity / quantum squeezing / impurity

1 Introduction

In recent years, non-Hermitian Hamiltonians have attracted considerable interest for their ability to host exotic physical phenomena absent in their Hermitian counterparts [1–32]. These effects have been extensively investigated across diverse classical and quantum platforms, including ultracold atoms [33], superconducting qubits [34], quantum [35] and electrical [36] circuits, optical waveguides [37], and acoustic metamaterials [38]. Among the defining features of non-Hermitian systems is the non-Hermitian skin effect (NHSE), characterized by the accumulation of a macroscopic number of eigenstates at the boundaries [39–83], giving rise to an anomalous spectral response that is extraordinarily sensitive to boundary conditions. The NHSE has been shown to give rise to a variety of critical phenomena [84–87] and anomalous dynamical behaviors [88–93] in one-dimensional (1D) systems, lying beyond the scope of the conventional Bloch band description. These effects are now well understood within the generalized Bloch framework, formulated in terms of the generalized Brillouin zone (GBZ) [39–42].

The NHSE has been thoroughly explored in 1D systems. When generalized to higher dimensions, however, recent studies have revealed an even richer set of unconventional phenomena. Notable examples include the geometry-dependent skin effect [94–97], where boundary accumulation is dictated by the system’s geometry, and the algebraic NHSE [42, 98, 99], in which skin modes exhibit quasi-long-range behavior with power-law decaying spatial profiles, in sharp contrast to the exponential localization characteristic of 1D systems. Uniquely emerging in higher-dimensional non-Hermitian systems, the algebraic NHSE leads to striking consequences, such as ultra spectral sensitivity to impurities, wherein the eigenenergy spectrum is dramatically altered by two spatially separated impurities with infinitesimally weak on-site potentials [100]. This sensitive effect arises from the formation of non-local bound states between the impurities and is closely associated with a divergence in the nonlocal back-and-forth Green’s function [100]. In marked contrast to the perturbative robustness inherent to Hermitian systems and 1D nonreciprocal non-Hermitian counterparts, this phenomenon reveals the genuinely new physics that arises only in higher-dimensional non-Hermitian Hamiltonians.

Non-Hermitian Hamiltonians have been widely investigated in both classical and quantum settings, typically through engineered couplings to external dissipative baths [1], which poses significant experimental challenges, particularly in quantum systems. As a complementary approach, bosonic quadratic Hamiltonians with quantum squeezing remain fully Hermitian yet can emulate effective non-Hermitian dynamics, thereby offering a genuine quantum framework to explore non-Hermitian phenomena without relying on external reservoirs [101–118]. In this setting, effective non-Hermiticity arises from the interplay of bosonic statistics and anomalous particle–hole coupling, rather than from an explicitly non-Hermitian microscopic Hamiltonian. Quantum squeezing therefore provides a coherent and tunable resource for engineering non-Hermitian Bogoliubov dynamics, with pairing amplitudes and phases controlled through external parametric drives. Importantly, this Hermitian-parent construction constrains the resulting dynamics through the intrinsic particle–hole and pseudo-Hermitian structures of the bosonic Bogoliubov–de Gennes BdG problem.

To date, investigations of the NHSE in Hermitian quadratic bosonic systems have focused predominantly on one-dimensional geometries. It therefore remains an important open question whether this constrained bosonic BdG mechanism can realize genuinely higher-dimensional non-Bloch phenomena. In particular, can a 2D Hermitian quadratic bosonic Hamiltonian support algebraic skin localization, rather than conventional exponential skin accumulation, and can its Bogoliubov spectrum display a nonperturbative response to asymptotically weak local or nonlocal perturbations? Addressing these questions would establish whether coherent quantum squeezing can reproduce characteristic higher-dimensional non-Hermitian responses within a Hermitian bosonic platform.

In this study, we investigate Hermitian quadratic many-body bosonic Hamiltonians that do not conserve particle number. The bosonic quadratic Hamiltonian is implemented through quantum squeezing in a 2D lattice model. By applying the BdG formalism, the many-body problem is mapped onto an effective single-particle Hamiltonian, whose quasiparticle excitations exhibit a complex eigenvalue spectrum, a hallmark of non-Hermitian physics. Unlike their one-dimensional counterparts, the 2D bosonic quadratic Hamiltonian on a square lattice hosts the algebraic NHSE, characterized by quasi-long-range power-law localization of excitations away from the boundary. Remarkably, although the complex eigenspectrum is insensitive to a single weak on-site impurity, it becomes extremely sensitive to perturbations involving either two infinitesimal on-site impurities or a weak long-range hopping impurity. The manuscript is structured as follows. In Section 2, we construct the 2D Hermitian bosonic quadratic Hamiltonian, incorporating both on-site and off-site quantum squeezing. In Section 3, we demonstrate the emergence of quasi-long-range localization in this 2D Hermitian bosonic quadratic model. Finally, in Section 4, we illustrate the ultra spectral sensitivity to two infinitesimal impurities and elucidate its origin using the Green’s function formalism.

2 Model

We consider a quadratic Hamiltonian describing bosons on a 2D lattice. As schematically illustrated in Fig. 1, the lattice is subjected to on-site and off-site parametric driving. The system Hamiltonian in real space is written as

(1)\begin{array}{rl}{\hat{H}}_{\text{R}}=& \sum _{x,y}{\omega }_0{\hat{a}}_{x,y}^{†}{\hat{a}}_{x,y}+\sum _{x,y}\left(J_{xy}{\hat{a}}_{x+1,y+1}^{†}{\hat{a}}_{x,y}+H.c.\right)\\ & +\sum _{x,y}\left(J_x{\hat{a}}_{x+1,y}^{†}{\hat{a}}_{x,y}+J_y{\hat{a}}_{x,y+1}^{†}{\hat{a}}_{x,y}+H.c.\right)\\ & +\sum _{x,y}\left({\mathrm{\Delta }}_0{\hat{a}}_{x,y}^{†}{\hat{a}}_{x,y}^{†}+{\mathrm{\Delta }}_x{\hat{a}}_{x+1,y}^{†}{\hat{a}}_{x,y}^{†}+H.c.\right),\end{array}

where {\hat{a}}_{x,y}^{†} ({\hat{a}}_{x,y}) denotes the bosonic creation (annihilation) operator at the lattice site (x,y), J_x and J_y represent the nearest-neighbor hopping amplitudes along the x and y directions, respectively, and J_{xy} corresponds to next-nearest-neighbor hopping along the anti-diagonal direction. The parameter {\mathrm{\Delta }}_0 characterizes on-site quantum squeezing, while {\mathrm{\Delta }}_x accounts for off-site quantum squeezing between adjacent sites along the x direction. These terms break global U(1) symmetry, violating particle number conservation. The Hamiltonian {\hat{H}}_{\text{R}} is Hermitian, and its parameters can, in general, take real or complex values. Unless otherwise specified, we set {\omega }_0=0 throughout, which corresponds to working in the rotating frame.

In experimental platforms, the pairing terms in Eq. (1) can be realized in systems such as quantum superconducting circuits [119–121], optomechanical setups [117, 122, 123] and nanophotonic platforms [124, 125]. As a classical analogue, parametrically modulated split-ring-resonator (SRR) lattices provide a possible route for emulating the effective bosonic BdG dynamics considered here, although such systems do not realize genuine quantum squeezing. SRR arrays are particularly suitable because both the magnitude and sign of the electromagnetic coupling can be controlled through the relative orientation, separation, and arrangement of neighboring resonators. In particular, the vertical and diagonal couplings demonstrated in coupled SRR lattices [126] provide the geometric flexibility required to implement the 2D hopping network of our model. The anomalous terms could, in principle, be synthesized by temporally modulating the capacitance, inductance, or resonance frequency of the SRRs. A modulation near twice the resonance frequency generates the classical analogue of an on-site parametric pairing term, whereas modulation of the coupling between two resonators at the sum of their resonance frequencies emulates an intersite two-mode pairing term. The amplitudes and phases of these modulations would determine the corresponding magnitudes and phases of the effective pairing parameters. After expressing the resonator amplitudes in a rotating frame, the coupled-mode equations acquire a Nambu-space structure analogous to the bosonic BdG dynamical equation studied here. Consequently, appropriately designed and modulated SRR arrays could reproduce the algebraic skin localization and impurity-induced spectral reconstruction at the level of classical wave dynamics. We emphasize, however, that this construction is a classical dynamical analogue: it emulates the linear BdG equations but does not generate quantum-squeezed states or their quantum correlations.

The eigenvalue problem associated with the quadratic bosonic Hamiltonian {\hat{H}}_{\text{R}} in Eq. (1) can be reformulated in terms of a BdG formalism [111, 115, 127], as

(2){\hat{H}}_{\text{R}}={\displaystyle \frac{1}{2}}{\hat{\mathrm{\Psi }}}^{†}{\mathcal{H}}_{\text{BdG}}\hat{\mathrm{\Psi }},

where \hat{\mathrm{\Psi }}=(\hat{\mathit{a}},{\hat{\mathit{a}}}^{†}{)}^{\mathrm{T}}, and \hat{\mathit{a}} is a column vector of bosonic annihilation operator, defined as

(3)\hat{\mathit{a}}=\left({\hat{a}}_{1,1},{\hat{a}}_{2,1}\dots ,{\hat{a}}_{L_x,1},{\hat{a}}_{1,2}\dots ,{\hat{a}}_{x,y},\dots ,{\hat{a}}_{L_x,L_y}\right).

Here, L_x and L_y denote the number of lattice sites along x and y directions (see Fig. 1), respectively. The BdG Hamiltonian {\mathcal{H}}_{\text{BdG}} takes the following matrix form with

(4){\mathcal{H}}_{\text{BdG}}=\left(\begin{array}{cc}\mathit{h}& \mathrm{\Delta }\\ {\mathrm{\Delta }}^{†}& {\mathit{h}}^{\text{T}}\end{array}\right),

where \mathit{h}={\mathit{h}}^{†} is a Hermitian matrix describing the single-particle hopping and on-site energy terms, and \mathrm{\Delta }={\mathrm{\Delta }}^{\text{T}} is a symmetric matrix that encodes the pairing interactions. This BdG Hamiltonian exhibits the bosonic analogy of particle-hole symmetry [105] with

(5){\tau }_x{\mathcal{H}}_{\text{BdG}}^{\ast }{\tau }_x={\mathcal{H}}_{\text{BdG}},

where {\tau }_{\mu } (\mu =x,y,z) are generalized Pauli matrices defined in Nambu space with

(6){\tau }_x=\left(\begin{array}{cc}\mathcal{O}& \mathit{I}\\ \mathit{I}& \mathcal{O}\end{array}\right),{\tau }_y=\mathrm{i}\left(\begin{array}{cc}\mathcal{O}& -\mathit{I}\\ \mathit{I}& \mathcal{O}\end{array}\right),{\tau }_z=\left(\begin{array}{cc}\mathit{I}& \mathcal{O}\\ \mathcal{O}& -\mathit{I}\end{array}\right).

Here, \mathcal{O} denotes the zero matrix and \mathit{I} is the identity matrix.

When periodic boundary conditions (PBCs) are imposed, the system can be transformed into momentum space and described by the Bloch Hamiltonian {\hat{H}}_{\text{B}}=\frac{1}{2}\sum \mathit{k}{\hat{\mathrm{\Psi }}}_{\mathit{k}}^{†}{\mathcal{H}}_{\text{BdG}}(\mathit{k}){\hat{\mathrm{\Psi }}}_{\mathit{k}}+C, where {\hat{\mathrm{\Psi }}}_{\mathit{k}}=({\hat{a}}_{\mathit{k}},{\hat{a}}_{-\mathit{k}}^{†}{)}^{\mathrm{T}}, and C=-\text{Tr}[{\mathcal{H}}_0(\mathit{k})]/2 is a constant energy offset. The momentum-space BdG Hamiltonian {\mathcal{H}}_{\text{BdG}}(\mathit{k}) is given by

(7){\mathcal{H}}_{\text{BdG}}(\mathit{k})=\left(\begin{array}{cc}{\mathcal{H}}_0(\mathit{k})& \mathrm{\Delta }(\mathit{k})\\ {\mathrm{\Delta }}^{\ast }(\mathbf{-}\mathit{k})& {\mathcal{H}}_0^{\ast }(\mathbf{-}\mathit{k})\end{array}\right).

Here, {\mathcal{H}}_0(\mathit{k}) is single-particle hopping term, with

(8){\mathcal{H}}_0(\mathit{k})=J_x{\mathrm{e}}^{-\mathrm{i}{k}_x}+J_y{\mathrm{e}}^{-\mathrm{i}{k}_y}+J_{xy}{\mathrm{e}}^{-\mathrm{i}(k_x+k_y)}+H.c.,

and \mathrm{\Delta }(\mathit{k})=2{\mathrm{\Delta }}_0+2{\mathrm{\Delta }}_x\cos k_x denotes the pairing term.

Within this formalism, although {\mathcal{H}}_{\text{BdG}}(\mathit{k}) is Hermitian, the system’s dynamics exhibit a subtle non-Hermitian character [111, 115]. This emergent non-Hermiticity becomes explicit in the Heisenberg dynamics of the bosonic Nambu spinor,

(9)\mathrm{i}{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }t}}{\hat{\mathrm{\Psi }}}_{\mathit{k}}(t)=[{\hat{\mathrm{\Psi }}}_{\mathit{k}}(t),{\hat{H}}_{\mathrm{B}}]={\mathcal{M}}_{\mathrm{B}}(\mathit{k}){\hat{\mathrm{\Psi }}}_{\mathit{k}}(t),

where the dynamical matrix is {\mathcal{M}}_{\mathrm{B}}(\mathit{k})={\sigma }_z{\mathcal{H}}_{\mathrm{B}\mathrm{d}\mathrm{G}}(\mathit{k}). The prefactor {\sigma }_z originates directly from the canonical bosonic commutation relations in Nambu space,

(10)[{\hat{\mathrm{\Psi }}}_{\mathit{k},\alpha },{\hat{\mathrm{\Psi }}}_{{\mathit{k}}^{\prime },\beta }^{†}]=({\sigma }_z{)}_{\alpha \beta }{\delta }_{\mathit{k},{\mathit{k}}^{\prime }},

with {\sigma }_{\mu } denoting the Pauli matrices acting in the particle-hole space. In contrast, for fermionic superconductors the canonical anticommutation relations give \{{\hat{\mathrm{\Psi }}}_{\mathit{k},\alpha },{\hat{\mathrm{\Psi }}}_{{\mathit{k}}^{\prime },\beta }^{†}\}={\delta }_{\alpha \beta }{\delta }_{\mathit{k},{\mathit{k}}^{\prime }}, so that no additional {\sigma }_z metric enters the fermionic BdG equation of motion.

3 Quantum squeezing-induced algebraic non-Hermitian skin effect

The non-Hermitian skin effect in 1D systems is well understood within the framework of non-Bloch band theory [39, 40]. In this framework, a complex deformation of the Bloch wavevector, k\to \tilde{k}=k+\mathrm{i}\mu (\mu \in \mathbb{R}), gives rise to the generalized Brillouin zone (GBZ), where \beta ={\mathrm{e}}^{\mathrm{i}\tilde{k}} encodes the non-Bloch nature of the eigenstates. Under open boundary conditions (OBCs), the eigenstates of 1D non-Hermitian systems take thus the form of modified Bloch waves dressed with an exponential factor {\mathrm{e}}^{\mu }. This results in the characteristic exponential accumulation of eigenstates toward the system boundaries.

A straightforward extension of the 1D concept to higher-dimensional nonreciprocal systems leads to exponentially localized eigenstates along multiple spatial directions, characterized by a single vectorial attenuation factor, e.g., \mathit{\mu }=({\mu }_x,{\mu }_y) in 2D system [81, 128]. However, recent studies have revealed an unusual form of algebraic NHSE in reciprocal non-Hermitian systems [98, 99], where eigenstates exhibit power-law rather than exponential localization. In such systems, the anomalous localization behavior renders the conventional GBZ description, based on a single attenuation vector \mathit{\mu }, inadequate. Instead, a distribution or set of attenuation factors may be needed to fully characterize the system [98, 99]. The algebraic localization further suggests the presence of long-range correlations, offering new opportunities for manipulating quantum states. Despite these recent developments, the manifestation of such algebraic NHSE in Hermitian systems remains largely unexplored. In the following, we investigate the emergence of algebraic NHSE in bosonic quadratic Hermitian systems.

To demonstrate the emergence of the algebraic NHSE in Hermitian systems with the bosonic quadratic Hamiltonian, we analyze the eigenspectrum of its elementary excitation (or quasiparticle), governed by the dynamical matrix {\mathcal{M}}_{\text{B}}={\tau }_z{\mathcal{H}}_{\text{BdG}} in real space. Figures 2(a) and (e) show the complex eigenenergies (orange dots) of {\mathcal{M}}_{\mathrm{B}} under OBCs for the square geometry with different parameters, while the gray regions indicate the spectra under PBCs. The localization behavior of the full eigenspectrum can be quantitatively characterized by evaluating the fractal dimension of the eigenstate {\psi }_i(\mathit{r}) under OBCs [6, 129], defined as

(11)\mathcal{D}[{\psi }_i]=-{\displaystyle \frac{\ln \left[{\displaystyle {\sum }_{\mathit{r}}\left(|{\psi }_{p,i}(\mathit{r}){|}^4+|{\psi }_{h,i}(\mathit{r}){|}^4\right)}\right]}{\ln \sqrt{L_x{L}_y}}},

where \mathit{r}=(x,y), {\psi }_{p,i}(\mathit{r}) and {\psi }_{h,i}(\mathit{r}), in the Nambu space, represent the particle and hole components of the normalized wavefunction {\psi }_i(\mathit{r}) of {\mathcal{M}}_{\text{B}}, and i labels specific eigenstate with eigenvalue E_i. Here, \mathcal{D}[{\psi }_i]=2 signifies a fully extended state, \mathcal{D}[{\psi }_i]=1 indicates strong localization along the boundary, and \mathcal{D}[{\psi }_i]=0 denotes a state localized well at a corner.

Figures 2(b) and (f) present the eigenenergy-resolved fractal dimension \mathcal{D}[{\psi }_i] corresponding to the parameters used in panels (a) and (e), respectively. In both cases, a small fraction of eigenstates correspond to well-localized 1D edge modes with \mathcal{D}[{\psi }_i]\simeq 1, while others form 2D extended states with \mathcal{D}[{\psi }_i]\simeq 2. Strikingly, the majority of eigenstates exhibit intermediate fractal dimension values between 1 and 2, indicating a hybrid localization character associated with quasi-long-range localization. The probability density distributions of eigenstates, given by |{\mathrm{\Psi }}_i(\mathit{r}){|}^2=|{\psi }_{p,i}(\mathit{r}){|}^2+|{\psi }_{h,i}(\mathit{r}){|}^2 with \mathit{r}=(x,y), are shown in Figs. 2(c) and (g). The chosen eigenstates have fractal dimension with values \mathcal{D}[{\psi }_i]\in [1.4,1.6], with their spatial profiles exhibiting long tails along the x direction. This quasi-long-range localization is further confirmed by examining the layer density of eigenstates along the x direction [94], defined as

(12){\mathcal{P}}_i(x)={\displaystyle \sum _y\left(|{\psi }_{p,i}(x,y){|}^2+|{\psi }_{h,i}(x,y){|}^2\right).}

Figures 2(d) and (h) display the layer-resolved densities {\mathcal{P}}_1(x) and {\mathcal{P}}_2(x) on a logarithmic scale (red dots) for the eigenstate with energy E_1 from Fig. 2(a) and E_2 from Fig. 2(e). For the eigenstate with E_1, a power-law decay fit to {\mathcal{P}}_1(x) reveals a quasi-long-range localization of algebraic form rather than exponential. For the eigenstate with E_2, a comparison between exponential (red line) and power-law (blue line) fits reveals a crossover behavior: near the boundary, the decay follows an exponential form, whereas farther away, the profile develops a power-law tail. This power-law tail signifies quasi-long-range localization. We note that the quasi-long-range localization remains robust against boundary disorder and geometric roughness, as demonstrated in Appendix A. Furthermore, the finite-size analysis presented in Appendix B shows that the observed slowly decaying tail is not merely a finite-size artifact.

To further examine the quasi-long-range localization in different lattice geometry, we calculate the fractal dimension \mathcal{D}[{\psi }_i] and probability density |\mathrm{\Psi }(\mathit{r}){|}^2 for typical eigenstates in two oblique-square lattices with different tilt angles, shown in Figs. 2(i, j) and (k, l). The intermediate fractal dimension values \mathcal{D}[{\psi }_i]\in [1.4,1.6] indicate that strong quasi-long-range localization persists for different lattice geometries. These results demonstrate the emergence of the algebraic NHSE independent of lattice geometry in bosonic quadratic Hermitian systems.

4 Ultra spectral sensitivity to impurities

4.1 Numerical results

We now demonstrate the ultra spectral sensitivity to impurities in bosonic quadratic Hermitian systems, a phenomenon that was previously identified only in higher-dimensional non-Hermitian systems [100].

The quasiparticle excitation spectrum is strongly influenced by boundary conditions, as seen by the contrast between Figs. 3(a, e) under PBCs along both x and y directions and (b, f) under OBC in x direction and PBC in y direction. Even more strikingly, the eigenspectrum of the bosonic quadratic Hermitian system shows remarkable sensitivity to impurities with sufficiently weak strength. To illustrate this, we introduce on-site impurities at positions {\mathit{r}}_1 and {\mathit{r}}_2, where the Hamiltonian takes the form \hat{H}={\hat{H}}_{\text{R}}+\hat{V}, with

(13)\hat{V}=V_1{\hat{a}}_{{\mathit{r}}_1}^{†}{\hat{a}}_{{\mathit{r}}_1}+V_2{\hat{a}}_{{\mathit{r}}_2}^{†}{\hat{a}}_{{\mathit{r}}_2},

where V_1 and V_2 denote the on-site potentials. The full Hamiltonian of the impurity lattice, expressed in the Nambu spinor basis, is written as

(14)\hat{H}={\hat{H}}_{\text{R}}+\hat{V}={\displaystyle \frac{1}{2}}{\hat{\mathrm{\Psi }}}^{†}\mathcal{H}\hat{\mathrm{\Psi }},

where \mathcal{H}={\mathcal{H}}_{\text{BdG}}+\mathcal{V}, with \mathcal{V} being the BdG representation of the impurity potential. Unless otherwise specified, here, we consider OBC along the x direction and PBC along the y direction.

We plot the quasiparticle excitation spectra for two sets of parameters in Figs. 3(c, d) and (g, h), corresponding to a single impurity with sufficiently weak strength [V_1=0.01, V_2=0, see Figs. 3(c, g)] and two distant impurities with sufficiently weak strength [V_1=V_2=0.01, see Figs. 3(d, h)]. Whereas a single impurity with sufficiently weak on-site potential leaves the eigenspectrum almost unchanged [see Figs. 3(c, g)], introducing two impurities with sufficiently weak strength produces a marked alteration of the quasiparticle excitations [see Figs. 3(d, h)], as evident from the comparison with the single-impurity case. Specifically, the two infinitesimal impurities not only create additional eigenstates [black arrows in Figs. 3(d, h)], but also eliminate existing ones, such as the state indicated by the red arrow in Fig. 3(d). This behavior demonstrates that the bosonic quadratic Hermitian system, studied here, exhibits an ultra spectral sensitivity to infinitesimal impurities, a feature that is generally absent in particle-conserving Hermitian systems.

In addition to the on-site impurity, the bosonic quadratic Hermitian system can also exhibit ultra spectral sensitivity to a long-range hopping impurity. To illustrate this, we consider the following form of long-range hopping impurity

(15)\hat{V}=t_{\mathrm{p}}({\hat{a}}_{x_1,y_1}^{†}{\hat{a}}_{x_2,y_2}+H.c.).

With OBC imposed along the x direction and PBC along the y direction, we compute the complex quasiparticle eigenenergies by fixing the impurity position at (x_1,y_1)=(1,1) and varying the second site (x_2,y_2) to explore different hopping ranges. The results are shown in Fig. 4(a) for (x_2,y_2)=(3L_x/5,L_y/2) and in Fig. 4(b) for (x_2,y_2)=(L_x,L_y/2), both obtained under a sufficiently weak hopping strength t_{\text{p}}. Despite the weakness of the hopping strength, the eigenspectrum exhibits impurity-induced states, indicated by the black arrows in Figs. 4(a, b), once the hopping range exceeds a critical length [see Fig. 4(a)].

4.2 Physical mechanism

To gain deeper insight into the ultra spectral sensitivity of bosonic quadratic Hermitian systems to impurities, we employ the Green’s function formalism, which offers a powerful framework for analyzing their response to perturbations [100].

The eigenvalue equation of the full system in the presence of impurities is written as

(16)(E{\tau }_0-\mathcal{M})\mathrm{\Psi }=0,

where {\tau }_0 is identity matrix, and \mathcal{M}={\tau }_z\mathcal{H}. The Green’s function in a bosonic BdG system can be defined as the inverse of this linear operator E{\tau }_0-\mathcal{M}. Explicitly, the full Green’s function in real space takes the form as

(17)\mathcal{G}=(E{\tau }_0-\mathcal{M}{)}^{-1}=({\tau }_0-{\mathcal{G}}_0{\mathcal{V}}_z{)}^{-1}{\mathcal{G}}_0,

where {\mathcal{V}}_z={\tau }_z\mathcal{V}, and {\mathcal{G}}_0=(E{\tau }_0-{\mathcal{M}}_{\text{BdG}}{)}^{-1} denotes the Green’s function of the pristine lattice in the absence of impurities. In this work, we restrict our analysis to energies E outside the spectrum of {\mathcal{M}}_{\text{BdG}}, for which {\mathcal{G}}_0 remains generally non-singular.

The effect of sufficiently weak impurities is captured by ({\tau }_0-{\mathcal{G}}_0{\mathcal{V}}_z{)}^{-1} in Eq. (17), which renormalizes the unperturbed propagator {\mathcal{G}}_0. When all eigenvalues of the response matrix {\mathcal{G}}_0{\mathcal{V}}_z remain much smaller than unity, \mathcal{G} stays close to {\mathcal{G}}_0 and the perturbed spectrum of \mathcal{M} is essentially unchanged. In contrast, once even a single eigenvalue approaches or exceeds unity, the system enters a non-perturbative regime in which both the Green’s function \mathcal{G} and the associated spectrum of \mathcal{M} exhibit strong departures from their unperturbed counterparts {\mathcal{G}}_0 and {\mathcal{M}}_{\text{BdG}}.

By defining the spectral radius \rho of {\mathcal{G}}_0{\mathcal{V}}_z [100], a quantitative criterion for the spectral stability against infinitesimal impurities can be formulated as

(18)\rho ({\mathcal{G}}_0{\mathcal{V}}_z)=max|\sigma ({\mathcal{G}}_0{\mathcal{V}}_z)|\ll 1,

where \sigma (\cdot ) denotes the set of eigenvalues. Violation of this criterion can lead to non-perturbative changes, even for weak impurities.

To obtain an analytical result for the Green’s function \mathcal{G} in Eq. (17) and gain an intuitive understanding of the ultra spectral sensitivity, we set J_x=0, J_y=\mathrm{i}{t}_y, and J_{xy}=\mathrm{i}{t}_{xy} with real pairing amplitudes {\mathrm{\Delta }}_0 and {\mathrm{\Delta }}_x, where t_y and t_{xy} are real. Under these conditions, the dynamical matrix {\mathcal{M}}_{\text{BdG}} can be unitarily transformed into a block-diagonal form (see details in Appendix C.1) with

(19)U{\mathcal{M}}_{\text{BdG}}{U}^{†}=\left(\begin{array}{cc}{\mathcal{M}}_p& \mathcal{O}\\ \mathcal{O}& {\mathcal{M}}_m\end{array}\right),

where U=({\tau }_0-\mathrm{i}{\tau }_x)/\sqrt{2} with U^{†}=U^{-1}, {\mathcal{M}}_p and {\mathcal{M}}_m represent the decoupled matrices in real space.

We now apply the same unitary transformation to the dynamical matrix \mathcal{M} in the presence of impurities, yielding \overline{\mathcal{M}}=U\mathcal{M}{U}^{†}. Under this transformation, the Green’s function \mathcal{G} in Eq. (17) becomes

(20)\overline{\mathcal{G}}=(E{\tau }_0-\overline{\mathcal{M}}{)}^{-1}=({\tau }_0-{\overline{\mathcal{G}}}_0{\overline{\mathcal{V}}}_z{)}^{-1}{\overline{\mathcal{G}}}_0,

where {\overline{\mathcal{G}}}_0=(E{\tau }_0-U{\mathcal{M}}_{\text{BdG}}{U}^{†}{)}^{-1}, and {\overline{\mathcal{V}}}_z=U{\mathcal{V}}_z{U}^{†}.

The unitary transformation of the Green’s function in Eq. (20) gives rise to a modified spectral stability criterion for infinitesimal impurities in Eq. (18), expressed via the reduced spectral radius with

(21)\rho ({\overline{\mathcal{G}}}_0{\overline{\mathcal{V}}}_z)\ll 1.

As shown in the Appendix C.2, for a single on-site impurity with \hat{V}=V_1{\hat{a}}_{{\mathit{r}}_1}^{†}{\hat{a}}_{{\mathit{r}}_1}, the spectral radius is obtained as

(22)\rho ({\overline{\mathcal{G}}}_0{\overline{\mathcal{V}}}_z)=|\sqrt{\xi }|=|V_1\sqrt{{\overline{\mathcal{G}}}_0^p({\mathit{r}}_1,{\mathit{r}}_1;E){\overline{\mathcal{G}}}_0^m({\mathit{r}}_1,{\mathit{r}}_1;E)}|,

where the unperturbed Green’s function {\overline{\mathcal{G}}}_0^{p/m}({\mathit{r}}_i,{\mathit{r}}_j;E) is given by

(23){\overline{\mathcal{G}}}_0^{p/m}({\mathit{r}}_i,{\mathit{r}}_j;E)=⟨{\mathit{r}}_i|(E-{\mathcal{M}}_{p/m}{)}^{-1}|{\mathit{r}}_j⟩.

We explicitly evaluate the Green’s function {\overline{\mathcal{G}}}_0^{p/m}({\mathit{r}}_i,{\mathit{r}}_j;E) under OBC along the x direction and PBC along the y direction for the cylindrical geometry. To capture the NHSE and perform the momentum-space integration in GBZ, we introduce the non-Bloch wavevector {\beta }_x={\mathrm{e}}^{\mathrm{i}{k}_x+{\mu }_x} along the x direction, while keeping the Bloch wavevector k_y along the y direction [100]. The Green’s function then reads (see Appendix D)

(24)\begin{array}{rl}& {\overline{\mathcal{G}}}_0^{p/m}({\mathit{r}}_i,{\mathit{r}}_j;E)\\ & ={\int }_0^{2\pi }\mathrm{d}{k}_y{\oint }_{{\beta }_x\in \text{GBZ}}\frac{\text{d}{\beta }_x}{{\beta }_x}\frac{{\mathrm{e}}^{\mathrm{i}{k}_y(y_i-y_j)}{\beta }_x^{x_i-x_j}}{E-E_{\pm }({\beta }_x,k_y)},\end{array}

where E_{\pm }({\beta }_x,k_y) reads

(25)\begin{array}{rl}{E}_{\pm }({\beta }_x,k_y)=& \pm \mathrm{i}{\mathrm{\Delta }}_x\left({\beta }_x+{\beta }_x^{-1}\right)+\mathrm{i}{t}_y\left({\mathrm{e}}^{\mathrm{i}{k}_y}-{\mathrm{e}}^{-\mathrm{i}{k}_y}\right)\\ & \pm 2\mathrm{i}{\mathrm{\Delta }}_0+\mathrm{i}{t}_{xy}\left({\beta }_x{\mathrm{e}}^{\mathrm{i}{k}_y}-{\beta }_x^{-1}{\mathrm{e}}^{-\mathrm{i}{k}_y}\right)\\ & +\mathrm{i}{t}_y\left({\mathrm{e}}^{\mathrm{i}{k}_y}-{\mathrm{e}}^{-\mathrm{i}{k}_y}\right).\end{array}

According to Eq. (24), the local Green’s functions {\overline{\mathcal{G}}}_0^{p/m}({\mathit{r}}_1,{\mathit{r}}_1;E) in Eq. (22) remain finite because the energy E is chosen off-resonant from all eigenstates of the unperturbed dynamical matrix {\mathcal{M}}_{\text{BdG}}. Consequently, the spectral radius in Eq. (22) vanishes in the limit V_1\to 0, corresponding to an infinitesimal impurity. This indicates that the system’s spectrum remains robust against a single sufficiently-weak impurity perturbation, with no emergence of additional eigenvalues, consistent with the numerical results in Figs. 3(c, g).

Furthermore, in the presence of double sufficiently-weak on-site impurities with \hat{V}=V_1{\hat{a}}_{{\mathit{r}}_1}^{†}{\hat{a}}_{{\mathit{r}}_1}+V_2{\hat{a}}_{{\mathit{r}}_2}^{†}{\hat{a}}_{{\mathit{r}}_2}, the spectral radius (see Appendix C.2) is written as

(26)\rho ({\overline{\mathcal{G}}}_0{\overline{\mathcal{V}}}_z)=max(|\sqrt{{\xi }_{\pm }}|),

where

(27){\xi }_{\pm }={\displaystyle \frac{1}{2}}\left({\mathcal{B}\mathcal{C}}_{11}+{\mathcal{B}\mathcal{C}}_{22}\pm \sqrt{({\mathcal{B}\mathcal{C}}_{11}-{\mathcal{B}\mathcal{C}}_{22}{)}^2+4{\mathcal{B}\mathcal{C}}_{12}{\mathcal{B}\mathcal{C}}_{21}}\right),

and \mathcal{B}\mathcal{C} is the 2\times 2 matrix, with its elements given by

(28)\begin{array}{rl}{\mathcal{B}\mathcal{C}}_{11}=& V_1^2{\overline{\mathcal{G}}}_0^p({\mathit{r}}_1,{\mathit{r}}_1;E){\overline{\mathcal{G}}}_0^m({\mathit{r}}_1,{\mathit{r}}_1;E)\\ & +V_1{V}_2{\overline{\mathcal{G}}}_0^p({\mathit{r}}_1,{\mathit{r}}_2;E){\overline{\mathcal{G}}}_0^m({\mathit{r}}_2,{\mathit{r}}_1;E),\end{array}

(29)\begin{array}{rl}{\mathcal{B}\mathcal{C}}_{12}=& V_1{V}_2{\overline{\mathcal{G}}}_0^p({\mathit{r}}_1,{\mathit{r}}_1;E){\overline{\mathcal{G}}}_0^m({\mathit{r}}_1,{\mathit{r}}_2;E)\\ & +V_2^2{\overline{\mathcal{G}}}_0^p({\mathit{r}}_1,{\mathit{r}}_2;E){\overline{\mathcal{G}}}_0^m({\mathit{r}}_2,{\mathit{r}}_2;E),\end{array}

(30)\begin{array}{rl}{\mathcal{B}\mathcal{C}}_{21}=& V_1^2{\overline{\mathcal{G}}}_0^p({\mathit{r}}_2,{\mathit{r}}_1;E){\overline{\mathcal{G}}}_0^m({\mathit{r}}_1,{\mathit{r}}_1;E)\\ & +V_1{V}_2{\overline{\mathcal{G}}}_0^p({\mathit{r}}_2,{\mathit{r}}_2;E){\overline{\mathcal{G}}}_0^m({\mathit{r}}_2,{\mathit{r}}_1;E),\end{array}

(31)\begin{array}{rl}{\mathcal{B}\mathcal{C}}_{22}=& V_1{V}_2{\overline{\mathcal{G}}}_0^p({\mathit{r}}_2,{\mathit{r}}_1;E){\overline{\mathcal{G}}}_0^m({\mathit{r}}_1,{\mathit{r}}_2;E)\\ & +V_2^2{\overline{\mathcal{G}}}_0^p({\mathit{r}}_2,{\mathit{r}}_2;E){\overline{\mathcal{G}}}_0^m({\mathit{r}}_2,{\mathit{r}}_2;E).\end{array}

As shown in Eqs. (26)−(31), unlike the single impurity case, the spectral radius for two impurities depend not only on the local Green’s function {\overline{\mathcal{G}}}_0^{p/m}({\mathit{r}}_i,{\mathit{r}}_i;E) (i=1,2), but also on the propagators {\overline{\mathcal{G}}}_0^{p/m}({\mathit{r}}_i,{\mathit{r}}_j;E) with i\ne j.

To carry out the contour integral for the propagator in Eq. (24), we solve characteristic equations E-E_{\pm }({\beta }_x,k_y)=0 in the subspace defined in Eq. (19), where E is outside the eigenspectrum of {\mathcal{M}}_{p/m}. The corresponding roots of the characteristic equations, {\beta }_{1,2}^{\pm }(k_y), are ordered such that |{\beta }_2^{\pm }(k_y)|>|{\beta }_1^{\pm }(k_y)|. Then, using the residue theorem and the thermodynamic limit, the propagators {\overline{\mathcal{G}}}_0^{p/m}({\mathit{r}}_i,{\mathit{r}}_j;E) (i\ne j) are approximated (see details in Appendix D) as

(32){\overline{\mathcal{G}}}_0^{p/m}({\mathit{r}}_2,{\mathit{r}}_1;E)\simeq {\displaystyle {\int }_0^{2\pi }\mathrm{d}{k}_y[{\beta }_1^{\pm }(k_y){]}^L{\mathrm{e}}^{\mathrm{i}{k}_y{\delta }_y},}

(33){\overline{\mathcal{G}}}_0^{p/m}({\mathit{r}}_1,{\mathit{r}}_2;E)\simeq {\displaystyle {\int }_0^{2\pi }\mathrm{d}{k}_y[{\displaystyle \frac{1}{{\beta }_2^{\pm }(k_y)}}{]}^L{\mathrm{e}}^{-\mathrm{i}{k}_y{\delta }_y},}

where L=x_2-x_1 and {\delta }_y=y_2-y_1.

Our goal is to estimate the spectral radius \rho ({\overline{\mathcal{G}}}_0{\overline{\mathcal{V}}}_z) in Eq. (26) for a large lattice size along the x direction, and to demonstrate that this quantity deviates from the spectral stability criterion in Eq. (21) in the presence of two infinitesimal on-site impurities, thereby revealing the ultra spectral sensitivity of Hermitian bosonic quadratic systems. For large L, the integrals in Eqs. (32) and (33) can be approximated by choosing values of {\delta }_y that cancel the oscillatory phase in the integrands (see Appendix D). These stationary-phase points provide the dominant contributions to the integrals.

For the new eigenvalue E with positive imaginary part, for instance, the eigenvalue E indicated by the orange circle in Fig. 3(h), the integrals in Eqs. (32) and (33) can be estimated as

(34)\{\begin{array}{l}{\overline{\mathcal{G}}}_0^p({\mathit{r}}_2,{\mathit{r}}_1;E)\sim {\mathrm{e}}^{{\mu }_{\text{max,1}}^{+}L},\\ {\overline{\mathcal{G}}}_0^m({\mathit{r}}_2,{\mathit{r}}_1;E)\sim {\mathrm{e}}^{{\mu }_{\text{max,1}}^{-}L}\sim 0,\\ {\overline{\mathcal{G}}}_0^p({\mathit{r}}_1,{\mathit{r}}_2;E)\sim {\mathrm{e}}^{-{\mu }_{\text{min,2}}^{+}L},\\ {\overline{\mathcal{G}}}_0^m({\mathit{r}}_1,{\mathit{r}}_2;E)\sim {\mathrm{e}}^{-{\mu }_{\text{min,2}}^{-}L}\sim 0,\end{array}

with

(35){\mu }_{\text{max,1}}^{\pm }=\underset{k_y\in [0,2\pi ]}{max}\ln |{\beta }_1^{\pm }(k_y)|,

and

(36){\mu }_{\text{min,2}}^{\pm }=\underset{k_y\in [0,2\pi ]}{min}\ln |{\beta }_2^{\pm }(k_y)|.

Here, {\mu }_{\text{max},1}^{+}>0, , , and {\mu }_{\text{min},2}^{-}>0 are typically satisfied in the parameter regime of interest [see red and blue dots in Figs. 5(a1) and (a3)], thereby facilitating the demonstration of the ultra spectral sensitivity.

Therefore, for double on-site impurities, the spectral radius associated with the induced new eigenvalue E with positive imaginary part can be approximated by

(37)\rho ({\overline{\mathcal{G}}}_0{\overline{\mathcal{V}}}_z)\sim V{\mathrm{e}}^{\frac{{\mu }_{\text{max,1}}^{+}L}{4}}{\mathrm{e}}^{-\frac{{\mu }_{\text{min,2}}^{+}L}{4}}.

Here, we consider the case where the two on-site impurities have identical on-site potentials with V_1=V_2=V.

For the new eigenvalue E with negative imaginary part, for instance, the eigenvalue E marked by the red circle in Fig. 3(h), the integrals in Eqs. (32) and (33) can be estimated as

(38)\{\begin{array}{l}{\overline{\mathcal{G}}}_0^p({\mathit{r}}_2,{\mathit{r}}_1;E)\sim {\mathrm{e}}^{{\mu }_{\text{max,1}}^{+}L}\sim 0,\\ {\overline{\mathcal{G}}}_0^m({\mathit{r}}_2,{\mathit{r}}_1;E)\sim {\mathrm{e}}^{{\mu }_{\text{max,1}}^{-}L},\\ {\overline{\mathcal{G}}}_0^p({\mathit{r}}_1,{\mathit{r}}_2;E)\sim {\mathrm{e}}^{-{\mu }_{\text{min,2}}^{+}L}\sim 0,\\ {\overline{\mathcal{G}}}_0^m({\mathit{r}}_1,{\mathit{r}}_2;E)\sim {\mathrm{e}}^{-{\mu }_{\text{min,2}}^{-}L},\end{array}

where , {\mu }_{\text{max},1}^{-}>0, {\mu }_{\text{min},2}^{+}>0, and are typically satisfied in the parameter regime of interest [see red and blue dots in Figs. 5(b1) and (b3)]. The corresponding spectral radius associated with the induced new eigenvalue E with negative imaginary part is

(39)\rho ({\overline{\mathcal{G}}}_0{\overline{\mathcal{V}}}_z)\sim V{\mathrm{e}}^{\frac{{\mu }_{\text{max,1}}^{-}L}{4}}{\mathrm{e}}^{-\frac{{\mu }_{\text{min,2}}^{-}L}{4}},

where V_1=V_2=V is assumed.

In contrast to the single weak on-site impurity case, the presence of two infinitesimal on-site impurities can cause the propagators to diverge for large lattice sizes [see Eqs. (34) and (38)]. Consequently, the spectral radius in Eqs. (37) and (39) deviates from the spectral stability criterion in Eq. (21), thereby revealing the ultra spectral sensitivity of Hermitian bosonic quadratic systems. The underlying mechanism is directly analogous to that uncovered in the non-Hermitian case discussed in Ref. [100]: the nonperturbative spectral response originates from the enhancement of nonlocal Green’s-function propagation between spatially separated impurities.

This ultraspectral response can be probed through linear-response spectroscopy. For a weak coherent probe and an approximately uniform damping rate \gamma, the retarded susceptibility is

(40){\chi }^R(\omega )={\left[(\omega +\mathrm{i}\gamma /2){\tau }_0-\mathcal{M}\right]}^{-1},

where \mathcal{M} is the bosonic dynamical matrix. The impurity-induced spectral reconstruction modifies the poles of {\chi }^R(\omega ) and can therefore manifest as the emergence, disappearance, or displacement of resonances, together with pronounced changes in transmission, reflection, spatial response, and phase-sensitive amplification. In superconducting-circuit, photonic, and optomechanical implementations, a direct protocol is to compare the pristine, single-impurity, and double-impurity configurations. Whereas one weak impurity produces only a perturbative response, adding a second distant impurity can strongly reconstruct the resonance spectrum and enhance the nonlocal response between the impurity positions. These changes are directly resolvable when they exceed the relevant linewidth; otherwise, they may remain detectable through the transmission phase or spatially resolved susceptibility. Spatially resolved measurements can also reveal the slowly decaying boundary profile associated with the algebraic NHSE.

5 Conclusions

In this work, we demonstrate the algebraic NHSE and non-perturbative spectral sensitivity to impurities in Hermitian bosonic quantum many-body systems. We construct a 2D Hermitian bosonic quadratic Hamiltonian on a square lattice, realized through both on-site and off-site quantum squeezing. By numerically investigating the bosonic excitations within the BdG framework, we reveal a quasi-long-range power-law localization of complex eigenstates, in stark contrast to the exponential localization observed in the 1D counterpart. Furthermore, the 2D Hermitian bosonic quadratic Hamiltonian is insensitive to infinitesimal single impurities, yet exhibits ultra spectral sensitivity to either double infinitesimal on-site impurities or a weak long-range hopping impurity. Specifically, for double infinitesimal on-site impurities, the eigenspectrum either develops new eigenstates or causes certain eigenstates to vanish. Using the Green’s function formalism, we analytically uncover the mechanism underlying this ultra spectral sensitivity, which arises from the formation of nonlocal bound states between impurities, leading to the divergence of the nonlocal Green’s function.

6 Appendix A: Effects of boundary disorder and roughness on quasi-long-range localization

To determine whether the algebraic non-Hermitian skin effect depends on an ideally clean and smooth boundary, we have performed additional numerical calculations for two representative classes of boundary imperfections: random boundary on-site potentials and geometric boundary roughness. In both cases, the Hamiltonian in the interior of the sample is left unchanged.

First, we introduce random on-site potentials on the outermost boundary,

(A1)\begin{array}{c}{\hat{H}}_{\mathrm{o}\mathrm{n}}={\displaystyle \sum _{\mathit{r}\in \mathrm{\partial }\mathrm{\Omega }}{V}_{\mathit{r}}{\hat{a}}_{\mathit{r}}^{†}{\hat{a}}_{\mathit{r}},V_{\mathit{r}}=W_{\mathrm{o}\mathrm{n}}{\xi }_{\mathit{r}},}\end{array}

where the independent random variables {\xi }_{\mathit{r}} are uniformly distributed over [-1,1], and W_{\mathrm{o}\mathrm{n}} controls the disorder strength.

Second, we model geometric boundary roughness by randomly removing sites from the outermost boundary. For each boundary site, we introduce a binary variable {\chi }_{\mathit{r}}\in \{0,1\} satisfying

(A2)\begin{array}{c}\mathrm{P}\mathrm{r}\mathrm{o}\mathrm{b}({\chi }_{\mathit{r}}=0)=p_{\mathrm{r}},\mathrm{P}\mathrm{r}\mathrm{o}\mathrm{b}({\chi }_{\mathit{r}}=1)=1-p_{\mathrm{r}},\end{array}

where p_{\mathrm{r}} is the removal probability. The hopping and pairing amplitudes are then modified according to

(A3)\begin{array}{c}{J}_{\mathit{r}{\mathit{r}}^{\prime }}\to {\chi }_{\mathit{r}}{\chi }_{{\mathit{r}}^{\prime }}{J}_{\mathit{r}{\mathit{r}}^{\prime }},{\mathrm{\Delta }}_{\mathit{r}{\mathit{r}}^{\prime }}\to {\chi }_{\mathit{r}}{\chi }_{{\mathit{r}}^{\prime }}{\mathrm{\Delta }}_{\mathit{r}{\mathit{r}}^{\prime }}.\end{array}

In the numerical implementation, each removed site and all bonds connected to it are deleted directly from the real-space dynamical matrix. This procedure produces a genuinely irregular physical boundary, rather than representing a missing site through an artificially large local potential.

For each disorder realization l, we construct the corresponding finite-size bosonic dynamical matrix and calculate N_{\mathrm{e}\mathrm{v}}=16 eigenmodes in the vicinity of the target complex energy

(A4)\begin{array}{c}{E}_{\mathrm{t}\mathrm{a}\mathrm{r}}=2.09+9.23\mathrm{i}.\end{array}

Because boundary disorder shifts and rearranges the complex eigenspectrum, following a fixed eigenvalue index of the clean system would not provide a reliable state-by-state comparison. We therefore select the mode whose eigenenergy is closest to the target value with

(A5)\begin{array}{c}{n}_l^{\ast }=\underset{n\in {\mathcal{S}}_l}{\text{argmin}}\left|E_{n,l}-E_{\mathrm{t}\mathrm{a}\mathrm{r}}\right|,\end{array}

where {\mathcal{S}}_l denotes the set of the 16 candidate modes in the lth disorder realization. The selected right eigenvector is normalized as

(A6)\begin{array}{c}⟨{\psi }_{n_l^{\ast }}^R|{\psi }_{n_l^{\ast }}^R⟩=1,\end{array}

consistent with the normalization adopted for the clean system. Its layer-resolved density is then calculated as

(A7)\begin{array}{c}{\mathcal{P}}_l(x)={\displaystyle \sum _y\left[|{\psi }_{p,l}(x,y){|}^2+|{\psi }_{h,l}(x,y){|}^2\right].}\end{array}

Because the skin modes accumulate near two opposite boundaries, we characterize their decay using the distance from the nearest boundary,

(A8)\begin{array}{c}d=min(x,L_x+1-x),\end{array}

and combine the densities of layers having the same value of d. For notational simplicity, this boundary distance is denoted by x in Fig. A1. The disorder-averaged profile is

(A9)\begin{array}{c}\overline{\mathcal{P}}(x)={\displaystyle \frac{1}{N_{\mathrm{d}\mathrm{i}\mathrm{s}}}}{\displaystyle \sum _{l=1}^{N_{\mathrm{d}\mathrm{i}\mathrm{s}}}{\mathcal{P}}_l(x),N_{\mathrm{d}\mathrm{i}\mathrm{s}}=100.}\end{array}

The error bars in Figs. A1(a, c) quantify the sample-to-sample fluctuations among the 100 independent disorder realizations and therefore represent physical fluctuations induced by boundary disorder, rather than numerical uncertainties of the eigensolver.

Figure A1(a) shows the disorder-averaged profiles for boundary on-site disorder with W_{\mathrm{o}\mathrm{n}}=0.3, 0.6, and 0.9. The disorder modifies the short-distance density and increases the realization-to-realization fluctuations. Nevertheless, even at the largest disorder strength considered, the profile retains a pronounced slowly decaying tail and deviates systematically from the straight-line behavior expected for purely exponential localization on the semilogarithmic plot. The corresponding representative real-space density distribution at W_{\mathrm{o}\mathrm{n}}=0.9, shown in Fig. A1(b), continues to exhibit boundary accumulation accompanied by substantial penetration into the interior.

Figure A1(c) presents the corresponding results for geometric boundary roughness with removal probabilities p_{\mathrm{r}}=0.01, 0.05, and 0.10. Increasing p_{\mathrm{r}} produces stronger spatial fluctuations and a more pronounced deformation of the disorder-averaged profile, particularly at large distances from the boundary. At p_{\mathrm{r}}=0.10, the large error bars reflect substantial realization-to-realization variations caused by the irregular boundary geometry. Despite these fluctuations, the profile does not cross over to conventional short-ranged exponential decay. Instead, an extended non-exponential tail remains visible. The representative density distribution at p_{\mathrm{r}}=0.10 in Fig. A1(d) likewise shows that the mode remains preferentially accumulated near the irregular boundary while maintaining appreciable weight deep inside the sample.

These results demonstrate that the long-ranged, non-exponential character of the skin localization is not an artifact of an ideally clean and geometrically smooth boundary. Boundary on-site disorder mainly perturbs the local density near the boundary, whereas geometric roughness causes stronger mode mixing, spatial inhomogeneity and sample-to-sample fluctuations.

7 Appendix B: Finite-size effect

To examine whether the slowly decaying tail arises merely from a finite-size crossover, we have performed an additional longitudinal finite-size analysis for L_x=80,100,120,160,\text{and}200, while keeping the transverse size fixed at L_y=40. The results are summarized in Fig. A2.

Figures A2(a–e) show the complex eigenspectra color coded by the corresponding fractal dimensions. The overall eigenenergy-resolved fractal-dimension structure remains qualitatively unchanged as L_x increases. To compare spatial profiles associated with the same spectral branch, for each size we select an eigenstate from a common neighborhood of the complex-energy plane, as marked by the black dots. The selected eigenenergies vary only weakly, from E=1.853+9.255\mathrm{i} at L_x=80 to E=1.848+9.242\mathrm{i} at L_x=200.

The corresponding layer-resolved densities are presented in Figs. A2(f–j). For every size considered, the profile displays two distinguishable spatial regimes. Close to the left boundary, the density initially follows an approximately exponential decay,

(B1)\begin{array}{c}\mathcal{P}(x)\propto {\mathrm{e}}^{-x/\xi },\end{array}

represented by the red straight lines on the semilogarithmic plots. At larger distances, the numerical profile deviates systematically from this exponential extrapolation and develops a substantially slower tail consistent with

(B2)\begin{array}{c}\mathcal{P}(x)\propto x^{-\alpha }.\end{array}

The blue curves show a common algebraic envelope with \alpha \simeq 0.12 over the resolved long-distance regime.

As L_x increases from 80 to 200, the available spatial window becomes progressively larger. Nevertheless, the numerical profiles do not approach the rapidly decreasing red exponential extrapolation. Instead, the slow tail persists and remains broadly consistent with the algebraic envelope. For the largest sizes, weak spatial oscillations are superimposed on the slowly varying background; these modulations do not restore conventional exponential decay.

These results provide a nontrivial finite-size check of the algebraic-skin interpretation. In particular, the persistence of the slowly decaying tail over an expanding longitudinal interval disfavors the possibility that it is solely a short-distance crossover caused by the smaller lattices. Together with the stable eigenenergy-resolved fractal-dimension pattern, the results are consistent with algebraic, quasi-long-range boundary localization over the system sizes accessible to our calculations.

8 Appendix C: Details on Green’s function and reduced spectral radius

8.1 Green’s function in bosonic BdG system

We begin by writing the full Hamiltonian in the presence of impurities, expressed in the Nambu basis as

(C1)\begin{array}{c}\hat{H}={\hat{H}}_{\text{R}}+\hat{V}={\displaystyle \frac{1}{2}}{\hat{\mathrm{\Psi }}}^{†}\mathcal{H}\hat{\mathrm{\Psi }},\end{array}

where \mathrm{\Psi }=(\hat{\mathit{a}},{\hat{\mathit{a}}}^{†}{)}^{\mathrm{T}}.

According to Eq. (9), the eigenvalue equation of the full system is written as

(C2)\begin{array}{c}(E{\tau }_0-\mathcal{M})\mathrm{\Psi }=0.\end{array}

The Green’s function of the bosonic BdG system is then defined as the inverse of the operator E{\tau }_0-\mathcal{M}. Explicitly, the full Green’s function reads

(C3)\begin{array}{c}\mathcal{G}=(E{\tau }_0-\mathcal{M}{)}^{-1}=({\tau }_0-{\mathcal{G}}_0{\mathcal{V}}_z{)}^{-1}{\mathcal{G}}_0,\end{array}

where {\mathcal{V}}_z={\tau }_z\mathcal{V}, and {\mathcal{G}}_0=(E{\tau }_0-{\mathcal{M}}_{\text{BdG}}{)}^{-1} denotes the Green’s function of the unperturbed system.

To obtain an analytical expression for the Green’s function \mathcal{G} and the spectral radius \rho, and thereby gain an intuitive understanding of the ultra spectral sensitivity, we set J_x=0, J_y=\mathrm{i}{t}_y, and J_{xy}=\mathrm{i}{t}_{xy} with real pairing amplitudes {\mathrm{\Delta }}_0 and {\mathrm{\Delta }}_x, where t_y and t_{xy} are real. Under these conditions, the momentum-space BdG Hamiltonian is written as

(C4)\begin{array}{c}\begin{array}{rl}{\mathcal{H}}_{\text{BdG}}(\mathit{k})=& [2t_y\sin k_y+2t_{xy}\sin (k_x+k_y)]{\sigma }_z\\ & +(2{\mathrm{\Delta }}_0+2{\mathrm{\Delta }}_x\cos k_x){\sigma }_x.\end{array}\end{array}

Consequently, the eigenvalues of the dynamical matrix {\mathcal{M}}_{\text{BdG}}(\mathit{k})={\sigma }_z{\mathcal{H}}_{\text{BdG}}(\mathit{k}), denoted E_{\pm }(\mathit{k}), are readily obtained as

(C5)\begin{array}{c}\begin{array}{rl}{E}_{\pm }(\mathit{k})=& \pm 2\mathrm{i}{\mathrm{\Delta }}_0\pm 2\mathrm{i}{\mathrm{\Delta }}_x\cos k_x+2t_y\sin k_y\\ & +2t_{xy}\sin (k_x+k_y).\end{array}\end{array}

Moreover, the dynamical matrix {\mathcal{M}}_{\text{BdG}}(\mathit{k}) can be unitarily transformed into a block-diagonal form with

(C6)\begin{array}{c}{U}_{\mathit{k}}{\mathcal{M}}_{\text{BdG}}(\mathit{k})U_{\mathit{k}}^{†}=\left(\begin{array}{cc}{E}_{+}(\mathit{k})& 0\\ 0& E_{-}(\mathit{k})\end{array}\right),\end{array}

where U_{\mathit{k}}=({\sigma }_0-\mathrm{i}{\sigma }_x)/\sqrt{2}.

In real space, the dynamical matrix {\mathcal{M}}_{\text{BdG}} can likewise be expressed in a block-diagonal form

(C7)\begin{array}{c}U{\mathcal{M}}_{\text{BdG}}{U}^{†}=\left(\begin{array}{cc}{\mathcal{M}}_p& \mathcal{O}\\ \mathcal{O}& {\mathcal{M}}_m\end{array}\right),\end{array}

where U=({\tau }_0-\mathrm{i}{\tau }_x)/\sqrt{2} with U^{†}=U^{-1}, {\mathcal{M}}_p and {\mathcal{M}}_m denote the decoupled matrices in real space.

We now apply the same unitary transformation to the dynamical matrix \mathcal{M} in the presence of impurities, yielding \overline{\mathcal{M}}=U\mathcal{M}{U}^{†}. Under this transformation, the Green’s function \mathcal{G} in Eq. (C3) becomes

(C8)\begin{array}{c}\overline{\mathcal{G}}=(E{\tau }_0-\overline{\mathcal{M}}{)}^{-1}=({\tau }_0-{\overline{\mathcal{G}}}_0{\overline{\mathcal{V}}}_z{)}^{-1}{\overline{\mathcal{G}}}_0,\end{array}

where {\overline{\mathcal{G}}}_0=(E{\tau }_0-U{\mathcal{M}}_{\text{BdG}}{U}^{†}{)}^{-1}, and {\overline{\mathcal{V}}}_z=U{\mathcal{V}}_z{U}^{†}.

8.2 Spectral radius for on-site impurities

Here, we focus on the spectral radius in a bosonic quadratic Hermitian system perturbed by on-site impurities \hat{V}=V_1{\hat{a}}_{{\mathit{r}}_1}^{†}{\hat{a}}_{{\mathit{r}}_1}+V_2{\hat{a}}_{{\mathit{r}}_2}^{†}{\hat{a}}_{{\mathit{r}}_2} in Eq. (13).

The unitary transformation of the Green’s function in Eq. (C8) gives rise to a modified spectral stability criterion for infinitesimal impurities in Eq. (18), expressed via the reduced spectral radius with

(C9)\begin{array}{c}\rho ({\overline{\mathcal{G}}}_0{\overline{\mathcal{V}}}_z)\ll 1.\end{array}

We now turn to a detailed discussion of the calculation of this spectral radius \rho ({\overline{\mathcal{G}}}_0{\overline{\mathcal{V}}}_z).

In position representation, the response matrix {\overline{\mathcal{G}}}_0{\overline{\mathcal{V}}}_z can be written as

(C10)\begin{array}{c}\begin{array}{rl}{\overline{\mathcal{G}}}_0{\overline{\mathcal{V}}}_z=& \left(\begin{array}{cccccccccc}0& 0& 0& \cdots & 0& \mathrm{i}{V}_1{\overline{\mathcal{G}}}_0^p({\mathit{r}}_1,{\mathit{r}}_1;E)& \mathrm{i}{V}_2{\overline{\mathcal{G}}}_0^p({\mathit{r}}_1,{\mathit{r}}_2;E)& 0& \cdots & 0\\ 0& 0& 0& \cdots & 0& \mathrm{i}{V}_1{\overline{\mathcal{G}}}_0^p({\mathit{r}}_2,{\mathit{r}}_1;E)& \mathrm{i}{V}_2{\overline{\mathcal{G}}}_0^p({\mathit{r}}_2,{\mathit{r}}_2;E)& 0& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & 0& \mathrm{i}{V}_1{\overline{\mathcal{G}}}_0^p({\mathit{r}}_N,{\mathit{r}}_1;E)& \mathrm{i}{V}_2{\overline{\mathcal{G}}}_0^p({\mathit{r}}_N,{\mathit{r}}_2;E)& 0& \cdots & 0\\ -\mathrm{i}{V}_1{\overline{\mathcal{G}}}_0^m({\mathit{r}}_1,{\mathit{r}}_1;E)& -\mathrm{i}{V}_2{\overline{\mathcal{G}}}_0^m({\mathit{r}}_1,{\mathit{r}}_2;E)& 0& \cdots & 0& 0& 0& 0& \cdots & 0\\ -\mathrm{i}{V}_1{\overline{\mathcal{G}}}_0^m({\mathit{r}}_2,{\mathit{r}}_1;E)& -\mathrm{i}{V}_2{\overline{\mathcal{G}}}_0^m({\mathit{r}}_2,{\mathit{r}}_2;E)& 0& \cdots & 0& 0& 0& 0& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ -\mathrm{i}{V}_1{\overline{\mathcal{G}}}_0^m({\mathit{r}}_N,{\mathit{r}}_1;E)& -\mathrm{i}{V}_2{\overline{\mathcal{G}}}_0^m({\mathit{r}}_N,{\mathit{r}}_2;E)& 0& \cdots & 0& 0& 0& 0& \cdots & 0\end{array}\right)\\ =& \left(\begin{array}{cc}\mathcal{O}& \mathcal{B}\\ \mathcal{C}& \mathcal{O}\end{array}\right),\end{array}\end{array}

where the element of the unperturbed Green’s function is given by

(C11)\begin{array}{c}{\overline{\mathcal{G}}}_0^{p/m}({\mathit{r}}_i,{\mathit{r}}_j;E)=⟨{\mathit{r}}_i|(E-{\mathcal{M}}_{p/m}{)}^{-1}|{\mathit{r}}_j⟩.\end{array}

As discussed in detail in Appendix C.3, the eigenvalues of the reduced response matrix {\overline{\mathcal{G}}}_0{\overline{\mathcal{V}}}_z in Eq. (C10) are determined by the matrix product \mathcal{B}\mathcal{C}, where the response matrix possesses two (four) non-zero eigenvalues for single on-site impurity (double on-site impurities) [see Appendix C.3].

When a single impurity is introduced (V_1\ne 0,V_2=0), the matrix \mathcal{B}\mathcal{C} has only one non-zero eigenvalue

(C12)\begin{array}{c}\xi =V_1^2{\overline{\mathcal{G}}}_0^p({\mathit{r}}_1,{\mathit{r}}_1;E){\overline{\mathcal{G}}}_0^m({\mathit{r}}_1,{\mathit{r}}_1;E),\end{array}

while the remaining N-1 eigenvalues are zero. Consequently, the response matrix {\overline{\mathcal{G}}}_0{\overline{\mathcal{V}}}_z has two non-zero eigenvalues, \pm \sqrt{\xi }, and 2N-2 zeros. Then, the spectral radius is derived as

(C13)\begin{array}{c}\rho ({\overline{\mathcal{G}}}_0{\overline{\mathcal{V}}}_z)=|V_1\sqrt{{\overline{\mathcal{G}}}_0^p({\mathit{r}}_1,{\mathit{r}}_1;E){\overline{\mathcal{G}}}_0^m({\mathit{r}}_1,{\mathit{r}}_1;E)}|.\end{array}

Furthermore, when two impurities are introduced (V_1,V_2\ne 0), the matrix \mathcal{B}\mathcal{C} has non-zero entries solely in the upper-left 2\times 2 block with

(C14)\begin{array}{c}\mathcal{B}\mathcal{C}=\left(\begin{array}{ccccc}{\mathcal{B}\mathcal{C}}_{11}& {\mathcal{B}\mathcal{C}}_{12}& 0& \cdots & 0\\ {\mathcal{B}\mathcal{C}}_{21}& {\mathcal{B}\mathcal{C}}_{22}& 0& \cdots & 0\\ 0& 0& 0& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & 0\end{array}\right),\end{array}

where four non-zero elements of the matrix \mathcal{B}\mathcal{C} are written as

(C15)\begin{array}{c}\begin{array}{rl}{\mathcal{B}\mathcal{C}}_{11}=& V_1^2{\overline{\mathcal{G}}}_0^p({\mathit{r}}_1,{\mathit{r}}_1;E){\overline{\mathcal{G}}}_0^m({\mathit{r}}_1,{\mathit{r}}_1;E)\\ & +V_1{V}_2{\overline{\mathcal{G}}}_0^p({\mathit{r}}_1,{\mathit{r}}_2;E){\overline{\mathcal{G}}}_0^m({\mathit{r}}_2,{\mathit{r}}_1;E),\end{array}\end{array}

(C16)\begin{array}{c}\begin{array}{rl}{\mathcal{B}\mathcal{C}}_{12}=& V_1{V}_2{\overline{\mathcal{G}}}_0^p({\mathit{r}}_1,{\mathit{r}}_1;E){\overline{\mathcal{G}}}_0^m({\mathit{r}}_1,{\mathit{r}}_2;E)\\ & +V_2^2{\overline{\mathcal{G}}}_0^p({\mathit{r}}_1,{\mathit{r}}_2;E){\overline{\mathcal{G}}}_0^m({\mathit{r}}_2,{\mathit{r}}_2;E),\end{array}\end{array}

(C17)\begin{array}{c}\begin{array}{rl}{\mathcal{B}\mathcal{C}}_{21}=& V_1^2{\overline{\mathcal{G}}}_0^p({\mathit{r}}_2,{\mathit{r}}_1;E){\overline{\mathcal{G}}}_0^m({\mathit{r}}_1,{\mathit{r}}_1;E)\\ & +V_1{V}_2{\overline{\mathcal{G}}}_0^p({\mathit{r}}_2,{\mathit{r}}_2;E){\overline{\mathcal{G}}}_0^m({\mathit{r}}_2,{\mathit{r}}_1;E),\end{array}\end{array}

(C18)\begin{array}{c}\begin{array}{rl}{\mathcal{B}\mathcal{C}}_{22}=& V_1{V}_2{\overline{\mathcal{G}}}_0^p({\mathit{r}}_2,{\mathit{r}}_1;E){\overline{\mathcal{G}}}_0^m({\mathit{r}}_1,{\mathit{r}}_2;E)\\ & +V_2^2{\overline{\mathcal{G}}}_0^p({\mathit{r}}_2,{\mathit{r}}_2;E){\overline{\mathcal{G}}}_0^m({\mathit{r}}_2,{\mathit{r}}_2;E).\end{array}\end{array}

The matrix \mathcal{B}\mathcal{C} has two non-zero eigenvalues, which are given by

(C19)\begin{array}{c}{\xi }_{\pm }={\displaystyle \frac{1}{2}}\left({\mathcal{B}\mathcal{C}}_{11}+{\mathcal{B}\mathcal{C}}_{22}\pm \sqrt{({\mathcal{B}\mathcal{C}}_{11}-{\mathcal{B}\mathcal{C}}_{22}{)}^2+4{\mathcal{B}\mathcal{C}}_{12}{\mathcal{B}\mathcal{C}}_{21}}\right).\end{array}

Therefore, the spectral radius of the response matrix for two on-site impurities is

(C20)\begin{array}{c}\rho ({\overline{\mathcal{G}}}_0{\overline{\mathcal{V}}}_z)=max(|\sqrt{{\xi }_{\pm }}|).\end{array}

8.3 Details on eigenvalues of response matrix

Here, we provide a detailed discussion of how the eigenvalues of the response matrix {\overline{\mathcal{G}}}_0{\overline{\mathcal{V}}}_z in Eq. (C10) are obtained, which is written in the off-diagonal block-matrix form as

(C21)\begin{array}{c}{\overline{\mathcal{G}}}_0{\overline{\mathcal{V}}}_z=\left(\begin{array}{cc}\mathcal{O}& \mathcal{B}\\ \mathcal{C}& \mathcal{O}\end{array}\right),\end{array}

where \mathcal{O} denotes the zero matrix, while \mathcal{B} and \mathcal{C} represent the upper-right and lower-left subblocks of the response matrix, respectively.

The eigenvalue problem of the response matrix is written as

(C22)\begin{array}{c}\left(\begin{array}{cc}\mathcal{O}& \mathcal{B}\\ \mathcal{C}& \mathcal{O}\end{array}\right)\left(\begin{array}{c}\mathit{x}\\ \mathit{y}\end{array}\right)=\lambda \left(\begin{array}{c}\mathit{x}\\ \mathit{y}\end{array}\right),\end{array}

which yields the coupled equations

(C23)\begin{array}{c}\{\begin{array}{l}\mathcal{B}\mathit{y}=\lambda \mathit{x},\\ \mathcal{C}\mathit{x}=\lambda \mathit{y}.\end{array}\end{array}

By eliminating either \mathit{x} or \mathit{y} from these equations, we obtain reduced eigenvalue problems with

(C24)\begin{array}{c}\mathcal{B}\mathcal{C}\mathit{x}={\lambda }^2\mathit{x},\mathcal{C}\mathcal{B}\mathit{y}={\lambda }^2\mathit{y}.\end{array}

This relation shows that the nonzero eigenvalues \lambda of the original block matrix must satisfy {\lambda }^2=\xi, where \xi is a nonzero eigenvalue of the matrix product \mathcal{B}\mathcal{C} (or \mathcal{C}\mathcal{B}). As a result, for each nonzero \xi, the response operator has a pair of eigenvalues \lambda =\pm \sqrt{\xi }.

If either \mathcal{B}\mathcal{C} or \mathcal{C}\mathcal{B} possesses zero eigenvalues (i.e., \xi =0), these translate directly into zero eigenvalues of the original off-diagonal block matrix in Eq. (C21). The total number of such zero eigenvalues is determined by the nullity (dimension of the kernel) of the matrices \mathcal{B} and \mathcal{C}. Specifically, \dim (\ker (\mathcal{B}))+\dim (\ker (\mathcal{C})) gives the number of independent zero eigenvectors of the block matrix. For the matrix \mathcal{B} and \mathcal{C} defined in Eq. (C10), it follows that \dim (\ker (\mathcal{B}))+\dim (\ker (\mathcal{C}))=2N-2 (\dim (\ker (\mathcal{B}))+\dim (\ker (\mathcal{C}))=2N-4) for one impurity (two impurities). Therefore, the block matrix has exactly 2N-2 (2N-4) zero eigenvalues for one impurity (two impurities). These zero modes arise from the localized nature of the on-site impurities, which affect only a restricted subset of the response matrix. For spectral analysis, however, the focus is on the nonzero eigenvalues. The physically relevant part of the spectrum is governed by the matrix \mathcal{B}\mathcal{C}, and all spectral radius calculations are therefore performed with respect to this matrix.

9 Appendix D: Analytical approximation of the propagator

We apply the GBZ Green-function method of Ref. [100] to the particle and hole blocks {\mathcal{M}}_{p/m}. Under OBC along x and PBC along y, we introduce {\beta }_x={\mathrm{e}}^{\mathrm{i}{k}_x+{\mu }_x} and obtain

(D1)\begin{array}{c}\begin{array}{rl}& {\overline{\mathcal{G}}}_0^{p/m}({\mathit{r}}_1,{\mathit{r}}_2;E)\\ & =⟨{\mathit{r}}_1|(E-{\mathcal{M}}_{p/m}{)}^{-1}|{\mathit{r}}_2⟩\\ & ={\int }_0^{2\pi }\mathrm{d}{k}_y{\oint }_{{\beta }_x\in \text{GBZ}}\frac{\text{d}{\beta }_x}{{\beta }_x}\frac{{\mathrm{e}}^{\mathrm{i}{k}_y(y_1-y_2)}{\beta }_x^{x_1-x_2}}{E-E_{\pm }({\beta }_x,k_y)}.\end{array}\end{array}

Here, the constant prefactor is omitted and consistently neglected in the following derivations, and E_{\pm }({\beta }_x,k_y) reads

(D2)\begin{array}{c}\begin{array}{rl}{E}_{\pm }({\beta }_x,k_y)=& \pm \mathrm{i}{\mathrm{\Delta }}_x\left({\beta }_x+{\beta }_x^{-1}\right)+\mathrm{i}{t}_y\left({\mathrm{e}}^{\mathrm{i}{k}_y}-{\mathrm{e}}^{-\mathrm{i}{k}_y}\right)\\ & \pm 2\mathrm{i}{\mathrm{\Delta }}_0+\mathrm{i}{t}_{xy}\left({\beta }_x{\mathrm{e}}^{\mathrm{i}{k}_y}-{\beta }_x^{-1}{\mathrm{e}}^{-\mathrm{i}{k}_y}\right)\\ & +\mathrm{i}{t}_y\left({\mathrm{e}}^{\mathrm{i}{k}_y}-{\mathrm{e}}^{-\mathrm{i}{k}_y}\right).\end{array}\end{array}

For the contour integral in Eq. (D1), we consider the new eigenvalue E of the dynamical matrix \mathcal{M}, in the presence of on-site impurities, outside the eigenspectrum (indicated by {\sigma }_{\mathrm{c}\mathrm{y}\mathrm{l}}) of {\mathcal{M}}_{p/m} (i.e., E\notin {\sigma }_{\mathrm{c}\mathrm{y}\mathrm{l}}) in the cylindrical geometry. To carry out this contour integral, we solve characteristic equations E-E_{\pm }({\beta }_x,k_y)=0 for the propagators {\overline{\mathcal{G}}}_0^{p/m}({\mathit{r}}_1,{\mathit{r}}_2;E) in the subspace defined in Eq. (C7). The corresponding roots of the characteristic equations, {\beta }_{1,2}^{\pm }(k_y), are ordered such that |{\beta }_2^{\pm }(k_y)|>|{\beta }_1^{\pm }(k_y)|. Then, the equation (D2) becomes

(D3)\begin{array}{c}E-E_{\pm }({\beta }_x,k_y)={\displaystyle \frac{\left[{\beta }_x-{\beta }_1^{\pm }(k_y)\right]\left[{\beta }_x-{\beta }_2^{\pm }(k_y)\right]}{{\beta }_x}}.\end{array}

Note that the GBZ is determined by standing-wave condition |{\beta }_1^{\pm }(k_y)|=|{\beta }_2^{\pm }(k_y)|.

Since {\mathcal{M}}_{p/m}^{\mathrm{T}}\ne {\mathcal{M}}_{p/m}, the two directional propagators are unequal. For x_2>x_1,

(D4)\begin{array}{c}\begin{array}{rl}& {\overline{\mathcal{G}}}_0^{p/m}({\mathit{r}}_2,{\mathit{r}}_1;E)={\int }_0^{2\pi }\mathrm{d}{k}_y{\mathrm{e}}^{\mathrm{i}{k}_y(y_2-y_1)}\\ & \times {\oint }_{{\beta }_x\in \text{GBZ}}\text{d}{\beta }_x\frac{{\beta }_x^{x_2-x_1}}{\left[{\beta }_x-{\beta }_1^{\pm }(k_y)\right]\left[{\beta }_x-{\beta }_2^{\pm }(k_y)\right]},\end{array}\end{array}

and

(D5)\begin{array}{c}\begin{array}{rl}& {\overline{\mathcal{G}}}_0^{p/m}({\mathit{r}}_1,{\mathit{r}}_2;E)={\int }_0^{2\pi }\mathrm{d}{k}_y{\mathrm{e}}^{\mathrm{i}{k}_y(y_1-y_2)}\\ & \times {\oint }_{{\beta }_x\in \text{GBZ}}\text{d}{\beta }_x\frac{1}{{\beta }_x^{x_2-x_1}\left[{\beta }_x-{\beta }_1^{\pm }(k_y)\right]\left[{\beta }_x-{\beta }_2^{\pm }(k_y)\right]}.\end{array}\end{array}

By applying Vieta’s formulas to the quadratic polynomial in Eq. (D3), we can verify that {\beta }_1^{\pm }(k_y) lies inside the GBZ contour, whereas {\beta }_2^{\pm }(k_y) lies outside it, as shown in Fig. A3. When x_2>x_1, the integrand of {\overline{\mathcal{G}}}_0^{p/m}({\mathit{r}}_2,{\mathit{r}}_1;E) therefore possesses one pole {\beta }_1^{\pm } inside the contour, and {\overline{\mathcal{G}}}_0^{p/m}({\mathit{r}}_1,{\mathit{r}}_2;E) has two poles {\beta }_1^{\pm } [see Fig. A3(a)], and 0 inside the contour [see Fig. A3(b)]. Applying the residue theorem along the GBZ contour, we obtain

(D6)\begin{array}{c}{\overline{\mathcal{G}}}_0^{p/m}({\mathit{r}}_2,{\mathit{r}}_1;E)={\displaystyle {\int }_0^{2\pi }\mathrm{d}{k}_y{\displaystyle \frac{[{\beta }_1^{\pm }(k_y){]}^{x_2-x_1}{\mathrm{e}}^{\mathrm{i}{k}_y(y_2-y_1)}}{{\beta }_1^{\pm }(k_y)-{\beta }_2^{\pm }(k_y)}},}\end{array}

and

(D7)\begin{array}{c}{\overline{\mathcal{G}}}_0^{p/m}({\mathit{r}}_1,{\mathit{r}}_2;E)={\displaystyle {\int }_0^{2\pi }\mathrm{d}{k}_y{\displaystyle \frac{[{\beta }_2^{\pm }(k_y){]}^{x_1-x_2}{\mathrm{e}}^{\mathrm{i}{k}_y(y_1-y_2)}}{{\beta }_1^{\pm }(k_y)-{\beta }_2^{\pm }(k_y)}}.}\end{array}

For L=x_2-x_1 and {\delta }_y=y_2-y_1, the nonexponential denominators can be omitted when estimating the leading large-L behavior:

(D8)\begin{array}{c}{\overline{\mathcal{G}}}_0^{p/m}({\mathit{r}}_2,{\mathit{r}}_1;E)\simeq {\displaystyle {\int }_0^{2\pi }\mathrm{d}{k}_y[{\beta }_1^{\pm }(k_y){]}^{\mathrm{L}}{\mathrm{e}}^{\mathrm{i}{k}_y{\delta }_y},}\end{array}

(D9)\begin{array}{c}{\overline{\mathcal{G}}}_0^{p/m}({\mathit{r}}_1,{\mathit{r}}_2;E)\simeq {\displaystyle {\int }_0^{2\pi }\mathrm{d}{k}_y[{\displaystyle \frac{1}{{\beta }_2^{\pm }(k_y)}}{]}^{\mathrm{L}}{\mathrm{e}}^{-\mathrm{i}{k}_y{\delta }_y}.}\end{array}

Here, \asymp denotes equality at leading exponential order. These propagators are subsequently inserted into \rho ({\overline{\mathcal{G}}}_0{\overline{\mathcal{V}}}_z) to determine the two-impurity spectral response.

For the new eigenvalue E with positive imaginary part, for instance, the eigenvalue E indicated by the orange circle in Fig. 3(h), the integrals in Eqs. (D8) and (D9) can be estimated as

(D10)\begin{array}{c}\{\begin{array}{l}{\overline{\mathcal{G}}}_0^p({\mathit{r}}_2,{\mathit{r}}_1;E)\sim {\mathrm{e}}^{{\mu }_{\text{max,1}}^{+}L},\\ {\overline{\mathcal{G}}}_0^m({\mathit{r}}_2,{\mathit{r}}_1;E)\sim {\mathrm{e}}^{{\mu }_{\text{max,1}}^{-}L}\sim 0,\\ {\overline{\mathcal{G}}}_0^p({\mathit{r}}_1,{\mathit{r}}_2;E)\sim {\mathrm{e}}^{-{\mu }_{\text{min,2}}^{+}L},\\ {\overline{\mathcal{G}}}_0^m({\mathit{r}}_1,{\mathit{r}}_2;E)\sim {\mathrm{e}}^{-{\mu }_{\text{min,2}}^{-}L}\sim 0,\end{array}\end{array}

with

(D11)\begin{array}{c}{\mu }_{\text{max,1}}^{\pm }=\underset{k_y\in [0,2\pi ]}{max}\ln |{\beta }_1^{\pm }(k_y)|,\end{array}

and

(D12)\begin{array}{c}{\mu }_{\text{min,2}}^{\pm }=\underset{k_y\in [0,2\pi ]}{min}\ln |{\beta }_2^{\pm }(k_y)|.\end{array}

Here, {\mu }_{\text{max},1}^{+}>0, , , and {\mu }_{\text{min},2}^{-}>0 are typically satisfied in the parameter regime of interest [see red and blue dots in Figs. 5(a1) and (a3)], thereby facilitating the demonstration of the ultra spectral sensitivity.

Therefore, for double on-site impurities, the spectral radius associated with the induced new eigenvalue E with positive imaginary part can be approximated by

(D13)\begin{array}{c}\rho ({\overline{\mathcal{G}}}_0{\overline{\mathcal{V}}}_z)\sim V{\mathrm{e}}^{\frac{{\mu }_{\text{max,1}}^{+}L}{4}}{\mathrm{e}}^{-\frac{{\mu }_{\text{min,2}}^{+}L}{4}}.\end{array}

Here, we consider the case where the two on-site impurities have identical on-site potentials with V_1=V_2=V.

For the new eigenvalue E with negative imaginary part, for instance, the eigenvalue E marked by the red circle in Fig. 3(h), the integrals in Eqs. (D8), and (D9) can be estimated as

(D14)\begin{array}{c}\{\begin{array}{l}{\overline{\mathcal{G}}}_0^p({\mathit{r}}_2,{\mathit{r}}_1;E)\sim {\mathrm{e}}^{{\mu }_{\text{max,1}}^{+}L}\sim 0,\\ {\overline{\mathcal{G}}}_0^m({\mathit{r}}_2,{\mathit{r}}_1;E)\sim {\mathrm{e}}^{{\mu }_{\text{max,1}}^{-}L},\\ {\overline{\mathcal{G}}}_0^p({\mathit{r}}_1,{\mathit{r}}_2;E)\sim {\mathrm{e}}^{-{\mu }_{\text{min,2}}^{+}L}\sim 0,\\ {\overline{\mathcal{G}}}_0^m({\mathit{r}}_1,{\mathit{r}}_2;E)\sim {\mathrm{e}}^{-{\mu }_{\text{min,2}}^{-}L},\end{array}\end{array}

where , {\mu }_{\text{max},1}^{-}>0, {\mu }_{\text{min},2}^{+}>0, and are typically satisfied in the parameter regime of interest [see red and blue dots in Figs. 5(b1) and (b3)]. The corresponding spectral radius associated with the induced new eigenvalue E with negative imaginary part is

(D15)\begin{array}{c}\rho ({\overline{\mathcal{G}}}_0{\overline{\mathcal{V}}}_z)\sim V{\mathrm{e}}^{\frac{{\mu }_{\text{max,1}}^{-}L}{4}}{\mathrm{e}}^{-\frac{{\mu }_{\text{min,2}}^{-}L}{4}},\end{array}

where V_1=V_2=V is assumed.

References

[1] Y. Ashida , Z. Gong , and M. Ueda , Non-Hermitian physics, Adv. Phys. 69(3), 249 (2020)

[2] F. Minganti , A. Miranowicz , R. W. Chhajlany , I. I. Arkhipov , and F. Nori , Hybrid-Liouvillian formalism connecting exceptional points of non-Hermitian Hamiltonians and Liouvillians via postselection of quantum trajectories, Phys. Rev. A 101(6), 062112 (2020)

[3] E. J. Bergholtz , J. C. Budich , and F. K. Kunst , Exceptional topology of non-Hermitian systems, Rev. Mod. Phys. 93(1), 015005 (2021)

[4] T. Micallo , C. Lehmann , and J. C. Budich , Correlation-induced sensitivity and non-Hermitian skin effect of quasiparticles, Phys. Rev. Res. 5(4), 043105 (2023)

[5] C. Leefmans , A. Dutt , J. Williams , L. Yuan , M. Parto , F. Nori , S. Fan , and A. Marandi , Topological dissipation in a time-multiplexed photonic resonator network, Nat. Phys. 18(4), 442 (2022)

[6] Y. Qin , K. Zhang , and L. Li , Geometry-dependent skin effect and anisotropic Bloch oscillations in a non-Hermitian optical lattice, Phys. Rev. A 109(2), 023317 (2024)

[7] H. Xiao and Q. B. Zeng , Coexistence of non-Hermitian skin effect and extended states in one-dimensional nonreciprocal lattices, Phys. Rev. B 110(2), 024205 (2024)

[8] M. Reisenbauer , H. Rudolph , L. Egyed , K. Hornberger , A. V. Zasedatelev , M. Abuzarli , B. A. Stickler , and U. Delić , Non-Hermitian dynamics and non-reciprocity of optically coupled nanoparticles, Nat. Phys. 20(10), 1629 (2024)

[9] K. Ochkan , R. Chaturvedi , V. Könye , L. Veyrat , R. Giraud , D. Mailly , A. Cavanna , U. Gennser , E. M. Hankiewicz , B. Büchner , J. van den Brink , J. Dufouleur , and I. C. Fulga , Non-Hermitian topology in a multi-terminal quantum Hall device, Nat. Phys. 20(3), 395 (2024)

[10] Y. M. Hu , H. Y. Wang , Z. Wang , and F. Song , Geometric origin of non-Bloch PT symmetry breaking, Phys. Rev. Lett. 132(5), 050402 (2024)

[11] S. Longhi , Erratic non-Hermitian skin localization, Phys. Rev. Lett. 134(19), 196302 (2025)

[12] K. Zhang , Z. Yang , and K. Sun , Edge theory of non-Hermitian skin modes in higher dimensions, Phys. Rev. B 109(16), 165127 (2024)

[13] B. Li , C. Chen , and Z. Wang , Universal non-Hermitian transport in disordered systems, Phys. Rev. Lett. 135(3), 033802 (2025)

[14] S. X. Wang and Z. Yan , Enhanced sensitivity in non-Hermitian systems at infernal points, Phys. Rev. Res. 7(2), L022037 (2025)

[15] J. Gliozzi , G. De Tomasi , and T. L. Hughes , Many-body non-Hermitian skin effect for multipoles, Phys. Rev. Lett. 133(13), 136503 (2024)

[16] T. Yoshida , S. B. Zhang , T. Neupert , and N. Kawakami , Non-Hermitian Mott skin effect, Phys. Rev. Lett. 133(7), 076502 (2024)

[17] A. Yang , Z. Fang , K. Zhang , and C. Fang , Tailoring bound state geometry in high-dimensional non-Hermitian systems, Commun. Phys. 8(1), 124 (2025)

[18] P. Molignini , O. Arandes , and E. J. Bergholtz , Anomalous skin effects in disordered systems with a single non-Hermitian impurity, Phys. Rev. Res. 5(3), 033058 (2023)

[19] H. Jiang and C. H. Lee , Dimensional transmutation from non-Hermiticity, Phys. Rev. Lett. 131(7), 076401 (2023)

[20] C. A. Li , B. Trauzettel , T. Neupert , and S. B. Zhang , Enhancement of second-order non-Hermitian skin effect by magnetic fields, Phys. Rev. Lett. 131(11), 116601 (2023)

[21] Z. Gong , M. Bello , D. Malz , and F. K. Kunst , Anomalous behaviors of quantum emitters in non-Hermitian baths, Phys. Rev. Lett. 129(22), 223601 (2022)

[22] Z. F. Cai , Y. C. Wang , Y. R. Zhang , T. Liu , and F. Nori , Versatile control of nonlinear topological states in non-Hermitian systems, Commun. Phys. 8(1), 360 (2025)

[23] L. H. Mo , Z. Xiao , R. Moessner , and H. Zhao , Non-Hermitian delocalization in one dimension via emergent compactness, Phys. Rev. B 111(23), 235412 (2025)

[24] D. C. Ohnmacht , V. Wilhelm , H. Weisbrich , and W. Belzig , Non-Hermitian topology in multiterminal superconducting junctions, Phys. Rev. Lett. 134(15), 156601 (2025)

[25] J. Shi and A. N. Poddubny , Chiral dissociation of bound photon pairs for a non-Hermitian skin effect, Phys. Rev. Lett. 134(23), 233602 (2025)

[26] A. Hashemi , E. L. Pereira , H. Li , J. L. Lado , and A. Blanco-Redondo , Observation of non-Hermitian topology from optical loss modulation, Nat. Mater. 24(9), 1393 (2025)

[27] K. Sun and H. Hu , Lyapunov formulation of band theory for disordered non-Hermitian systems, (2025)

[28] D. Nakamura , K. Shiozaki , K. Shimomura , M. Sato , and K. Kawabata , Non-Hermitian origin of detachable boundary states in topological insulators, Phys. Rev. Lett. 135(9), 096601 (2025)

[29] K. Zhou , Z. Yang , B. Zeng , and Y. Hu , Critical non-Hermitian edge modes, (2025)

[30] B. B. Wang , Z. Cheng , H. Y. Zou , Y. Ge , K. Q. Zhao , Q. R. Si , S. Q. Yuan , H. X. Sun , H. Xue , and B. Zhang , Observation of disorder-induced boundary localization, Proc. Natl. Acad. Sci. USA 122(19), e2422154122 (2025)

[31] J. P. Esparza and V. Juričić , Exceptional magic angles in non-Hermitian twisted bilayer graphene, Phys. Rev. Lett. 134(22), 226602 (2025)

[32] T. R. Liu , T. Liu , and M. Xiao , Anomalous non- Hermitian skin effect of chiral boundary states, Phys. Rev. B 112(8), L081112 (2025)

[33] Z. Ren , D. Liu , E. Zhao , C. He , K. K. Pak , J. Li , and G. B. Jo , Chiral control of quantum states in non- Hermitian spin–orbit-coupled fermions, Nat. Phys. 18(4), 385 (2022)

[34] M. Naghiloo , M. Abbasi , Y. N. Joglekar , and K. W. Murch , Quantum state tomography across the exceptional point in a single dissipative qubit, Nat. Phys. 15(12), 1232 (2019)

[35] R. Shen , T. Chen , B. Yang , and C. H. Lee , Observation of the non-Hermitian skin effect and Fermi skin on a digital quantum computer, Nat. Commun. 16, 1340 (2025)

[36] S. Imhof , C. Berger , F. Bayer , J. Brehm , L. W. Molenkamp , T. Kiessling , F. Schindler , C. H. Lee , M. Greiter , T. Neupert , and R. Thomale , Topolectrical-circuit realization of topological corner modes, Nat. Phys. 14(9), 925 (2018)

[37] R. El-Ganainy , K. G. Makris , M. Khajavikhan , Z. H. Musslimani , S. Rotter , and D. N. Christodoulides , Non-Hermitian physics and PT symmetry, Nat. Phys. 14(1), 11 (2018)

[38] J. Wu , Y. Hu , Z. He , K. Deng , X. Huang , M. Ke , W. Deng , J. Lu , and Z. Liu , Hybrid-order skin effect from loss-induced nonreciprocity, Phys. Rev. Lett. 134(17), 176601 (2025)

[39] S. Yao and Z. Wang , Edge states and topological invariants of non-Hermitian systems, Phys. Rev. Lett. 121(8), 086803 (2018)

[40] K. Yokomizo and S. Murakami , Non-Bloch band theory of non-Hermitian systems, Phys. Rev. Lett. 123(6), 066404 (2019)

[41] K. Zhang , Z. Yang , and C. Fang , Correspondence between winding numbers and skin modes in non-Hermitian systems, Phys. Rev. Lett. 125(12), 126402 (2020)

[42] Y. Xiong , Z. Y. Xing , and H. Hu , Non-Hermitian skin effect in arbitrary dimensions: non-Bloch band theory and classification, (2024)

[43] T. Liu , Y. R. Zhang , Q. Ai , Z. Gong , K. Kawabata , M. Ueda , and F. Nori , Second-order topological phases in non-Hermitian systems, Phys. Rev. Lett. 122(7), 076801 (2019)

[44] F. K. Kunst , E. Edvardsson , J. C. Budich , and E. J. Bergholtz , Biorthogonal bulk-boundary corre- spondence in non-Hermitian systems, Phys. Rev. Lett. 121(2), 026808 (2018)

[45] W. Z. Ling , Z. F. Cai , and T. Liu , Interaction-induced second-order skin effect, Phys. Rev. B 111(20), 205418 (2025)

[46] Y. Li , Z. F. Cai , T. Liu , and F. Nori , Dissipation and interaction-controlled non-Hermitian skin effects, (2024)

[47] Z. Gong , Y. Ashida , K. Kawabata , K. Takasan , S. Higashikawa , and M. Ueda , Topological phases of non-Hermitian systems, Phys. Rev. X 8(3), 031079 (2018)

[48] S. Yao , F. Song , and Z. Wang , Non-Hermitian Chern bands, Phys. Rev. Lett. 121(13), 136802 (2018)

[49] F. Song , S. Yao , and Z. Wang , Non-Hermitian skin effect and chiral damping in open quantum systems, Phys. Rev. Lett. 123(17), 170401 (2019)

[50] J. Y. Lee , J. Ahn , H. Zhou , and A. Vishwanath , Topological correspondence between Hermitian and non-Hermitian systems: Anomalous dynamics, Phys. Rev. Lett. 123(20), 206404 (2019)

[51] K. Kawabata , T. Bessho , and M. Sato , Classification of exceptional points and non-Hermitian topological semimetals, Phys. Rev. Lett. 123(6), 066405 (2019)

[52] J. Zhang , B. Peng , Ş. K. Özdemir , K. Pichler , D. O. Krimer , G. Zhao , F. Nori , Y. Liu , S. Rotter , and L. Yang , A phonon laser operating at an exceptional point, Nat. Photonics 12(8), 479 (2018)

[53] Z. Y. Ge , Y. R. Zhang , T. Liu , S. W. Li , H. Fan , and F. Nori , Topological band theory for non-Hermitian systems from the Dirac equation, Phys. Rev. B 100, 054105 (2019)

[54] F. Minganti , A. Miranowicz , R. W. Chhajlany , and F. Nori , Quantum exceptional points of non-Hermitian Hamiltonians and Liouvillians: The effects of quantum jumps, Phys. Rev. A 100(6), 062131 (2019)

[55] H. Zhao , X. Qiao , T. Wu , B. Midya , S. Longhi , and L. Feng , Non-Hermitian topological light steering, Science 365(6458), 1163 (2019)

[56] K. Kawabata , K. Shiozaki , M. Ueda , and M. Sato , Symmetry and topology in non-Hermitian physics, Phys. Rev. X 9(4), 041015 (2019)

[57] D. S. Borgnia , A. J. Kruchkov , and R. J. Slager , Non-Hermitian boundary modes and topology, Phys. Rev. Lett. 124(5), 056802 (2020)

[58] V. M. Martinez Alvarez , J. E. Barrios Vargas , and L. E. F. Foa Torres , Non-Hermitian robust edge states in one dimension: Anomalous localization and eigenspace condensation at exceptional points, Phys. Rev. B 97(12), 121401 (2018)

[59] Y. T. Zhang , S. Jiang , Q. Li , and Q. F. Sun , An analytical solution for quantum scattering through a PT-symmetric delta potential, Front. Phys. (Beijing) 16(4), 43503 (2021)

[60] Y.-W. Liang , Z.-F. Cai , T. Liu , and X. Wei , Quantum squeezing-induced critical non-Hermitian skin effects with z₂ skin modes, Front. Phys. (Beijing) 21, 095207 (2026)

[61] I. I. Arkhipov , A. Miranowicz , F. Minganti , and F. Nori , Liouvillian exceptional points of any order in dissipative linear bosonic systems: Coherence functions and switching between and anti-symmetries, Phys. Rev. A 102(3), 033715 (2020)

[62] N. Okuma , K. Kawabata , K. Shiozaki , and M. Sato , Topological origin of non-Hermitian skin effects, Phys. Rev. Lett. 124(8), 086801 (2020)

[63] T. Liu , J. J. He , Z. Yang , and F. Nori , Higher-order Weyl-exceptional-ring semimetals, Phys. Rev. Lett. 127(19), 196801 (2021)

[64] H. Zhang , Z. Guo , Y. Li , Y. Yang , Y. Chen , and H. Chen , A universal non-Hermitian platform for bound state in the continuum enhanced wireless power transfer, Front. Phys. (Beijing) 19(4), 43209 (2024)

[65] Y. Li , C. Liang , C. Wang , C. Lu , and Y. C. Liu , Gain-loss-induced hybrid skin-topological effect, Phys. Rev. Lett. 128(22), 223903 (2022)

[66] Y. Wang , J. Lin , and P. Xu , Transmission-reflection decoupling of non-Hermitian photonic doping epsilon-near-zero media, Front. Phys. (Beijing) 19(3), 33206 (2023)

[67] G. Sun , J. C. Tang , and S. P. Kou , Biorthogonal quantum criticality in non-Hermitian many-body systems, Front. Phys. (Beijing) 17(3), 33502 (2021)

[68] Y. P. Wu , G. Q. Zhang , C. X. Zhang , J. Xu , and D. W. Zhang , Interplay of nonreciprocity and nonlinearity on mean-field energy and dynamics of a Bose-Einstein condensate in a double-well potential, Front. Phys. (Beijing) 17(4), 42503 (2021)

[69] C. Chen , Y. Liu , L. Zhao , X. Hu , and Y. Fu , Asymmetric nonlinear-mode-conversion in an optical waveguide with PT symmetry, Front. Phys. (Beijing) 17(5), 52504 (2022)

[70] K. Li and Y. Xu , Non-Hermitian absorption spectroscopy, Phys. Rev. Lett. 129(9), 093001 (2022)

[71] N. Zeng , T. Liu , K. Xia , Y.-R. Zhang , and F. Nori , Non-Hermitian sensing from the perspective of post-selected measurements, (2025)

[72] Z.-F. Cai , X. Wang , Z.-X. Liang , T. Liu , and F. Nori , Chiral-extended photon-emitter dressed states in non-Hermitian topological baths, Phys. Rev. A 111, L061701 (2025)

[73] J. R. Li , C. Jiang , H. Su , D. Qi , L. L. Zhang , and W. J. Gong , Parity-dependent skin effects and topological properties in the multilayer nonreciprocal Su−Schrieffer−Heeger structures, Front. Phys. (Beijing) 19(3), 33204 (2023)

[74] R. Lin , T. Tai , L. Li , and C. H. Lee , Topological non-Hermitian skin effect, Front. Phys. (Beijing) 18(5), 53605 (2023)

[75] K. Zhang , C. Fang , and Z. Yang , Dynamical degeneracy splitting and directional invisibility in non-Hermitian systems, Phys. Rev. Lett. 131(3), 036402 (2023)

[76] Y. Y. Zou , Y. Zhou , L. M. Chen , and P. Ye , Detecting bulk and edge exceptional points in non-Hermitian systems through generalized Petermann factors, Front. Phys. (Beijing) 19(2), 23201 (2023)

[77] H. Wang , J. Liu , T. Liu , and W. B. Ju , Observation of impurity-induced scale-free localization in a disordered non-Hermitian electrical circuit, Front. Phys. (Beijing) 20(1), 014203 (2025)

[78] X. Li , J. Liu , and T. Liu , Localization-delocalization transitions in non-Hermitian Aharonov−Bohm cages, Front. Phys. (Beijing) 19(3), 33211 (2024)

[79] W. W. Jin , J. Liu , X. Wang , Y. R. Zhang , X. Huang , X. Wei , W. Ju , Z. Yang , T. Liu , and F. Nori , Anderson delocalization in strongly coupled disordered non-Hermitian chains, Phys. Rev. Lett. 135(7), 076602 (2025)

[80] M. Parto , C. Leefmans , J. Williams , F. Nori , and A. Marandi , Non-Abelian effects in dissipative photonic topological lattices, Nat. Commun. 14(1), 1440 (2023)

[81] H. Y. Wang , F. Song , and Z. Wang , Amoeba formulation of non-Bloch band theory in arbitrary dimensions, Phys. Rev. X 14(2), 021011 (2024)

[82] C. X. Guo , L. Su , Y. Wang , L. Li , J. Wang , X. Ruan , Y. Du , D. Zheng , S. Chen , and H. Hu , Scale-tailored localization and its observation in non-Hermitian electrical circuits, Nat. Commun. 15(1), 9120 (2024)

[83] H. Q. Liang , Z. Ou , L. Li , and G. F. Xu , Intrinsic perturbation induced anomalous higher-order boundary states in non-Hermitian systems, Phys. Rev. B 111(24), L241112 (2025)

[84] L. Li , C. H. Lee , S. Mu , and J. Gong , Critical non-Hermitian skin effect, Nat. Commun. 11(1), 5491 (2020)

[85] B. Li , H. R. Wang , F. Song , and Z. Wang , Scale-free localization and symmetry breaking from local non-PT hermiticity, Phys. Rev. B 108(16), L161409 (2023)

[86] Z. F. Cai , T. Liu , and Z. Yang , Non-Hermitian skin effect in periodically driven dissipative ultracold atoms, Phys. Rev. A 109(6), 063329 (2024)

[87] K. Kawabata , T. Numasawa , and S. Ryu , Entanglement phase transition induced by the non-Hermitian skin effect, Phys. Rev. X 13(2), 021007 (2023)

[88] W. T. Xue , Y. M. Hu , F. Song , and Z. Wang , Non-Hermitian edge burst, Phys. Rev. Lett. 128(12), 120401 (2022)

[89] W. T. Xue , F. Song , Y. M. Hu , and Z. Wang , Non-Bloch edge dynamics of non-Hermitian lattices, (2025)

[90] M. Yang and C. H. Lee , Percolation-induced PT symmetry breaking, Phys. Rev. Lett. 133, 136602 (2024)

[91] C. Yuce and H. Ramezani , Non-Hermitian edge burst without skin localization, Phys. Rev. B 107(14), L140302 (2023)

[92] L. Xiao , W. T. Xue , F. Song , Y. M. Hu , W. Yi , Z. Wang , and P. Xue , Observation of non-Hermitian edge burst in quantum dynamics, Phys. Rev. Lett. 133(7), 070801 (2024)

[93] T. H. Yang and C. Fang , Real-time edge dynamics of non-Hermitian lattices, Phys. Rev. Lett. 135(18), 186401 (2025)

[94] K. Zhang , Z. Yang , and C. Fang , Universal non-Hermitian skin effect in two and higher dimensions, Nat. Commun. 13(1), 2496 (2022)

[95] W. Wang , M. Hu , X. Wang , G. Ma , and K. Ding , Experimental realization of geometry-dependent skin Phys. effect in a reciprocal two-dimensional lattice, Phys. Rev. Lett. 131(20), 207201 (2023)

[96] Q. Zhou , J. Wu , Z. Pu , J. Lu , X. Huang , W. Deng , M. Ke , and Z. Liu , Observation of geometry- dependent skin effect in non-Hermitian phononic crystals with exceptional points, Nat. Commun. 14(1), 4569 (2023)

[97] T. Wan , K. Zhang , J. Li , Z. Yang , and Z. Yang , Observation of the geometry-dependent skin effect and dynamical degeneracy splitting, Sci. Bull. (Beijing) 68(20), 2330 (2023)

[98] K. Zhang , C. Shu , and K. Sun , Algebraic non- Hermitian skin effect and generalized Fermi surface formula in arbitrary dimensions, Phys. Rev. X 15(3), 031039 (2025)

[99] K. Ding , M. Li , and J. Lin , Algebraic skin effect in two-dimensional non-Hermitian metamaterials, (2025)

[100] C. Shu , K. Zhang , and K. Sun , Ultraspectral sensitivity and nonlocal bound states in algebraic non-Hermitian skin effect, Phys. Rev. B 112(23), 235152 (2025)

[101] Y. X. Wang and A. A. Clerk , Non-Hermitian dynamics without dissipation in quantum systems, Phys. Rev. A 99(6), 063834 (2019)

[102] A. McDonald , T. Pereg-Barnea , and A. A. Clerk , Phase-dependent chiral transport and effective non-Hermitian dynamics in a bosonic Kitaev-Majorana chain, Phys. Rev. X 8(4), 041031 (2018)

[103] V. P. Flynn , E. Cobanera , and L. Viola , Deconstructing effective non-Hermitian dynamics in quadratic bosonic Hamiltonians, New J. Phys. 22(8), 083004 (2020)

[104] Q. Wang , C. Zhu , Y. Wang , B. Zhang , and Y. D. Chong , Amplification of quantum signals by the non-Hermitian skin effect, Phys. Rev. B 106(2), 024301 (2022)

[105] N. Okuma , Boundary-dependent dynamical instability of bosonic Greens function: Dissipative Bogoliubov–de Gennes Hamiltonian and its application to non-Hermitian skin effect, Phys. Rev. B 105(22), 224301 (2022)

[106] X. W. Luo , C. Zhang , and S. Du , Quantum squeezing and sensing with pseudo-anti-parity−time symmetry, Phys. Rev. Lett. 128(17), 173602 (2022)

[107] C. Lv and Q. Zhou , Complexity geometry in Hermitian and non-Hermitian quantum dynamics, Phys. Rev. D 110(8), 084039 (2024)

[108] V. P. Flynn , E. Cobanera , and L. Viola , Topology by dissipation: Majorana bosons in metastable quadratic Markovian dynamics, Phys. Rev. Lett. 127(24), 245701 (2021)

[109] L. L. Wan and X. Y. Lü , Quantum-squeezing-induced point-gap topology and skin effect, Phys. Rev. Lett. 130(20), 203605 (2023)

[110] C. C. Wanjura , M. Brunelli , and A. Nunnenkamp , Topological framework for directional amplification in driven-dissipative cavity arrays, Nat. Commun. 11(1), 3149 (2020)

[111] S. Lieu , Topological symmetry classes for non- Hermitian models and connections to the bosonic Bogoliubov–de Gennes equation, Phys. Rev. B 98(11), 115135 (2018)

[112] Z. Yang , Non-perturbative breakdown of Bloch’s theorem and Hermitian skin effects, (2020)

[113] O. Arandes and E. J. Bergholtz , Quantum sensing with driven-dissipative Su−Schrieffer−Heeger lattices, Phys. Rev. Res. 7(1), 013309 (2025)

[114] K. B. Estake , T. R. Vishnu , and D. Roy , From chiral topological dynamics to chiral topological amplification: Real vs imaginary parameters in a Hermitian bosonic chain, (2025)

[115] K. Yokomizo and S. Murakami , Non-Bloch band theory in bosonic Bogoliubov–de Gennes systems, Phys. Rev. B 103(16), 165123 (2021)

[116] M. Bestler , A. Dikopoltsev , and O. Zilberberg , Non-Hermitian topology and skin modes in the continuum via parametric processes, (2025)

[117] J. J. Slim , C. C. Wanjura , M. Brunelli , J. del Pino , A. Nunnenkamp , and E. Verhagen , Optomechanical realization of the bosonic Kitaev chain, Nature 627(8005), 767 (2024)

[118] A. Pocklington , Y. X. Wang , and A. A. Clerk , Dissipative pairing interactions: Quantum instabilities, topological light, and volume-law entanglement, Phys. Rev. Lett. 130(12), 123602 (2023)

[119] B. Abdo , A. Kamal , and M. Devoret , Nondegenerate three-wave mixing with the Josephson ring modulator, Phys. Rev. B 87, 014508 (2013)

[120] J. H. Busnaina , Z. Shi , A. McDonald , D. Dubyna , I. Nsanzineza , J. S. C. Hung , C. W. S. Chang , A. A. Clerk , and C. M. Wilson , Quantum simulation of the bosonic Kitaev chain, Nat. Commun. 15, 3065 (2024)

[121] M. Delbecq , A. Sarlette , M. Mirrahimi , P. Campagne-Ibarcq , T. Kontos , and Z. Leghtas , Dynamically enhancing qubit-photon interactions with antisqueezing, PRX Quantum 5, 020306 (2024)

[122] A. H. Safavi-Naeini , S. Gröblacher , J. T. Hill , J. Chan , M. Aspelmeyer , and O. Painter , Squeezed light from a silicon micromechanical resonator, Nature 500(7461), 185 (2013)

[123] S. Marti , U. von Lüpke , O. Joshi , Y. Yang , M. Bild , A. Omahen , Y. Chu , and M. Fadel , Quantum squeezing in a nonlinear mechanical oscillator, Nat. Phys. 20(9), 1448 (2024)

[124] V. D. Vaidya , B. Morrison , L. G. Helt , R. Shahrokshahi , D. H. Mahler , M. J. Collins , K. Tan , J. Lavoie , A. Repingon , M. Menotti , N. Quesada , R. C. Pooser , A. E. Lita , T. Gerrits , S. W. Nam , and Z. Vernon , Broadband quadrature-squeezed vacuum and nonclassical photon number correlations from a nanophotonic device, Sci. Adv. 6(39), eaba9186 (2020)

[125] R. Nehra , R. Sekine , L. Ledezma , Q. Guo , R. M. Gray , A. Roy , and A. Marandi , Few-cycle vacuum squeezing in nanophotonics, Science 377(6612), 1333 (2022)

[126] Z. Guo , J. Jiang , H. Jiang , J. Ren , and H. Chen , Observation of topological bound states in a double Su−Schrieffer−Heeger chain composed of split ring resonators, Phys. Rev. Res. 3(1), 013122 (2021)

[127] R. Shindou , R. Matsumoto , S. Murakami , and J. Ohe , Topological chiral magnonic edge mode in a magnonic crystal, Phys. Rev. B 87(17), 174427 (2013)

[128] K. Yokomizo and S. Murakami , Non-Bloch bands in two-dimensional non-Hermitian systems, Phys. Rev. B 107(19), 195112 (2023)

[129] N. Macé , F. Alet , and N. Laflorencie , Multifractal scalings across the many-body localization transition, Phys. Rev. Lett. 123(18), 180601 (2019)

Rights & permissions

Higher Education Press