Surface-Wave Propagation on Non-Hermitian Metasurfaces with Extreme Anisotropy

The problem geometry is schematically illustrated in Fig. 1a. We consider, in free space, a metasurface of infinite extent in the $x-y$ plane, featuring a 1-D periodic modulation of the surface conductivity along the $y$-direction, with alternating values $\sigma_{a}$ and $\sigma_{b}$ (and sub-periods $d_{a}$ and $d_{b}$, respectively). Assuming an implicit $\exp(-i\omega t)$ time-harmonic dependence, the generally complex-valued conductivities are written as

with the prime and double-prime symbols tagging the real and imaginary parts, respectively. Throughout the study, we focus on the parameter regime

In view of the assumed time-dependence, the first condition implies that the “$a$”- and “$b$”-type constituents exhibit loss and gain, respectively. The second condition instead guarantees that both constituents are either of inductive ($\sigma^{\prime\prime}_{a,b}>0$) or capacitive ($\sigma^{\prime\prime}_{a,b}<0$) nature; this rules out the possibility of a hyperbolic response, which is not of interest here since it has already been studied [19, 20].

We highlight that the surface-conductivity parameterization in (1) and (2) is especially suited for 2-D materials (e.g., graphene), and can be readily related to more conventional constitutive parameters. For instance, it can be obtained from the surface impedance as [42]

with the factor 2 accounting for the two-faced character of the sheet. Moreover, for a very thin dielectric layer of thickness $\Delta\ll\lambda$ (with $\lambda=2\pi c/\omega$ denoting the free-space wavelength and $c$ the corresponding wavespeed), it can be approximately related to the relative permittivity via [9, 43]

where $k=\omega/c=2\pi/\lambda$ and $\eta=\sqrt{\mu_{0}/\varepsilon_{0}}$ denote the free-space wavenumber and characteristic impedance, respectively.

In [44], it was shown that, in the limit of deeply subwavelength modulation periods $d=d_{a}+d_{b}\ll\lambda$, a structure like the one in Fig. 1a could be effectively modeled by a homogeneous, uniaxially anisotropic effective surface conductivity with relevant components

with $f_{a}=d_{a}/d$ and $f_{b}=1-f_{a}$ denoting the filling fractions. These mixing formulae closely resemble those occurring in the EMT modeling of multilayered volumetric metamaterials [45].

We are interested in studying the propagation of surface waves along this metasurface. In mathematical terms, this entails finding nontrivial source-free solutions $\propto\exp\left[i\left(k_{x}x+k_{y}y+k_{z}z\right)\right]$ which are evanescent along the $z$-direction and propagating in the $x-y$ plane (see Fig. 1). By assuming for now the EMT model in (5), it can be shown (see [5, 6, 7] for details) that the wavenumbers must satisfy the dispersion equation

The numerical solution of the nonlinear system in (6b) can be efficiently carried out following the approach described in [19]. The resulting modes are generally hybrid, and contain as special cases the transverse-electric and transverse-magnetic modes that can be supported by isotropic, inductive and capacitive metasurfaces, respectively. Moreover, it is readily verified that the dispersion equation in (6a) is invariant under the duality transformations

The above relationships somehow resemble those observed for self-complementary metasurfaces [46, 47, 48], which exploit the Babinet’s principle. In our case, rather than geometrical self-complementarity, we rely on ${\cal PT}$ symmetry.

In what follows, we will study the effects of the gain-loss interplay in the surface-wave propagation, via the approximate EMT modeling and full-wave numerical simulations.

Looking at the mixing formulae in (5), an interesting question is whether there exist specific combinations of the constituents (in terms of $\sigma_{a}$, $\sigma_{b}$ and $f_{a}$) that render the effective parameters purely reactive. Within the limitations of the EMT, this would imply a perfect balancing of the loss and gain effects, which in the language of ${\cal PT}$-symmetry, corresponds to the symmetric phase [28]. By substituting the complex-valued conductivities (1) in the mixing formulae (5), and enforcing the zeroing of the effective-parameter real parts, we obtain the conditions (see Appendix A for details)

In view of the assumption $\sigma^{\prime}_{a,b}>0$, the condition in (8a) is always feasible ($0\leq f_{a}\leq 1$). On the other hand, we must enforce the additional constraint

in order to ensure that ${\sigma^{\prime\prime}_{a}}$ in (8b) is consistently real-valued. Moreover, in view of our assumption $\sigma^{\prime\prime}_{a}\sigma^{\prime\prime}_{b}>0$, the sign determination in (8b) must be chosen consistently with the sign of $\sigma^{\prime\prime}_{b}$. By enforcing the conditions (8) in (5), after some algebra, we obtain

which identify a class of parameter combinations that yield purely imaginary effective parameters, i.e., loss-gain balance. This result is not necessarily surprising, as a balanced amount of gain and loss may be expected to compensate each other, but we will show hereafter that the gain-loss interplay may have much deeper implications on the modes supported by this metasurface.

By inspection of (10), it can be observed that the conductivity components can differ substantially in the limit $\sigma^{\prime}_{a}\rightarrow\sigma^{\prime}_{b}$. From (8), this yields $f_{a}=f_{b}=0.5$ and (with the sign determination of interest here) ${\sigma^{\prime\prime}_{a}}={\sigma^{\prime\prime}_{b}}$. In other words, by removing irrelevant superscripts, we obtain

which, with a suitable choice of the reference-system origin, corresponds to the aforementioned ${\cal PT}$-symmetry condition

It can be shown (see Appendix A for details) that, under these conditions, the effective parameters reduce to

The remarkably simple expressions in (13) clearly show that the gain-loss interplay can significantly affect the effective-parameter anisotropy. In particular, it is evident that the conditions

would lead to $\left|{\bar{\sigma}}_{xx}\right|\ll\left|{\bar{\sigma}}_{yy}\right|$, i.e., extreme anisotropy. For a quantitative illustration, Fig. 2 shows the anisotropy ratio $\left|{\bar{\sigma}}_{yy}\right|/\left|{\bar{\sigma}}_{xx}\right|$ as a function of the gain-loss parameter ${\bar{\sigma}}^{\prime}\eta$ for various values of the normalized susceptance $\sigma^{\prime\prime}\eta$; note that results do not depend on the sign of $\sigma^{\prime\prime}$, i.e., on the inductive/capacitive character. However, for a given anisotropy ratio, the dispersion characteristics do depend on the reactive character of the metasurface.

For example, Fig. 3 shows some representative dispersion characteristics, in terms of equi-frequency contours (EFCs), for an inductive ($\bar{\sigma}^{\prime\prime}>0$) scenario, by fixing $\sigma^{\prime\prime}\eta=0.1$ and varying the gain-loss parameter $\sigma^{\prime}\eta$. As can be observed, the field exhibits a propagating character within a spectral wavenumber region that can be estimated analytically from (6b) (see Appendix A for details), viz.

with the approximate equality holding in the limit $\sigma^{\prime\prime}\eta\ll 1$. Quite interestingly, this spatial bandwidth is essentially controlled by the normalized susceptance $\sigma^{\prime\prime}\eta$. Outside this region, the field is evanescent, with a purely imaginary $k_{y}$ (not shown for brevity). The gain-loss parameter $\sigma^{\prime}\eta$ controls instead the anisotropy degree. Specifically, starting from the trivial ($\sigma^{\prime}=0$) isotropic case (Fig. 3a), the anisotropy becomes increasingly pronounced by increasing $\sigma^{\prime}$ (Figs. 3b and 3c), and the EFCs approach a limiting curve for $\sigma^{\prime}\eta\gg 1$ (Fig. 3d).

The responses for capacitive ($\bar{\sigma}^{\prime\prime}>0$) scenarios can be in principle obtained via the duality transformations in (7). However, for completeness, they are also exemplified in Fig. 4 directly in terms of the constituent parameters $\sigma^{\prime}$ and $\sigma^{\prime\prime}$. In this case, the propagating spectral region is given by (see Appendix A for details)

where the approximate equality holds in the asymptotic limit $\eta\left(\sigma^{\prime}\right)^{2}\gg\left|\sigma^{\prime\prime}\right|$. We observe that, unlike the inductive case, the spatial bandwidth now depends on both $\sigma^{\prime}$ and $\sigma^{\prime\prime}$; the representative values considered in Fig. 4 are chosen so as to maintain the spatial bandwidth $k_{x}^{(max)}\approx 20k$, progressively moving from perfect isotropy (Fig. 4a) to extreme anisotropy (Fig. 4d).

Canalization effects, intended as the diffractionless transfer of subwavelength features over distances of several wavelengths, have been observed in hyperbolic [18, 19] and extreme-anisotropy [23] metasurfaces, and loss-induced canalization has also been demonstrated [49].

These effects can be intuitively understood by looking at the examples in Figs. 3 and 4 with higher anisotropy (e.g., Figs. 3c, 3d, 4c and 4d). It is apparent that a significant fraction of high-$k_{x}$ spectral components can propagate as unattenuated surface waves in the $x-y$ plane and, in view of the pronounced flatness of the EFCs, their group velocities (normal to the EFCs) are predominantly $y$-directed.

For illustration, Fig. 5 shows the surface-wave propagation along an inductive, non-Hermitian metasurface (with parameters as in Fig. 3c), excited by a $z$-directed elementary dipole. Specifically, Figs. 5a and 5b show the numerically computed (see Appendix B for details) field maps at a close distance from the metasurface, by considering the EMT model and actual conductivity modulation, respectively. Results are in fair agreement, with the differences attributable to nonlocal effects (see Sec. III-E below). In particular, they clearly display the aforementioned canalization effect, with unattenuated, diffractionless propagation of sub-wavelength features. This is markedly different from the conventional cylindrical wavefronts that would be observed in a homogeneous, isotropic case. As a further quantitative evidence, Fig. 5c shows the spatial spectrum (2-D Fourier transform) of the EMT field map, which is essentially peaked around the theoretical EFC. Similar results, not shown for brevity, are observed for capacitive scenarios as well.

Canalization effects like those in Fig. 5 are of great interest for applications to high-resolution imaging. Although the above phenomena may appear qualitatively similar to the canalization effects observed in hyperbolic metasurfaces [18, 19], we stress that the underlying physics is completely different. While in the hyperbolic case these effects are induced by the dual character (inductive/capacitive) of the two constituents, in our case there is no contrast between the reactive parts, and the extreme anisotropy is solely induced by the gain-loss interplay.

To better understand the crucial role played by ${\cal PT}$ symmetry in establishing the extreme-anisotropy response, it is insightful to explore different configurations that do not fulfill such condition. Within this framework, it is also worth highlighting that the unavoidable material dispersion dictates (via causality) that the ${\cal PT}$-symmetry condition may occur only at isolated frequencies [35, 50], and therefore the above phenomena are inherently narrowband.

Figure 6 shows the EFCs pertaining to four representative non-${\cal PT}$-symmetric configurations of interest. Specifically, Figs. 6a and 6b show the (real and imaginary, respectively) results for the inductive parameter configuration in Fig. 3c, but in the absence of gain ($\sigma_{a}\eta=1+i0.1$, $\sigma_{b}\eta=i0.1$). As also evident from the effective parameters (${\bar{\sigma}}_{xx}\eta=0.5+i0.1$, ${\bar{\sigma}}_{yy}\eta=0.019+i0.196$), the extreme anisotropy is now lost and, most importantly, the propagation is severely curtailed in view of the substantial values of $\mbox{Im}\left(k_{y}\right)$. This is clearly visible in the corresponding field map shown in Fig. 7a. Such a stark difference from the ${\cal PT}$-symmetric case suggests interesting possibilities for dynamically reconfiguring the surface-wave response, e.g., switching between a propagating, canalized regime and a strong damping by (de)activating the gain.

Figures 6c and 6d illustrate another interesting example with parameters $\sigma_{a}\eta=0.8+i0.608$, $\sigma_{b}\eta=-1+i0.1$, $f_{a}=0.556$, $f_{b}=0.444$ satisfying the gain-loss balance conditions in (8) but not the ${\cal PT}$-symmetry condition. As expected, since the effective parameters are purely reactive (inductive), the EFCs exhibit propagating and evanescent regions as in Fig. 3. However, although the gain constituent is the same as in Fig. 3c, the spatial bandwidth and the anisotropy are less pronounced. As can be observed in the corresponding field map shown in Fig. 7b, this results in a general deterioration of the canalization effects. Generally speaking, it can be verified numerically that, for a fixed gain constituent, among the infinite parameter configurations that fulfill the gain-loss balance conditions in (8), the ${\cal PT}$-symmetry condition in (12) guarantees the maximum anisotropy.

As previously mentioned, when material dispersion is taken into account, the ${\cal PT}$-symmetry condition can only be perfectly satisfied at isolated frequencies. It is therefore instructive to look at the effects of moderate mismatches that can occur at close-by frequencies. For illustration, Fig. 6e and 6f show the EFCs for parameters as in Fig. 3c when the susceptance are perfectly balanced, but there is a moderate gain-loss mismatch ($\sigma_{a}\eta=1.2+i0.1$, $\sigma_{b}\eta=-1+i0.1$). Interestingly, the resulting effective parameters (${\bar{\sigma}}_{xx}\eta=0.1+i0.1$, ${\bar{\sigma}}_{yy}\eta=-5.95+i6.15$) exhibit simultaneously loss and gain along different directions. Such “indefinite“ non-Hermitian character is also observed in volumetric metamaterials [51, 41]. The EFCs differ substantially from the perfectly ${\cal PT}$-symmetric reference case in Fig. 3c, especially in light of an imaginary part which assumes small negative values at lower wavenumbers and increasingly positive values for higher wavenumbers. Figure 7c shows the corresponding field map, from which we still observe some canalization effects, though with a more visible spreading and attenuation by comparison with Fig. 5.

For the same parameter configuration, Figs. 6g and 6h show instead the effects of the imbalance in the susceptances (while maintaining the inductive character), but assuming now perfect gain-loss balance ($\sigma_{a}=1+i0.1$, $\sigma_{b}=-1+i0.15$). In spite of the sizable imbalance, the EFC departures from the perfectly ${\cal PT}$-symmetric reference case in Fig. 3c seem less dramatic. As can be observed, there is still an extended spectral region characterized by a small positive imaginary part [$\mbox{Im}\left(k_{y}\right)\lesssim 0.02k$] and pronounced anisotropy, although the spatial bandwidth is moderately reduced. Quite interestingly, outside this region we observe $\mbox{Im}\left(k_{y}\right)<0$ which, assuming the conventional choice $\mbox{Re}\left(k_{y}\right)>0$ for the branch-cut, would imply amplification. However, this is not necessarily the case, as in structures mixing gain and loss an a priori choice of the branch-cut is not obvious in the presence of unbounded domains. In fact, a counterintuitive sign flip in the propagation constant has been observed at certain critical incidence conditions [52] (see also the discussion in [41]). In our specific case, the numerical simulations in Fig. 7d indicate the presence of canalization effects accompanied by a mild attenuation, thereby suggesting that $\mbox{Re}\left(k_{y}\right)<0$ may be the proper choice for the branch-cut in the regions where $\mbox{Im}\left(k_{y}\right)<0$.

To sum up, the above examples indicate that, for moderate departures from the ${\cal PT}$-symmetry conditions, canalization effects are still attainable but with reduced resolution and propagation distance. In this respect, the gain-loss (im)balance turns out to be more critical than that in the susceptance.

As previously mentioned, the EMT approximation in (5) is generally accurate for deeply subwavelength modulation periods $d\ll\lambda$. For a more quantitative assessment, Fig. 8 compares its predictions with the rigorous EFCs obtained via full-wave numerical simulations (see Appendix B), for different values of the modulation period $d$. We observe that the results for $d=0.001\lambda$ are hardly distinguishable from the EMT prediction, whereas some visible differences progressively appear for $d=0.01\lambda$ and $d=0.02\lambda$, with changes in the local curvature and moderate reductions in the spatial bandwidth. Interestingly, the propagation constant remain purely real, thereby indicating that the ${\cal PT}$-symmetric-induced gain-loss compensation extends beyond the range of validity of the EMT approximation.

The observed departures from the EMT predictions indicate that nonlocal effects (i.e., spatial dispersion) are no longer negligible. In principle, these effects can be captured by introducing in the effective parameters some wavevector-dependent correction terms, e.g., along the lines of [53].

Although this study is essentially focused on exploring the basic phenomenology, and a practical implementation requires further investigations, one might wonder to what extent the parameter configurations required for ${\cal PT}$-symmetry-induced extreme anisotropy are feasible. Within this framework, the most critical element is the gain constituent, whose implementation varies with the operational frequency of interest. For instance, at microwave frequencies, active metasurfaces typically rely on negative-resistance elements based on amplifiers [54] or tunnel diodes [55]. At terahertz frequencies, an optically pumped graphene monolayer may be a viable option, as it can support population inversion and a negative dynamic conductivity [56]. Specifically, for typical parameters found in the literature [35, 39], the real part can reach values $\sigma_{g}^{\prime}\approx-0.02/\eta$, while the imaginary part $\sigma_{g}^{\prime\prime}$ can be tuned between positive and negative values (ranging approximately from $-0.05/\eta$ to $0.05/\eta$) by acting on the frequency and quasi-Fermi energy. These figures allow in principle to attain the large anisotropy ratios of interest here.

At optical wavelengths, gain media are typically obtained by doping host media with organic dyes [57, 58] or quantum dots [59, 60, 61]. To derive some basic quantitative estimates, we consider the thin-dielectric-layer model in (4), and invert for the complex-valued relative permittivity

Accordingly, the conditions in (14) for extreme anisotropy translate into

For instance, assuming $k\Delta=0.1$, in the capacitive case, a relative permittivity $\varepsilon=1.01-i2$ would yield a normalized conductivity $\sigma\eta=-0.2-i0.001$, like the one considered in the extreme-anisotropy case in Fig. 4d. These permittivity values are in line with those attainable at infrared wavelengths by doping transparent conductive oxides (such as indium tin oxide) with lanthanides [41, 62, 63, 64]. Similar considerations hold for the inductive case too.

As for the practical realizability of the required gain-loss spatial modulation, one possibility could be to rely on high-resolution selective optical pumping [65], possibly based on digital spatial light modulators [66]. Alternatively, one could think of relying on a uniform optical pumping, and patterning a thin layer of gain material with thin, lossy strips, so as to suitably overcompensate the gain in certain selected regions.

The above considerations suggest that the gain-loss configurations of interest are within reach with current or near-future technologies, although further studies are needed to develop some practical designs. In particular, the presence of a substrate should also be taken into account, and a simple EMT modeling may not be applicable, thereby requiring extensive numerical optimization driven by full-wave simulations. These aspects are beyond the scope of the present investigation, and will be the subject of forthcoming studies.

An important issue in non-Hermitian configurations featuring gain is the potential onset of instability, manifested as self-oscillations supported by the system [67]. In previous studies dealing with planar [33] and cylindrical [34] non-Hermitian metasurfaces, the stability issue was addressed by introducing a physical, dispersive model for the gain constituent, and by looking at the poles of a relevant transfer function analytically continued in the complex frequency plane. Here, we follow a similar approach, and consider for the conductivities in the loss and gain regions a Lorentzian (standard and inverted, respectively) model [68]

where $A_{L,G}$, $\omega_{L,G}$ and $\Gamma_{L,G}$ denote some dimensional constants, peak radian frequencies, and damping factors, respectively. Figure 9 shows the above dispersion laws with parameters chosen so as to attain at an operational radian frequency $\omega_{0}$ the nominal values $\sigma_{a}\eta=1+i0.1$, $\sigma_{b}\eta=-1+i0.1$ considered in the canalization example of Fig. 5.

We then consider the EMT dispersion equation in (6b), which can be viewed as a pole condition for the reflection/transmission coefficient under plane-wave illumination from free space, and substitute the dispersive effective parameters computed from (10) with (19). Finally, for fixed values of the wavenumbers, we look for roots of (6b) (i.e., poles) in the complex ${\tilde{\omega}}=\omega^{\prime}+i\omega^{\prime\prime}$ plane. In view of the assumed time-harmonic convention, instability corresponds to complex poles lying in the upper half-plane $\omega^{\prime\prime}>0$.

Figure 10 shows stability maps for representative values of the wavenumbers on the nominal EFC in Fig. 3c. We observe the expected presence of the surface-wave pole on the real axis at the operational radian frequency $\omega^{\prime}=\omega_{0}$, and the absence of poles in the upper half-plane $\omega^{\prime\prime}>0$, which indicates that the system is unconditionally stable for any temporal excitation. The above examples only serve to illustrate that it is possible, in principle, to attain stability via suitable dispersion engineering. Clearly, different parameter choices in the dispersion models and/or the anisotropy ratios may induce the transition of some poles to the upper half-plane $\omega^{\prime\prime}>0$, thereby causing instabilities.