Universal Route for the Emergence of Exceptional Points in PT-Symmetric Metamaterials with Unfolding Spectral Symmetries

We introduce a class of Parity-Time symmetric elastodynamic metamaterials (Ed-MetaMater) whose Hermitian counterpart exhibits a frequency spectrum with unfolding (fractal) symmetries. Our study reveals a scale-free formation of exceptional points (EP) whose density is dictated by the fractal dimension of their Hermitian spectra. Demonstrated in a quasi-periodic Aubry-Harper, a geometric H-tree-fractal, and an aperiodic Fibonacci Ed-MetaMater, the universal route for EP-formation is established via a coupled mode theory model with controllable fractal spectrum. This universality will enable the rational design of novel Ed-MetaMater for hypersensitive sensing and elastic wave control.

A PT-symmetric elastodynamic metamaterial (Ed-MetaMater) has recently been realized by embedding a gain and a lossy mechanical resonators in an elastic medium that facilitates coupling between them [22]. When the intensity of the equal gain/loss and/or the elastic coupling strength between the two resonators of such an Ed-MetaMater are varied, a branch-point singularity forms where the eigenmodes and eigenvalues of the system coalesce. Such degeneracy is known as an exceptional point (EP) and it is the most intriguing feature of PT-symmetric systems: It signifies a transition from a parameter domain where the eigenfrequencies are real and the corresponding eigenmodes of the system respect the PT-symmetry (exact PT-symmetric phase) To further estimate the non-Hermitian that enforces an EP degeneracy for a specific pair of modes, we plot "# vs. the frequency difference between the corresponding mode pairs when = 0 (Δ ' ≡ ∆ "# | ()' ) for each EP found in the spectrum. All EPs (< 25 ) in a Gen-10 PTsymmetric Aubry-Harper Ed-MetaMater are shown in Fig.2(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 "# which coalesce those modes to form an EP. In other words, the non-Hermitian perturbation strength "# 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 Psymmetric Ed-MetaMater. Therefore, a statistical analysis of "# 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 "# .
We evaluate the probability density function (PDF) (Δ ' ) of those EPs. For better statistical processing of these data, we refer to the integrated distribution (Δ ' ) = ∫ ( ) * + % whose derivative (Δ ' ) = − (Δ ' )/ Δ ' determines the PDF of the frequency split of the EPpairs and therefore the PDF ( "# ) of the gain/loss intensity that is necessary for inducing an EP degeneracy. We find that where the best fit parameter is found to be = 0.81 ≈ . Thus the PDF for the gain/loss intensity scales as ( "# ) ∼ "# ,(./0) , which is represented by the flat spread in (Δ ' ) × Δ ' vs. Δ ' (Fig.2(b)).

The PT-symmetric H-tree geometric fractal Ed-MetaMater
We investigate another class of Ed-MetaMater whose fractal spectrum is originating from a geometric fractality in configuration space [42]. This metamaterial is made of two identical planar components in H-motifs ( Fig.1(b)). The first generation of this fractal contains two Hshaped structures (indicated in red and blue)-each made of 3 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). 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 ( Fig.1(b)). The analysis of the correlation-dimension of the frequency spectrum indicates that its converges to 0.80 (Fig.S3 [45]). As previously, the PT-symmetry is created by introducing equal gain/loss ( = 0.001 − 1) to the left and right planar components ( Fig.1(b)). The P-symmetric H-tree-fractal Ed-MetaMater is harmonically excited (0-50 kHz) similar to Aubry-Harper ( Fig.S4 [45]). When is increased, similar to Aubry-Harper, all four generations of the PT-symmetric H-tree-fractal Ed-MetaMater show the emergence of numerous EPs at S "# (%) T which are proportional to the Δ ' associated with those specific EP-pairs ( Fig.2(c)). We evaluated the (Δ ' ) which allows estimating the PDF ( "# ) for the critical gain/loss intensity "# . Fig.2(d) shows the integrated distribution by the variable (Δ ' ) × Δ ' as a function of Δ ' . We find that for = 0.77 ≈ , the data demonstrate a flat spread, leading to the conclusion that ( "# ) ∼ "# ,(./0) . This finding again demonstrates the intimate relation between the emerging EPs and the fractality of the metamaterial's spectrum.

