Universal route for the emergence of exceptional points in PT-symmetric metamaterials with unfolding spectral symmetries

The finite element model of the quasi-periodic Aubry-Harper Ed-MetaMater (figure 1(a)) is generated by several vertical beams (length: 6 cm; radius: 0.5 cm) coupled by a horizontal rod with the beam positions given by x_s = sa + R sin(sθ), where $s\in \mathrm{Z}$ is the rod index, a (5 cm) is the distance between the centers of neighboring projection circles, and R (2 cm < a/2) is the radius of the circle. When the projection parameter θ ∈ (0, 1) is an irrational number, the period of the impedance profile of the structure is incommensurate with the lattice period. To this end, we use the ratio of two adjacent numbers in a Fibonacci sequence for θ = p/q, so that the impedance profile of the structure becomes commensurate with the lattice of rods with period q, defining the generation of this Ed-MetaMater (e.g. 7th generation (Gen-7) corresponds to θ = 5/8, which is composed of the 6th and 7th numbers in Fibonacci sequence: 0, 1, 1, 2, 3, 5, 8, 13, ...). The incommensurate limit associated with a truly quasi-periodic structure is investigated via a scaling procedure and it is reached when q → ∞.

The P-symmetric (γ = 0) Aubry-Harper Ed-MetaMater is harmonically excited (0–25 kHz) using a prescribed axial displacement at the left end of the horizontal coupling rod generating a longitudinal wave in the rod, and the corresponding sinusoidal axial reaction force at its fixed-right-end is measured, simulating a steady-state elastic wave dynamics in the metamaterial. The fractal dimension D of its frequency spectra (see figure S1 in supplementary material (https://stacks.iop.org/NJP/23/063079/mmedia)) is calculated using the correlation-dimension method [47] and found to be D = 0.83 ± 0.02 (see figure S2 in supplementary material) (a standard box-counting method results in the same D—albeit the correlation-dimension method converges faster to the value of D). We also find that the D of the P-symmetric Ed-MetaMater is the same as that of the corresponding Ed-MetaMater without being coupled to its mirror image. The D remains the same even in the PT-symmetric case—albeit in this case, it refers only to the real part of the frequencies. This robustness of the D of the real part of frequency spectrum against P or PT-symmetries was checked for all the systems we studied.

The PT-symmetry is created by introducing at each P-symmetric part of the Ed-MetaMater equal amount of energy gain and loss (−/+γ), characterized by an amplification/attenuation rate γ. When γ is increased, several pairs of modes interact and coalesce to form a cascade of EPs at different critical values {${\gamma }_{\mathrm{E}\mathrm{P}}^{(n)}$} and a square-root behavior typical of order-two EPs can be observed near the EP. An example of a typical EP in Gen-10 PT-symmetric Aubry-Harper Ed-MetaMater is shown in figure 2. At a critical gain/loss intensity γ_EP = 0.00251, the eigenvalues (figure 2(a)) and the corresponding eigenmodes (figure 2(c)) coalesce to form an EP. The square-root behavior typical of order-two EPs can be observed in the inset of figure 2(a); Δ_EP is the frequency difference between the corresponding mode pairs. The eigenmode components of the vector u (x_s ) in (figure 2(b)–(d)) describe the displacement field of the transverse tip-deflections of the cross beams at x_s . While in our computation we evaluate the complex displacement fields, the physically relevant information is their amplitude and the phase difference between the components of the mode. We expect that for γ ⩽ γ_EP (exact PT-symmetric phase) these eigenmodes are also eigenmodes of the PT-operator, i.e. the complex field remains invariant under the combined parity (x_s → x_−s ) and time (i → −i) operations (figures 2(b) and (c)). In contrast, for γ > γ_EP (broken PT-symmetric phase) the modes are not any more eigenmodes of the PT-operator. They are rather mapped to one another once the joint PT-operation is applied (figure 2(d)).

To further estimate the non-Hermitian perturbation strength γ that enforces an EP degeneracy for a specific pair of modes, we plot γ_EP vs Δ₀ for each EP found in the spectrum. Here, ${\left.{{\Delta}}_{0}\equiv {{\Delta}}_{\mathrm{E}\mathrm{P}}\right\vert }_{\gamma =0}$ is the frequency difference between the corresponding mode pairs when γ = 0. All EPs (<25 kHz) in a Gen-10 PT-symmetric Aubry-Harper Ed-MetaMater are shown in figure 3(a). Their linear relation demonstrates the intimate relation between the initial (i.e. when γ = 0) frequency split of these two interacting modes Δ₀ and the critical gain/loss intensity γ_EP which coalesce those modes to form an EP. In other words, the non-Hermitian perturbation strength γ_EP that is needed for enforcing a degeneracy between an EP-pair must be of the same order as the frequency split of those modes in the P-symmetric Ed-MetaMater. Therefore, a statistical analysis of γ_EP reduces to the statistical description of these Δ₀, which are associated with the specific mode pairs that eventually form EPs. The latter is easier to evaluate numerically since it does not require a high-resolution parametric evaluation of the modes—as opposed to the precise determination of γ_EP.

Zoom In Zoom Out Reset image size

We evaluate the probability density function (PDF) $\mathcal{P}\left({{\Delta}}_{0}\right)$ of those EPs. For better statistical processing of these data, we refer to the integrated distribution $\mathcal{N}\left({{\Delta}}_{0}\right)={\int }_{{{\Delta}}_{0}}^{\infty }\mathcal{P}(x)\mathrm{d}x$ whose derivative $\mathcal{P}\left({{\Delta}}_{0}\right)=-\mathrm{d}\mathcal{N}\left({{\Delta}}_{0}\right)/\mathrm{d}{{\Delta}}_{0}$ determines the PDF of the frequency split of the EP-pairs and therefore the PDF $\mathcal{P}({\gamma }_{\mathrm{E}\mathrm{P}})$ of the gain/loss intensity that is necessary for inducing an EP degeneracy. We find that

$$\begin{equation}\mathcal{N}\left({{\Delta}}_{0}\right)={\int }_{{{\Delta}}_{0}}^{\infty }\mathcal{P}(x)\mathrm{d}x~{{\Delta}}_{0}^{-k},\end{equation} \tag{ 1 }$$

where the best fit parameter k is found to be k = 0.81 ≈ D. Thus, the PDF for the gain/loss intensity scales as $\mathcal{P}\left({\gamma }_{\mathrm{E}\mathrm{P}}\right)\sim {\gamma }_{\mathrm{E}\mathrm{P}}^{-\left(1+D\right)}$, which is represented by the flat spread in $\mathcal{N}\left({{\Delta}}_{0}\right){\times}{{\Delta}}_{0}^{k}$ vs Δ₀ (figure 3(b)).

We investigate another class of Ed-MetaMater whose fractal spectrum is originating from a geometric fractality in configuration space [44]. The H-tree geometric fractal Ed-MetaMater is made of two identical planar components made of H-motifs (figure 4(a)). The first generation of this fractal contains two H-shaped structures—each made of three identical cylindrical beams (length: 11.6 cm, radius: 2.38 mm)—coupled by a passive (zero gain/loss) horizontal elastic rod (coupling length: 11.6 cm; length of exterior side ledges: 5.8 cm) such that they form a P-symmetric system. Each subsequent generation adds H-motifs scaled down in length by a factor of 2 (constant diameter) to each tip of the prior H-structure. The analysis of the correlation-dimension of the frequency spectrum indicates that its D converges to 0.80 (see figure S3 in supplementary material).

Zoom In Zoom Out Reset image size

As previously, the PT-symmetry is created by introducing equal gain/loss (γ = 0.001–1) to the left and right planar components (figure 4(a)). The P-symmetric H-tree-fractal Ed-MetaMater is harmonically excited (0–50 kHz) using a prescribed axial displacement at the left end of the horizontal coupling rod generating a longitudinal wave in the rod, and the corresponding sinusoidal axial reaction force at its fixed-right-end is measured, simulating a steady-state elastic wave dynamics in the metamaterial (see figure S4 in supplementary material). When γ is increased, similar to Aubry-Harper Ed-MetaMater, all four generations of the PT-symmetric H-tree-fractal Ed-MetaMater show the emergence of numerous EPs at $\left\{{\gamma }_{\mathrm{E}\mathrm{P}}^{(n)}\right\}$ which are proportional to the Δ₀ associated with those specific EP-pairs (figure 4(b)). We evaluated the $\mathcal{N}\left({{\Delta}}_{0}\right)$ which allows estimating the PDF $\mathcal{P}\left({\gamma }_{\mathrm{E}\mathrm{P}}\right)$ for the critical gain/loss intensity γ_EP. Figure 4(c) shows the integrated distribution by the variable $\mathcal{N}\left({{\Delta}}_{0}\right){\times}{{\Delta}}_{0}^{k}$ as a function of Δ₀. We find that for k = 0.77 ≈ D, the data demonstrate a flat spread, leading to the conclusion that $\mathcal{P}\left({\gamma }_{\mathrm{E}\mathrm{P}}\right)\sim {\gamma }_{\mathrm{E}\mathrm{P}}^{-(1+D)}$. This finding again demonstrates the intimate relation between the emerging EPs and the fractality of the metamaterial's spectrum.

To further verify the universal nature of equation (1) we studied another class of Ed-MetaMater with unfolding spectral symmetries—an aperiodic system based on Fibonacci substitutional rule (figure 5(a), details in supplementary material). It is created by equally spaced (3 cm) vertical beams (lengths of A: 8 cm, B: 6 cm; radius: 0.5 cm) organized based on Fibonacci substitutional rule (further details in supplementary material) that are coupled by a horizontal rod (radius: 0.5 cm; exterior ledges: 3 cm) and then mirrored to form a P-symmetric system. The system is harmonically excited (0–50 kHz) using a prescribed axial displacement at the left end of the horizontal coupling rod generating a longitudinal wave in the rod, and the corresponding sinusoidal axial reaction force at its fixed-right-end is measured, simulating a steady-state elastic wave dynamics in the metamaterial (see figure S5 in supplementary material). The correlation-dimension analysis indicates that the frequency spectrum of this system is characterized by D = 0.80 (see figure S6 in supplementary material). The relation between Δ₀ and γ_EP is found to be linear again (figure 5(b)). The $\mathcal{N}\left({{\Delta}}_{0}\right){\times}{{\Delta}}_{0}^{k}$ vs Δ₀ demonstrates a flat spread with k = 0.80 ≈ D (figure 5(c)), concluding that $\mathcal{P}\left({\gamma }_{\mathrm{E}\mathrm{P}}\right)\sim {\gamma }_{\mathrm{E}\mathrm{P}}^{-(1+D)}$.

Zoom In Zoom Out Reset image size

The intimate relation between $\mathcal{P}\left({\gamma }_{\mathrm{E}\mathrm{P}}\right)$ and the spectral fractal dimension of an Ed-MetaMater at γ = 0 implies the existence of an underlying universal route for the creation of EPs in systems with fractal spectrum. To this end, we develop a CMT-based model that utilizes on-site resonant modes that follow an aperiodic Fibonacci substitutional rule (details in figure 6(a) and supplementary material). The CMT Fibonacci model is described by the Hamiltonian:

$$\begin{equation}H=\sum\limits _{n}\vert n\rangle {V}_{n}\langle n\vert +\sum\limits _{n}\vert n+1\rangle \langle n\vert +\mathrm{c}.\mathrm{c}.,\end{equation} \tag{ 2 }$$

where the coupled resonant frequencies V_n take only two values ±V arranged in a Fibonacci sequence and {|n⟩} is the local mode basis. This system is known to have a Cantor-set spectrum with zero Lebesgue measure for all V > 0 [48]. A benefit of the Fibonacci CMT modeling is that its spectral fractal dimension D can be tuned by varying on-site resonance of the model, V, thus giving us the possibility to scrutinize the relation $\mathcal{P}\left({\gamma }_{\mathrm{E}\mathrm{P}}\right)\sim {\gamma }_{\mathrm{E}\mathrm{P}}^{-\left(1+D\right)}$ for a variety of spectral fractal dimensions. It turns out that the PDF of the nearest level spacing s_n ≡ ω_n+1 − ω_n of such family of systems follows a scale-free distribution whose power-law behavior is dictated by the fractality of the spectrum, i.e. $\mathcal{P}\left(s\right)\sim {s}^{-\left(1+D\right)}$ [49–51]. This power law is a signature of level clustering and it is distinct from the PDF $\mathcal{P}\left(s\right)$ of chaotic or integrable systems [52, 53]. We point out that the realization of this class of systems is not confined only to aperiodic systems like the Fibonacci chain model in equation (2), but also applicable to quasi-periodic systems with metal-insulator transition at some critical value of the on-site resonance (e.g. the Aubry-Harper model) [34, 54, 55], or wave systems with a chaotic classical limit as the kicked Harper model [53]. Therefore, our CMT model represents a typical example of a whole class of systems with fractal spectrum.

Zoom In Zoom Out Reset image size

The P-symmetric Fibonacci model is implemented by coupling the Hamiltonian of equation (2) with its mirror image. The corresponding effective CMT Hamiltonian takes the form:

$$\begin{align}{H}_{\mathrm{P}}& =\left(\sum\limits _{n=-N}^{-1}\left\vert n,\rangle ,{V}_{n},\langle ,n\right\vert +\sum\limits _{n=-N}^{-2}\left\vert n,+,1,\rangle ,\langle ,n\right\vert +\mathrm{c}.\mathrm{c}.\right)+\left({\sum }_{n=1}^{N}\left\vert n,\rangle ,{\overline{V}}_{n},\langle ,n\right\vert +{\sum }_{n=1}^{N-1}\left\vert n,+,1,\rangle ,\langle ,n\right\vert +\mathrm{c}.\mathrm{c}.\right)\\ & \quad +\left(\vert 1\rangle \langle -1\vert +\mathrm{c}.\mathrm{c}.\right),\end{align} \tag{ 3 }$$

where {$\left.{\bar{V}}_{n}\right\}$ is the mirror-symmetric Fibonacci sequence of {V_n } and t describes the coupling between two Fibonacci chains. We found that the P-symmetric variant has a fractal frequency spectrum with the same D as the one of the systems of equation (2). Finally, a PT-symmetric CMT model H_PT is implemented by introducing uniform gain/loss to the left/right portions of the system in equation (3), i.e. V_n → V_n − iγ and ${\bar{V}}_{n}\to {\bar{V}}_{n}+i\gamma$. Because of the simplicity in its structure, this model allows reaching higher generations for more accurate numerical analyses compared to the computationally costly finite element models in previous three examples.

Consider the parametric evolution of frequencies of the P-symmetric model as the coupling constant t that connects the two Fibonacci sub-systems increases. For t = 0, we have two replicas of the same Fibonacci chain in equation (2) and, therefore, the spectrum consists of pairs of degenerate modes. As the coupling t increases, the degeneracy is lifted ${\omega }_{n}^{{\pm}}={\omega }_{n}{\pm}t$. Simple degenerate perturbation theory with respect to t indicates that the new eigenstates are a linear symmetric/antisymmetric combination of the eigenstates of the Fibonacci Hamiltonian in equation (2). The above perturbative framework is applicable as long as the t is smaller than the distance between nearby frequencies s_n = ω_n+1 − ω_n of the uncoupled Hamiltonian H in equation (2). The frequency clustering occurring for fractal spectra, however, enforces a rapid breakdown of the perturbation theory, even for infinitesimal t. Nevertheless, the eigenstates of the Hamiltonian H_P(t) are still eigenfunctions of the P-symmetric operator and therefore are symmetric or anti-symmetric with respect to the mirror axis of the total chain. The frequency spacing of nearby levels, however, is not dictated by t but the fractal nature of the spectrum.

We treat the inclusion of a small non-Hermitian element ±γ perturbatively. In this case the total Hamiltonian H_PT can be written as H_PT = H_P(t) + iγΓ where the 2N × 2N perturbation matrix Γ has elements Γ_nm = δ_nm for n ⩽ −1 and Γ_nm = −δ_nm for n ⩾ 1. Finite γ leads to level shifts proportional to γ² since the first-order correction vanishes due to the P-symmetry of the corresponding unperturbed eigenmodes of H_P(t). For γ = γ_EP ≃ s = Δ₀, the perturbation theory breaks down, signaling level crossing and the appearance of pairs of complex frequencies. It is still intriguing the fact that the non-Hermitian perturbation operator iγΓ couples the nearby levels of H_p(t) in the case of fractal spectrum where the validity of level spacing, and therefore of perturbation theory, is 'blurred'—specifically in the thermodynamic N → ∞ limit. Nevertheless, our detailed numerical investigations confirmed the linear relation γ_EP ∼ Δ₀ for a variety of V-values and find that the linear relation holds with a good approximation in all cases (figure 6(b)). In case of finite system sizes N, some frequency differences Δ₀ are still dictated by t, though their weight goes to zero at the thermodynamic limit N → ∞. The above analysis allows us to associate the PDF of the gain/loss intensity that results in EP degeneracy with the distribution of level spacings, leading to the conclusion that $\mathcal{P}\left({\gamma }_{\mathrm{E}\mathrm{P}}\right)\sim {\gamma }_{\mathrm{E}\mathrm{P}}^{-\left(1+D\right)}$. We tested the validity of the above arguments numerically using the Fibonacci CMT model for a variety of potentials V and corresponding fractal dimensions D(V) and in all cases we find an excellent agreement with the above theoretical results (figure 6(c)).

The figure 7 comprehensively presents the relationship between the spectral fractal dimensions of all aforementioned P-symmetric systems and the power exponents corresponding to the EPs in the PT-symmetric Fibonacci CMT model with different on-site resonances (indicated by circles; further details see figures S7 and S8 in supplementary material), the PT-symmetric Aubry-Harper Ed-MetaMater (blue star), the PT-symmetric H-tree-fractal Ed-MetaMater (yellow triangle), and the PT-symmetric Fibonacci Ed-MetaMater (purple square). The universality in the relations between the emergence of EPs in these metamaterials and the fractality of their initial spectra is evident in figure 7. The linear fit (black line) with a slope ∼1 signifies the universality of this equality relationship, i.e., the power-law exponent describing a scale-free PDF $\mathcal{P}\left({\gamma }_{\mathrm{E}\mathrm{P}}\right)\sim {\gamma }_{\mathrm{E}\mathrm{P}}^{-(1+D)}$ in a PT-symmetric Ed-MetaMater with unfolding spectral symmetries can directly be obtained from its spectral fractal dimension D. This enables a universal route for effectively predicting the emergence of EPs by the initial spectrum itself.

Zoom In Zoom Out Reset image size

The linear relation found between the critical gain/loss required for creating EPs and the initial split between the mode pairs that coalesce, shows that the high-signal-quality hypersensitive sensors that exploit EPs in PT-symmetric metamaterials can be engineered by appropriate interacting mode pairs that facilitate experimentally realizable low gain/loss. Such systems can be realized via active materials, for example by using piezoelectric elements embedded in the crossbeams that are controlled by non-Foster circuits to provide balanced gain/loss [56] or in combination with passive materials with highly tunable loss [57]. Gain in ultrasonic frequency regime can also be achieved through electroacoustic amplification via phonon-electron interaction in piezoelectric semiconductor [58–60]. The classical noise in Ed-MetaMater (in contrast to quantum noise) may not be a significant concern towards achieving high signal-to-noise ratio as it has been indicated by recent studies on non-Foster circuits [61] that the noise is in same level as in typical diabolic degeneracies.