Abnormal Frequency Response Determined by Saddle Points in Non-Hermitian Crystal Systems

Abstract

In non-Hermitian crystal systems under open boundary condition (OBC), it is generally believed that the OBC modes with frequencies containing positive imaginary parts, when excited by external driving, will experience exponential growth in population, thereby leading to instability. However, our work challenges this conventional understanding. In such a system, we find an anomalous response that grows exponentially with the frequency aligned with those of saddle points. The frequencies of these saddle points on the complex plane are below the maximum imaginary part of OBC spectrum, but they can lie within or beyond the OBC spectrum. We derive general formulas of excitation–response relationships and find that this anomalous response can occur because the excitation of OBC modes eventually evolve toward these saddle points at long times. Only when the frequencies of all these saddle points are below the real axis do the non-Hermitian crystal systems remain stable under periodic excitation. Thus our results also provide new insights on the stability criterion of non-Hermitian crystal systems.

Introduction.—The study of response to excitation is a central issue in many physical models, ranging from forced vibration in classic mechanics [1] and transient response in circuits [2] to the dynamic evolution of wavefunction in quantum physics [3]. For a linear control system subjected to periodic excitation, it is generally believed that the response frequency is intimately connected to the excitation frequency, the system’s natural frequencies or a combination of both[1, 2]. This understanding also uncovers the stability criteria [4, 5]. For a linear control system to be stable, the natural frequencies may have imaginary parts but these imaginary parts must not be greater than zero; otherwise, the natural response would grow exponentially over time, thereby undermining the system’s stability.

In recent years, with the rapid development of non-Hermitian physics [6, 7, 8, 9, 10, 10, 11, 12, 13, 14, 15, 16], the issue of excitation–response relationships in non-Hermitian crystal systems, as shown in Fig. 1, has attracted extensive attentions [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 20, 23, 27]. The Green’s function is a crucial tool for studying this issue, and the explicit formulas of its non-Hermitian version [28, 29, 20, 30, 31, 32, 33] recently proposed by Xue et. al. [20] revealed a potential application of non-Hermitian systems as directional amplifiers. This non-Hermitian version primarily focused on the response from the external driving, also known as forced response, which is appropriate when large on-site dissipation is applied to suppress the contribution from the natural response [20]. However, without the constraint of Hermiticity, the natural frequencies of non-Hermitian crystal systems under OBC, i.e., OBC spectrum, can possess positive imaginary parts even in lossy lattices [27]. In such a case, a recent quantum walk experiment [34] indicated that the dynamic evolution of probabilities of the photon is dominated by the frequency with the largest imaginary part in the OBC spectrum. However, another study [35] by Stefano Longhi pointed out that the bulk dynamic evolution of initial states is limited by an upper bound determined by spectrum under periodic boundary condition. There is still some controversy regarding the excitation–response relationships of non-Hermitian systems, and to the best of our knowledge, a unified description has yet to be established.

Figure 1: A schematic diagram of a non-Hermitian crystal system undergoing periodic excitation $e^{-i\omega t}$. Here, the non-Hermiticity is characterized by on-site dissipation $i\kappa$ and non-conjugated hopping $t_{n}\neq t_{-n}^{*}$, where $t_{n}$ represents the hopping strength from site $i+n$ to site $i$ (for any $i$).

In this letter, for a non-Hermitian crystal system subjected to periodic excitation, we observe an anomalous exponentially growing response, with its frequency aligned with frequency of a saddle point beyond the range of OBC spectrum. By utilizing the Laplace transform in the dynamic equation, we derive general formulas of excitation–response relationships, successfully explaining this phenomenon and clarifying the aforementioned controversy. Our formulas demonstrate that the complete response can be divided into forced response by external excitation and natural response by natural frequencies, but the long-time behavior of complete response is determined by their competition. Intriguingly, the frequency components of this natural response are initially composed of the OBC spectrum but eventually evolve to the frequencies of saddle points. These saddle points’ frequencies can lie within or beyond the OBC spectrum. From this, we infer that for a non-Hermitian crystal system to be stable undergoing excitation, all of these saddle points must be located below the real axis.

Figure 2: A concrete model with parameters $\kappa=0$, $t_{3}=3,t_{2}=t_{1}=t_{-1}=2,t_{-2}=-2$, and $t_{-3}=1$ is illustrated, with the numerical calculation of OBC spectrum for $N=1600$ sites shown in (a), and GBZ and aGBZ presented in (b). In these diagrams, the saddle points are represented by pentagram, while the uppermost saddle point (USP) with the largest imaginary part is highlighted by purple pentagram in (a) and the corresponding wavevector is shown in (b). Under excitation at site 800, the evolution of amplitude $\left|a(t)\right|$ at site 750 is depicted with different on-site dissipation: (c) $\kappa=-4.6$, $\kappa=-3.9$ and (d) $\kappa=-3.2$. The amplitude in (d) exhibits exponential growth with frequency aligned with the frequency of the USP in (a).

Abnormal frequency response.—The schematic diagram of a non-Hermitian crystal system subjected to periodic excitation $e^{-i\omega t}$ with $\omega$ being the excitation frequency is shown in Fig. 1. We first attempt to exhibit counterintuitive results in Fig. 2 using a concrete non-Hermitian model, which takes a non-Bloch Hamiltonian

$$h(\beta)=3\beta^{3}+2\beta^{2}+2\beta+i\kappa+2\beta^{-1}-2\beta^{-2}+\beta^{-3},$$ (1)

where $\kappa$ is the on-site uniform dissipation. Here, $\beta$ is the allowed wavevector under OBC and the corresponding frequencies are the OBC spectrum $E_{OBC}$. The allowed wavevector and OBC spectrum of non-Hermitian system are dramatically different from those under periodic boundary condition due to their striking sensitivity to the boundary [36, 37, 38, 39]. When $\kappa=0$, Fig. 2(a) shows the complex OBC spectrum, which is symmetric about the real axis. The modes below the real axis have negative imaginary parts that denote a finite lifetime [40, 41]. In contrast, the modes with positive imaginary parts typically imply their population grows exponentially over time [42].

The OBC spectrum $E_{OBC}$ becomes continuous for infinite-size non-Hermitian crystal systems. While to form continuum spectrum, the allowed wavevector i.e., the roots of characteristic function $E_{OBC}-h(\beta)=-3\prod_{i}\left(\beta-\beta_{i}\right)/\beta^{3}=0$ must fulfill the conditions $\left|\beta_{3}\right|=\left|\beta_{4}\right|$, where $\beta_{i}$ and $\beta_{3,4}$ label the $i$-th root and $3,4$-th root sorted by their moduli, respectively [9, 8]. The collection of these allowed wavevector on the complex plane is the so-called generalized Brillouin zone (GBZ), as shown in Fig. 2(b). A standard numerical method to calculate two roots with equal modulus typically yields the auxiliary GBZ (aGBZ) $\left|\beta_{i}\right|=\left|\beta_{j}\right|$ (for any $i,j$) that comprises the GBZ [10, 11], as shown in Fig. 2(b).

Now we turn to the response as the periodic excitation is applied. The dynamic evolution for model described by Hamiltonian Eq. (1) is governed by [43]

$$i\frac{d\mathbf{a}}{dt}=H\mathbf{a}+\mathbf{b}_{in}e^{-i\omega t}$$ (2)

where $H$ is the matrix form of non-Hermitian Hamiltonian. Here, $\mathbf{a}=(a_{1}(t),a_{2}(t),\cdots,a_{N}(t))^{T}$ is the amplitude vector of the lattice and $\mathbf{b}_{in}$ is the input amplitude. Under excitation at site 800, the response amplitude $\left|a(t)\right|$ at site 750 with different dissipation are shown in Fig. 2(c-d). Only when $\kappa=-3.2$, the amplitude grows exponentially over time. A similar trend in the dynamic evolution of amplitude can be found at other sites as shown in Fig. S1 in Supplemental Materials. With large on-site uniform dissipation $\kappa=-4.6$, which ensures all the OBC spectrum lies below the real axis (Fig. 2(a)), Fig. 2(c) shows that the evolution of $\left|a(t)\right|$ remains unchanged after a small initial increase. From Eq. 2, the excitation $\mathbf{b}_{in}e^{-i\omega t}$ can be decomposed by a superposition of eigenvector, thus leading to the activation of a series of the OBC spectrum. However, this natural response decay with time due to their negative imaginary energy while the forced response survives under long-time evolution. Altering the dissipation to $\kappa=-3.9$ renders the maximum imaginary value of OBC spectrum (denoted by $\operatorname{Im}\left(E_{m}\right)$) becoming greater than zero. This change leads to a rapid amplification of response amplitude, but which eventually stabilize. These results suggest that the OBC modes with positive imaginary parts are indeed excited for a short period but then disappear.

However, when $\kappa=-3.2$, Fig. 2(d) shows an exponential growth of amplitude with specified frequency after a short-term growth and oscillation. This growth rate is precisely equal to the imaginary part of the uppermost saddle point (USP, represented by purple pentagram in Fig. 2(a)), which has maximum imaginary part among saddle points below the $\operatorname{Im}\left(E_{m}\right)$. Mathematically, the saddle points satisfy $\partial h(\beta)/\partial\beta|_{\beta_{s}}=0$, which also corresponds to the occurrence of multiple roots [44]. Therefore, they are always located on the aGBZ, but they may not be found on the GBZ [44]. The frequency of USP neither lies within the OBC spectrum nor alignes with the external driving frequency. Therefore, this exponential growth amplitude is referred to as abnormal frequency response.

General formulas of excitation–response relationships.—To gain deeper insight into the abnormal frequency response mentioned above, we turn to derive the excitation–response relationships for a general non-Hermitian Hamiltonian $h(\beta)=\sum_{n=-l}^{r}t_{n}\beta^{n}$. Considering the possibility of exponential growth in amplitude, we first introduce the Laplace transform in Eq. 2

$$\mathbf{a}(p)=\frac{1}{ip-H}\frac{1}{p+i\omega}\mathbf{b}_{in},$$ (3)

where $p$ is a complex variable that contains real $s>s_{0}$ to ensure convergence with $s_{0}$ being the abscissa of absolute convergence. The relationship between response and excitation can be defined as the Laplace transform of Green function

$$G(p)=\frac{1}{ip-H}\frac{1}{p+i\omega},$$ (4)

with the matrix element $[G(p)]_{jk}$ denoting the response at site $j$ to the excitation at site $k$. For an infinite-size system, such matrix can be given by inverting Toeplitz matrices as proposed by Xue et. al. [45, 20, 46]

$$[G(p)]_{jk}=\frac{1}{2\pi i}\oint_{\beta=GBZ}\frac{\beta^{j-k-1}}{ip-h(\beta)}\frac{1}{p+i\omega}\mathrm{d}\beta,$$ (5)

where the integral loop is proven to the GBZ [20]. Then the time evolution of system can be obtained through inverse Laplace transform

$$[G(t)]_{jk}=$$

$$\left(\frac{1}{2\pi i}\right)^{2}\int_{s-i\infty}^{s+i\infty}\oint_{\beta=GBZ}\frac{\beta^{j-k-1}}{ip-h(\beta)}\frac{1}{p+i\omega}e^{pt}\mathrm{d}\beta\mathrm{d}p.$$ (6)

This integral can be evaluated by calculating the residue at pole $ip-h(\beta)=0$, which requires expressing the root as a function of $p$. For a simple Hatano-Nelson model [47], the explicit analytical solution of dynamic evolution is derived, which is consistent with the numerical calculations as detailed in Supplemental Materials. However, when $h(\beta)$ involves long-range couplings beyond nearest neighbors, expressing the root as a function of $p$ can be challenging or impossible. To address this challenge, we change the order of integration, which simplifies the evaluation of this integral and makes the physical interpretation explicit, as discussed next.

Clearly, there are two kinds of poles for variable $p$ that contribute to the response. The contribution from pole $p+i\omega=0$ arises from the external driving

$$[G(t)]_{jk}(E)=\frac{1}{2\pi i}\oint_{\beta=GBZ}\frac{\beta^{j-k-1}}{\omega-h(\beta)}e^{-i\omega t}\mathrm{d}\beta,$$ (7)

which corresponds to the forced response. By further investigating the root $w-h(\beta)=\prod_{i}-t_{r}\left(\beta-\beta_{i}(\omega)\right)/\beta^{l}=0$, we have

$$[G(t)]_{jk}(E)$$

$$= \left\{\begin{aligned} &\sum_{n=1}^{l}\frac{-\beta_{n}(\omega)^{j-k+l-1}}{t_{r}\prod_{k\neq n}[\beta_{n}(\omega)-\beta_{k}(\omega)]}e^{-i\omega t},j-k+l-1\geqslant 0,\\ &\sum_{n=l+1}^{l+r}\frac{\beta_{n}(\omega)^{j-k+l-1}}{t_{r}\prod_{k\neq n}[\beta_{n}(\omega)-\beta_{k}(\omega)]}e^{-i\omega t},j-k+l-1<0.\end{aligned}\right.$$ (8)

For the response sites far from the excitation site, the above equation can be simplified as $[G(t)]_{jk}(E)\sim\beta_{l}(\omega)^{j-k}$ when $j>>k$ and $[G(t)]_{jk}(E)\sim\beta_{l+1}(\omega)^{j-k}$ when $j<<k$. These results return to the explicit Green’s functions, as detailed in previous studies [20, 30]. Now we turn to investigate the contribution from pole $ip-h(\beta)=0$

$$[G(t)]_{jk}(N)=\frac{1}{2\pi i}\oint_{\beta=GBZ}-\frac{\beta^{j-k-1}}{\omega-h(\beta)}e^{-ih(\beta)t}\mathrm{d}\beta,$$ (9)

which represents the response arising from the natural frequencies, dubbed natural response. The natural response includes the contribution from all the frequencies of system. This term may result in the exponential growth of amplitude since $h(\beta)$ can have positive imaginary parts. It should be noted that all the frequencies of system are not equivalent to the OBC spectrum, though the integral zone is GBZ in Eq. (9), since the frequencies may alter during time evolution. From the OBC spectrum, the upper limit of this integral can be determined

$$\absolutevalue{[G(t)]_{jk}(N)}\leqslant e^{\operatorname{Im}\left(E_{m}\right)t}\frac{1}{2\pi}\oint_{\beta=\mathrm{GBZ}}\absolutevalue{\frac{\beta^{j-k-1}}{\omega-h(\beta)}}\mathrm{d}\beta,$$ (10)

since $\operatorname{Im}\left(h(\beta)\right)\leq\operatorname{Im}\left(E_{m}\right)$.

It should be noted the Taylor expansion of the exponential term in Eq. (9) gives rise to infinite essential singularities, posing a challenge to residue theorem. At long times, the asymptotic expansion of Eq. (9) can be obtained using the method of steepest descent from textbook [48, 49]. According to this method, the integral path can be designed to pass through all the saddle points, and the asymptotic expansion is mainly contributed by the surrounding regions of these saddle points. Although deforming the integration path may across some pole $h(\beta)=\omega$, this term introduces an oscillatory error with constant magnitude and can therefore be safely neglected. Consequently, the response by the natural frequencies becomes the sum of contributions from all these saddle points, yielding the asymptotic form of the time-domain Green’s function as follows:

$$[G(t)]_{jk}(N)\sim\sum_{s=1}^{S}\frac{({\beta_{s}})^{j-k-1}}{h(\beta_{s})-\omega}e^{-ih(\beta_{s})t}\sqrt{\frac{i}{2\pi h^{\prime\prime}(\beta_{s})t}}.$$ (11)

Here, $h^{\prime\prime}(\beta_{s})$ is the second derivative [48] at saddle points and $S$ represents the number of saddle points below $\operatorname{Im}\left(E_{m}\right)$. For simplicity, we assume that all the involved saddle points are second order and $h(\beta_{s})-\omega\neq 0$ to avoid resonance. Eq. (Abnormal Frequency Response Determined by Saddle Points in Non-Hermitian Crystal Systems) and Eq. (11) indicate the frequency components of natural response initially consist of the OBC spectrum and eventually evolve towards saddle points’ frequencies at long times. It is expected that the USP will dominate the long-time dynamic behavior because it has maximum positive imaginary part among all the saddle points. As pointed out earlier, such saddle points must lie on the aGBZ, but they may not necessarily lie on the GBZ.

Figure 3: Spatial distribution of response near the excitation site 800 at $t=100$ for model Eq. 1 with different on-site dissipation (a) $\kappa=-3.9$ and (b) $\kappa=-3.2$. Inset in (a) shows the two roots $\beta_{3}(\omega),\beta_{4}(\omega)$ of forced frequency $\omega$ close to GBZ. Besides these two roots, inset in (b) also shows $\beta_{s}$ at USP.

After deriving the response-excitation relationships, we can now explain the result in Fig. 2 in detail. When $\kappa=-4.9$ and $\kappa=-3.9$, the saddle points’ frequencies below $\operatorname{Im}\left(E_{m}\right)$ have only negative imaginary parts, regardless of whether $\operatorname{Im}\left(E_{m}\right)$ is above the real axis. Consequently, the natural response decays and the forced response survives. However, when $\kappa=-3.2$, the frequency of USP have a positive imaginary part, then the amplitude is exponentially growing characterized by this frequency.

The forced response and the natural response also exhibit significant difference in spatial distribution as shown in Fig. 3. When $\kappa=-4.9$, the forced response dominates. The response at right side of excitation site $[G(t)]_{ij}(E)\sim\beta_{4}^{i-j}$, while at the left side $[G(t)]_{ij}(E)\sim\beta_{3}^{i-j}$, in agreement with the theoretical prediction by Eq. (Abnormal Frequency Response Determined by Saddle Points in Non-Hermitian Crystal Systems), demonstrating a potential application as directional amplifiers. However, when $\kappa=-3.2$, the natural response governs the dynamic evolution. The numerical calculation of spatial distribution of response on both sides of the excitation point follows same scaling relation $[G(t)]_{ij}(I)\sim\beta_{s}^{i-j}$, consistent with Eq. (11). These results also confirm the excitation–response relationships from another perspective.

Figure 4: (a) The OBC spectrum of corrected model Eq. 1 with each parameter multiplied by $e^{-i\frac{\pi}{4}}$. (b) The dynamic evolution of amplitude aligns with the frequency of new UPS in (a). (c) The OBC spectrum and (d) the variation of augment of amplitude with time for corrected Hatano-Nelson model $h(\beta)=1.8\beta+\beta^{-1}+0.4i$. The numerical result is denoted by red line and the theoretical result is denoted by blue line.

The competition among saddle points.— The general formulas of excitation–response relationships also explain the observation of $\operatorname{Im}\left(E_{m}\right)$-dominated dynamic evolution of probabilities of the photon experimentally since where $\operatorname{Im}\left(E_{m}\right)$ corresponds to the frequency of USP [34, 44]. Such a scenario has its counterpart in corrected model Eq. (1) by multiplying same factor $e^{-i\frac{\pi}{4}}$ on each parameter, as shown in Fig. 4(a). Here, the OBC spectrum is rotated by 45 degree and the imaginary part of USP is $\operatorname{Im}\left(E_{m}\right)$. Considering the same excitation and response site as in Fig. 2, the response amplitude is exponential growth with frequency matching that of this new USP as shown in Fig. 4(b). This result reflects the competition among saddle points as source of response frequency.

We further investigate the scenario where the OBC spectrum possesses the same positive imaginary parts as exemplified by corrected Hatano-Nelson model in Fig. 4(c). Here, the two ends of OBC spectrum are the frequencies of saddle points and both saddle points are USP. Unlike the previous focus on the temporal growth of the amplitude, here we investigate the argument of response amplitude, as shown in Fig. 4(d). The theoretical results involving the superposition of contribution from these two saddle points given by Eq. (11) are essentially consistent with the numerical results in oscillation period. Fig. 4(d) also shows phase difference between theoretical and numerical results. This difference do not affect our conclusions since the asymptotic form in Eq. (11) only provides the long-time behavior without involving any initial phase.

Conclusions.—In conclusion, we theoretically derive general formulas of excitation–response relationships for non-Hermitian crystal systems under periodic excitation. Our results indicate the response contains the forced response from external driving and natural response from frequencies of system. Intriguingly, the latter initially consists of OBC spectrum but evolves towards the saddle point below the $\operatorname{Im}\left(E_{m}\right)$ at long times. These saddle points can lie within or beyond the OBC spectrum. Consequently, the long-time behavior of response is determined by the competition between forced response and natural response from these saddle points. When all these saddle point lie below the real axis, the forced response dominates; otherwise, the system’s amplitude exhibits exponential growth with the frequency of USP, leading to instability. Thus, our theory provides a comprehensive picture of the dynamic evolution and offers new insights into stability criterion of non-Hermitian crystal systems.

This work is supported by the Natural Science Foundation of Hunan Province (2024JJ6011).