A universal mathematical model for PT-symmetric fractal metamaterials
The intimate relation between ( "# ) 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 Fig.1(d) and [45]). The CMT Fibonacci model is described by the Hamiltonian: (2) where the coupled resonant frequencies % take only two values ± arranged in a Fibonacci sequence and {| ⟩} is the local mode basis. This system is known to have a Cantor-set spectrum with zero Lebesgue measure for all > 0 [47]. Moreover, its spectral fractal dimension can be tuned by varying on-site resonance of the model, . It turns out that the PDF of the nearest level spacing % ≡ %/. − % of such family of systems follows a scale-free distribution whose powerlaw behavior is dictated by the fractality of the spectrum, i.e. ( ) ∼ ,(./0) [48][49][50]. This power law is a signature of level clustering and it is distinct from the PDF ( ) of chaotic or integrable systems [51,52]. We point out that the realization of this class of systems is not confined only to aperiodic systems like the Fibonacci chain model in Eq.(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) [32,53,54], or wave systems with a chaotic classical limit as the kicked Harper model [52]. Therefore, our CMT model represents a typical example of a whole class of systems with fractal spectrum.  2) with its mirror image. The corresponding effective CMT Hamiltonian takes the form: where {`%} is the mirror-symmetric Fibonacci sequence of { % } and describes the coupling between two Fibonacci chains. We found that the P-symmetric variant has a fractal frequency spectrum with the same as the one of the systems of Eq.(2). Finally, a PT-symmetric CMT model #3 is implemented by introducing uniform gain/loss to the left/right portions of the system in Eq.(3), i.e. % → % − and `% → `% + . 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 that connects the two Fibonacci sub-systems increases. For =0, we have two replicas of the same Fibonacci chain in Eq.(2) and, therefore, the spectrum consists of pairs of degenerate modes. As the coupling increases, the degeneracy is lifted % ± = % ± . Simple degenerate perturbation theory with respect to indicates that the new eigenstates are a linear symmetric/antisymmetric combination of the eigenstates of the Fibonacci Hamiltonian in Eq.(2). The above perturbative framework is applicable as long as the is smaller than the distance between nearby frequencies % = %/. − % of the uncoupled Hamiltonian in Eq.(2). The frequency clustering occurring for fractal spectra, however, enforces a rapid breakdown of the perturbation theory, even for infinitesimal . Nevertheless, the eigenstates of the Hamiltonian # ( ) are still eigenfunctions of the P-symmetric operator and therefore are symmetric or antisymmetric with respect to the mirror axis of the total chain. The frequency spacing of nearby levels, however, is not dictated by but the fractal nature of the spectrum. We treat the inclusion of a small non-Hermitian element ± perturbatively. In this case the total Hamiltonian #3 can be written as #3 = # ( ) + Γ where the 2 × 2 perturbation matrix Γ has elements Γ %5 = %5 for ≤ −1 and Γ %5 = − %5 for ≥ 1. Finite leads to level shifts proportional to 2 since the first-order correction vanishes due to the P-symmetry of the corresponding unperturbed eigenmodes of # ( ). For = "# ≃ = Δ ' , the perturbation theory breaks down, signaling level crossing and the appearance of pairs of complex frequencies. We tested the linear relation "# ∼ Δ ' for a variety of -values and find that the linear relation holds with a good approximation in all cases (Fig.3(a)). In case of finite system sizes , some frequency differences Δ ' are still dictated by , though their weight goes to zero at the thermodynamic limit → ∞. 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 ( "# ) ∼ "# ,(./0) . We tested the validity of the above arguments numerically using the Fibonacci CMT model for a variety of potentials and corresponding fractal dimensions ( ) and in all cases we find an excellent agreement with the above theoretical results (Fig.3(b)).
The Fig.4 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 in Fig.S7 and Fig.S8 [45]), 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 Fig.4. 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 ( "# ) ∼ "# ,(./0) in a PTsymmetric Ed-MetaMater with unfolding spectral symmetries can directly be obtained from its spectral fractal dimension . This enables a universal route for effectively predicting the emergence of EPs by the initial spectrum itself.

Conclusions
In summary, we designed three PT-symmetric metamaterials with fractal frequency spectrum-a quasi-periodic Aubry-Harper Ed-MetaMater, an H-tree geometric fractal Ed-MetaMater, and an aperiodic Fibonacci Ed-MetaMater-and investigated them using steady-state dynamic finite element approach. The scale-free emergence of numerous EPs is seen in all metamaterials, showing an intimate relation between the scale-free distribution of critical gain/loss intensities and the spectral fractal dimension of the corresponding Hermitian spectra. Particularly, 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. We further verified those findings from the specific classes of quasi-periodic, fractal, and aperiodic metamaterials and generalized them to a universal law using a CMT-based PT-symmetric fractal mathematical model. The universal relations among the creation of EPs, the scale-free probability distribution of critical gain/loss intensity, and the fractal dimension of the underlying Hermitian spectrum in these PT-symmetric Ed-MetaMater provide a powerful and convenient tool for predicting the emergence of EPs. Our findings are applicable beyond the elastodynamic realm to PT-symmetric metamaterials in acoustic, optical, microwave, and radiofrequency domains as well.

Finite Element Modeling
We used a commercial finite element platform (Abaqus Simulia) to computationally model the steady-state dynamics of the Ed-MetaMater. The material properties of aluminum (Young's modulus: 68.9 GPa, Poisson's ratio: 0.33, density: 2700 kg/m 3 ) is assumed for all components of the Ed-MetaMater which are modeled with cylindrical beam elements (B32). The P-symmetric Ed-MetaMater is harmonically excited using a prescribed axial displacement at the left end of the horizontal coupling rod and the corresponding sinusoidal axial reaction force at its fixed-right-end is measured, simulating a steady-state elastic wave propagation in the metamaterial. The PT-symmetry is created by introducing at each P-symmetric part of the Ed-MetaMater equal amount of energy amplification and attenuation, characterized by an amplification/attenuation rate. These gain/loss mechanisms have been modeled by introducing a structural anti-damping/damping coefficient in Abaqus with its magnitude varying from = 0.001 to 1.

Fibonacci Substitutional Rule and Coupled-mode-theory Modeling for PT-symmetric Fractal Metamaterials
Based on Fibonacci sequence, the number of on-site potentials in the Gen-(n+2) follows: !"# = !"$ + ! with % = 0 and $ = 1. In this mathematical model, we use ' ' and ' ' to represent two different on-site potentials − and + that form the Hamiltonian (correspond to solid and hollow circles in Fig.1(d)).
In order to make the model PT-symmetric, a 'mirrored' chain is coupled with its original counterpart and = ± and = − ± are assumed where ∈ [0, 1] is the on-site potential and the imaginary part is the gain/loss intensity. The sign of is negative for and that are in the upper half the matrix (represents the gain component), and positive for and that are in the lower half the matrix (represents the loss component).
The first five generations of such PT-symmetric chain are shown below:

Chain Sequence
Mirrored Chain Sequence Substitutional Rule → , → → , → = 1 = 2 = 3 = 4 = 5 Fig.1(d) shows an illustration of this coupled model. In matrix notation, these two mirrored chains are orderly located along the diagonal with 1 as the neighboring element on either side and all other elements of the matrix are 0 (here, 1 represents the nearest neighbor coupling). For example, the Hamiltonian of the generations = 3, 4, 5 PT-symmetric fractal models are shown below: