Academia.eduAcademia.edu

Stability of active photonic metasurface pairs

2021, New Journal of Physics

https://doi.org/10.1088/1367-2630/ac37ac

Adding active components to a photonic device may dramatically enrich and improve its performance but, at the same time, creates the risk of instability, namely, occurrence of unwanted self-oscillations. Stability considerations are not always given the attention they deserve when setups employing gain media are investigated; thus, the desired effects or reported regimes may not be achievable. In this work, a generic electromagnetic configuration comprising a pair of planar impedance metasurfaces is examined and analytical stability conditions for its operation are derived. The obtained results for the analyzed basic module can shed light on the stability conditions of more complex active systems that incorporate such components and serve a broad range of applications from imaging and polarization engineering to invisibility cloaking and wavefront transformations.

You may also like PAPER • OPEN ACCESS Stability of active photonic metasurface pairs To cite this article: Constantinos Valagiannopoulos and Sergei A Tretyakov 2021 New J. Phys. 23 113045 View the article online for updates and enhancements. - Review and comparative assessment of engineering correlations for the fireinduced air inflow rate through a compartment opening Dimitrios Pikiokos and Dionysios Kolaitis - Visible spectroscopy of highly charged tungsten ions with the JET charge exchange diagnostic Alex Thorman, Edward Litherland-Smith, Sheena Menmuir et al. - Thermodynamic analysis of Al clusters formation over aluminum melt Alexey Zhokh, Peter Strizhak, Maksym Goryuk et al. This content was downloaded from IP address 178.91.253.97 on 26/11/2021 at 13:52 New J. Phys. 23 (2021) 113045 https://doi.org/10.1088/1367-2630/ac37ac PAPER Stability of active photonic metasurface pairs O P E N AC C E S S Constantinos Valagiannopoulos1 , ∗ R E C E IVE D 14 August 2021 1 2 R E VISE D 10 October 2021 AC C E PTE D FOR PUBL IC ATION 8 November 2021 ∗ and Sergei A Tretyakov2 Department of Physics, Nazarbayev University, Nur-Sultan, KZ-010000, Kazakhstan Department of Electronics and Nanoengineering, Aalto University, FI-00076, Finland Author to whom any correspondence should be addressed. E-mail: [email protected] and [email protected] Keywords: stability, metasurface, active media PUBL ISHE D 26 November 2021 Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. Abstract Adding active components to a photonic device may dramatically enrich and improve its performance but, at the same time, creates the risk of instability, namely, occurrence of unwanted self-oscillations. Stability considerations are not always given the attention they deserve when setups employing gain media are investigated; thus, the desired effects or reported regimes may not be achievable. In this work, a generic electromagnetic configuration comprising a pair of planar impedance metasurfaces is examined and analytical stability conditions for its operation are derived. The obtained results for the analyzed basic module can shed light on the stability conditions of more complex active systems that incorporate such components and serve a broad range of applications from imaging and polarization engineering to invisibility cloaking and wavefront transformations. 1. Introduction Active media, namely substances that are able to pump extra energy into a device, have revolutionized numerous photonic setups by giving them unique utilities via loss compensation, dynamically controlled tunability, and configurability [1]. Indeed, gain can be incorporated in electromagnetic structures to compensate dissipative losses by igniting virtually loss-free operation and amplification [2] or to provide a playground to explore novel topological phases by influencing the corresponding properties of bulk systems [3]. Interesting and important new effects have been reported based on the utilization of active materials like the field enhancement in silicon nanocrystals loaded with quantum dots [4] and surface plasmon boosting by stimulated emission of radiation [5]. Importantly, metamaterials that use particles with gain are able to achieve imaging with increased resolution via negative refraction [6], real-time manipulation of THz waves [7] and gate-controlled switching with graphene [8]. Lasers form a separate category of photonic designs employing active media; for example, lasing can be achieved via injection of ultrashort pulses into gain matter [9] or via electrical drive of individual active nanowires [10]. In particular, plasmonic lasers offer a possibility of exploring extreme interactions between waves and substances [11], enhancing the spontaneous emission rate in strongly coupled nanocavity arrays [12] and realizing spasing action with huge confinement of light [13]. Another class of optical configurations involving active materials is exploiting the so-called parity-time (PT) symmetry concerning a spatial balance between passivity and activity [14]. PT-symmetric configurations, with overall compensation of loss and gain, may offer several counter-intuitive responses like asymmetric light propagation [15], tailored transverse energy flow [16], and optical solitons [17]. All these fascinating effects occurring in active electromagnetic layouts have a hidden enemy: potential instability. Indeed, once gain media are around, the output can increase exponentially even in the absence of any input or for a bounded input. That is why stability analysis accompanies numerous treatises regarding a broad range of phenomena from anomalous scattering [18] and plasmonic resonances [19] to unidirectional cloaking [20], planar focusing [21], and invisibility [22]. Similar stability considerations are used in the framework of devices like saturable active dimers [23], photonic magnifying lenses [24], optically injected semiconductor lasers [25], and injection-locked photonic oscillators [26]. Importantly, © 2021 The Author(s). Published by IOP Publishing Ltd on behalf of the Institute of Physics and Deutsche Physikalische Gesellschaft New J. Phys. 23 (2021) 113045 C Valagiannopoulos and S A Tretyakov Figure 1. (a) The studied physical configuration. A pair of homogeneous metasurfaces characterized by sheet impedances Z1 and Z2 are positioned at distance d apart. They can be excited by both polarizations corresponding to incident plane waves traveling along directions defined by angle θ. The sketched excitation wave contains a TM and a TE part simultaneously. If at least one of the two metasurfaces is active, instabilities may occur. (b) A sketch representing the half plane of the complex frequency ω = Re[ω] + i Im[ω] for Re[ω] > 0. It indicates stability (Im[ω] < 0) or instability (Im[ω] > 0) regions of a photonic system having loss and gain, in relation to the poles of its source-free response. The scenario of lasing (Im[ω] = 0) is also shown. stability of negative non-foster reactive elements in lumped-distributed networks has been examined [27] while the influence of negative capacitors dispersion on the time-evolving response of radiating structures is identified [28]. Finally, as has been proven [29], PT-symmetric structures with exact compensation of gain and loss can be only marginally stable. However, stability considerations are very frequently omitted in many studies dealing with metamaterials-enabled designs, e.g. [30–33]; such a weak point leaves open the question of their own realizability. The objective of this paper is to systematically examine the stability conditions of a basic ubiquitous photonic setup defined by two coupled planar metasurfaces, which has not been yet scrutinized. In general, metasurface configurations are preferred since they may serve a variety of aims while keeping solutions simple. These impedance sheets introduce abrupt changes in optical properties and can be modeled by appropriate boundary conditions or circuit models [34]. Similar paired structures can be used in many applications including tensorial holography [35], reflection control [36], polarization engineering [37], wireless information transfer [38], teleportation of electromagnetic waves [39], waveguiding with tailored dispersion [40], and active cloaking [41]. In this work, we analytically derive and numerically verify the stability conditions of the generic layout formed by two parallel infinite metasurfaces for both field polarizations and any direction of incidence and propagation. The ranges of non-dispersive surface impedances of the two flat metasurfaces that secure conditional stability are determined and the reflectivity variation with respect to the characteristics of the incoming light is identified. A commonly utilized setup of paired sheets is employed [42, 43] and the poles of transfer functions [18, 19] across the complex frequency plane, that lead to instability, are analytically derived. The problem of a bulk slab filled with a single homogeneous medium is also analytically solved, and it is shown that an active layer alone, without a lossy counterpart to absorb the pumped power, inevitably turns out to be unstable. Our results are directly usable for the stability analysis of any integrated photonic system containing this pervasive metasurface component and thus crucial in the extensive assortment of respective optical applications. 2. Problem statement The setup under consideration is depicted in figure 1(a), where the used Cartesian coordinate system (x, y, z) is also defined. It comprises a pair of homogeneous photonic metasurfaces with uniform surface impedances Z1 and Z2 (measured in Ohm), placed into vacuum at a distance d apart. We assume that both metasurfaces are thin sheets maintaining only electric surface current, and that both impedances are real-valued scalars (purely resistive nature, where the resistance can take negative values). This is a generic setup that can be used for amplification of waves in a broad range of photonic devices from laser layouts to PT-symmetric configurations. 2 New J. Phys. 23 (2021) 113045 C Valagiannopoulos and S A Tretyakov If at least one of the two sheet resistances is negative, our active structure can be used to amplify reflected or transmitted waves since the magnitude of the reflection coefficient for illuminating beams can be larger than unity. Such a plane-wave amplifier is stable only if exponentially growing solutions are impossible for all allowed values of the propagation constant along the sheet plane. Let us now regard the operation of such systems as lasers, namely, as sources of plane waves in a simple scenario. In particular, assume that Z2 = 0 while Z1 = −η0 , where η 0 is the free-space impedance; for that selection of parameters, the input impedance seen at the plane of Z1 equals to (−η 0 ) at the frequencies where the distance d is an integral multiple of the operational wavelength. This means that waves traveling in the normal direction, along the axis z, constitute eigensolutions of this system while the reflection coefficient increases without bound at this frequency. Thus, the considered setup can in principle act as a laser, radiating a plane wave in the negative direction of the axis z. However, this functionality can be realized only if such a self-oscillation is stable with respect to small variations of the feeding field amplitude and frequency, meaning that the generation regime is stable. In this paper, we do not consider this stability condition which is determined by the nonlinear properties of the sheets, namely, by the derivatives of the sheet resistances with respect to the field amplitude; it can be studied by conventional methods known for generators, based on negative-resistance elements. Importantly, stability at a given pole at the real frequency axis is not enough for the overall stability; one should additionally ensure that self-oscillations at any other mode are not possible. For usual negative-resistance generators also this condition can be studied by known means, but the regarded system is fundamentally different, because the sheets have an infinite extent across the xy plane, and, thus, excitation of plane waves in any direction (not only along the normal z one) can, in principle, occur. Therefore, we must consider the dispersion equation for plane waves with arbitrary values of the propagation constant in the transverse plane and make sure that, for all other plane-wave modes, exponentially growing solutions are not permitted. Such a consideration is necessary for any other application of metasurfaces containing gain elements. Let us also stress that we consider effectively homogeneous sheets, where the meta-atom size and distances between the elements are much smaller than the oscillation wavelength. Moreover, we examine the situations where the sheet impedance is non-resonant, purely resistive; accordingly, the operational frequency of possible lasing configurations is far from the resonant frequencies of isolated particles. In fact, the lasing frequency is determined not by the resonances of the meta-atoms forming the metasurfaces but by the cavity between the two sheets, similarly to the first lasers and masers. These structures should not be confused with arrays of active resonant particles (meta-atoms), each of which, in the presence of gain, may lase and contribute to extremely complicated collective dynamics, e.g. [44, 45]. Stability of such arrays of resonant meta-atoms needs to be studied by different means than the theory of this paper. 3. Stability conditions Given the isotropy of the impedance sheets, we can consider field solutions that propagate as waves along an arbitrary axis at the sheet plane (we denote this axis as y) and do not depend on the coordinates along its orthogonal direction x. We define the two polarizations with respect to the axis x: the fields with a sole electric component E parallel to the x axis are called TM and the ones with a single magnetic component H parallel to the x axis are called TE. We note that this definition does not restrict the generality of the study as the orientation of the coordinate axes x, y at the isotropic plane is arbitrary. We assume a wave vector k = ky ŷ + kz ẑ constrained to the yz plane and a suppressed harmonic time dependence of the form exp(−iωt). We model the two metasurfaces by their sheet impedances; in particular, we look into the case of resistive sheets, formed by thin layers of conducting materials, which correspond to purely resistive, non-dispersive surface admittances. In terms of the relative permittivity ε of the used material, the sheet η0 , where η 0 , k = ω/c and c are impedance of a layer with thickness h can be approximated by Z = −i (ε−1)kh the impedance, wavenumber and speed of light into vacuum, respectively [46]; such a conjecture is valid when kh ≪ 1 and |ε| ≫ 1. For conductive layers, substituting ε ∼ = 1 + i ωεσ 0 , where σ is the bulk medium conductivity, gives 1/Z ∼ = σh, which is the model of a non-dispersive resistive sheet employed in this work. Note that in case of a dielectric layer with a positive and large real part of ε, the impedance is predominantly reactive and capacitive, while for plasmonic-media layers (negative permittivity behaving as ε∼ = 1 − ωp2 /ω 2 , where ω p is the plasma frequency), the sheet is inductive. In the following, we regard passive or active metasurfaces without reactances (Im[Z1 ] = Im[Z2 ] = 0) in order to serve the dispersion-free conjecture. In other words, we assume that both metasurfaces are thin sheets to maintain only electric surface current, and that both surface impedances are real-valued scalars exhibiting resistive behavior (Z1 , Z2 ∈ R). 3 New J. Phys. 23 (2021) 113045 C Valagiannopoulos and S A Tretyakov With respect to the assumption of negligible dispersion in the frequency range of interest, we note that, according to the Kramers–Kronig relations [47], a frequency-independent real part of passive surface impedance Re[Z] demands a vanishing reactive part Im[Z] = 0. Thus, for passive sheets an assumption of absence of dispersion leads to a conclusion that the sheet must be purely resistive (which is the model that we use in this manuscript). On the other hand, for active elements, the Kramers–Kronig relations do not apply, and the assumption of a non-dispersive negative resistance does not exclude presence of non-dispersive reactance. To the best of our knowledge, no realizations of such exotic properties are known, and we limit our study to purely resistive sheets also in the active case. Our analysis assumes that inevitable dispersive effects take place in the frequency ranges where the active sheet already loses its amplification ability. Presence of resonances at other frequencies does not compromise our conclusions because no instabilities are possible there. To analyze stability and find the self-oscillation conditions, we follow the conventional approach known for bulk-components systems, generalizing it to open, distributed, and multi-modal structures. Starting from linear time-domain differential equations for electromagnetic fields in the absence of sources, we look for a fundamental solution in the Laplace form exp(st), write the eigenvalue equations and study their solutions in the complex plane of s ∈ C. For stable systems, all the roots have negative real parts. To write the characteristic equation, we can use the respective complex phasors of electromagnetic fields and consider the frequency variable in the complex plane, writing s = −iω = Im[ω] − i Re[ω]. In the literature, this method has been extended to distributed systems [28,48], where the connection of an active element with a load via a transmission line is studied by writing the characteristic equation at the beginning of the transmission line. It is important to note that potential instability (exponential growth) of voltage and current at this incoming port (corresponding to a complex frequency value) does not mean that the respective quantities diverge along the transmission line, even though the input impedance seen by the active element is expressed as waves propagating parallel to the transmission line with propagation constant proportional to the frequency. In this work, we use this generalization to account for the open nature of the studied system; however, contrary to the aforementioned earlier studies, our system has infinitely many modes, because waves can propagate in space along any direction. One can state that this is equivalent to infinitely many lines connected to infinitely many active components and thus one should study possible instabilities of all these modes with arbitrary propagation directions to determine the stability conditions. To proceed with that, we consider plane-wave modes in surrounding space whose y-dependence is given by exp(+iβy), where propagation constant is taken real (β ∈ R). It is easy to show (for example, using a developed theory [40]) that for the system of two resistive sheets with arbitrary values of sheet resistances, there are no solutions to Maxwell’s equations at any complex frequency ω with a positive real part (Re[ω > 0]) as long as |β/k| > 1. Thus, we can restrict the analysis to |β/k| < 1 and write β = k sin θ for −90◦ < θ < 90◦ . The conditions for the source-free fields to exist can be written as:   Z2 1+2 cos θ ≡ GTM , η0 (1)    Z1 Z2 e2ikd cos θ = 1 + 2 sec θ 1 + 2 sec θ ≡ GTE , η0 η0 (2) 2ikd cos θ e Z1 = 1+2 cos θ η0  corresponding to the TM and TE waves, respectively. We note that only for dispersionless sheets the right-hand sides of these equations do not depend on the frequency, which allows for fully analytical solutions that we formulate next. Up to this point of derivation, it was assumed that ω > 0 and accordingly β ∈ R; nonetheless, as explained above, to understand the stability of the photonic setup of figure 1(a), we will allow the frequency to take complex values (ω ∈ C), working in the Laplace domain. Note that this analysis excludes possible instabilities in form of inhomogeneous plane waves in space that would exponentially grow along the metasurface plane; however, these solutions cannot be excited in the studied system of homogeneous and infinite sheets. Obviously, instability can only occur in the presence of gain; as noted by figure 1(b), the response increases without bound for Im[ω] > 0 (unstable), dies out for Im[ω] < 0 (stable) and performs undamped oscillations for Im[ω] = 0 (lasing). An alternative way to put it is that stability of generation at a given pole of the real axis is not enough for the stability of the generator. In addition, we must ensure that self-oscillations at any other frequency are not possible; in other words that, for given values of Z1 and Z2 , the transfer function, namely the transmission coefficient, has no poles across the upper half plane of the complex frequency ω, as indicated in the sketch of figure 1(b). 4 New J. Phys. 23 (2021) 113045 C Valagiannopoulos and S A Tretyakov It is straightforward to show that the poles of the system, namely the complex frequencies ω that satisfy the conditions (1) and (2), are given by:      ω 1 2mπ + Arg GTM/TE , (3) Re = ω0 2k0 d cos θ   ω 1 ln |GTM/TE |, (4) Im =− ω0 2k0 d cos θ where ω0 > 0 is a reference frequency used for normalization purposes, k0 = ω0 /c is the associated wavenumber and Arg[z] the primary argument of complex number z ∈ C. The parameter ω 0 is chosen at will and can render the results reported in this study relevant, regardless of the frequency zone that the respective gain media are operated. The periodic waves in the cavity between the two paired metasurfaces create infinite poles since (3) applies for m ∈ Z with Re[ω/ω0 ] > 0. The system is stable when: |GTM/TE |2 > 1. (5) It should be stressed that if the metasurface had a capacitive or inductive component, namely if Im[Z1 ]Im[Z2 ] = 0, we would have frequency-dependent impedances Z1 , Z2 and thus GTM/TE = GTM/TE (ω). Such a feature would render the expressions (3) and (4) together with the conditions (5) null and void. On the contrary, the analytical determination of the correct relations would call for much more complicated algebra if not be totally intractable. Note also that if just one of the two metasurfaces (either Z 1 or Z 2 ) possesses an impedance with value Z = −η 0 sec θ/2 for TM modes or Z = −η 0 cos θ/2 for TE modes, the response of the device is always unstable, no matter what is the value of the other impedance (Z 2 or Z 1 respectively). To elaborate this interesting property, let us consider the case when Z1 = −η0 /2, setting β = 0 (assuming normal incidence, θ = 0) for brevity. We can use an equivalent circuit where the impedance Z2 is connected to a semi-infinite transmission line with the impedance η 0 , modeling the empty half-space across z > 0. On the other side, impedance Z2 is connected to the input impedance of the transmission line of length d with the impedance η 0 loaded at the end (the characteristic impedance of the empty half-space across z < 0) and the impedance of the first sheet, namely Z1 = −η 0 /2. At this special value of Z1 , the input impedance of the considered transmission line does not depend on d and is equal to (−η0 ). We can readily understand that the equivalent circuit contains a loop whose total impedance is zero, since it is a parallel connection of impedances (−η 0 ), Z2 , and η 0 ; thus, the system allows lasing for any value of Z2 . Note that this consideration remains valid regardless of the value of the tangential propagation constant β. This result means that if at least one of the sheet resistances is negative, one can always find such a value of the propagation angle θ for which the system is unstable. Therefore, in realizations it is necessary to limit the possible span of allowed propagation directions; in most cases it is practically done by placing the device inside a waveguide, e.g. [38, 50]. Once the impedances move away from these singular values, the system becomes more stable regardless of the wave propagation direction and, naturally, is always stable once both sheets are passive. However, once they are both active, instability does not necessarily occur; remarkably, when active Z1 , Z2 < 0 become more and more negative beyond the limiting values indicated above, the paired metasurfaces tend to be operated farther from the instability. Furthermore, the conditions (1) and (2) can be used for simpler configurations, too; indeed, if we consider a sheet impedance Z1 = Z above a PEC ground plane, we assign Z2 → 0 in (1) and (2) and find that instability occurs for −sec θ < 2Z/η0 < 0 (TM waves) and −cos θ < 2Z/η0 < 0 (TE waves). The two double inequalities differ from each other because the electric field in TM polarization is always parallel to the metasurfaces while in TE polarization the same quantity is dependent on the incidence angle θ. Since the boundary conditions involve these parts of the field vectors that are tangential to the sheets and the surface impedance Z represents electric response, it is natural the θ-dependence in each stability constraint to be dissimilar. In the one-dimensional scenario (θ = 0), where both polarizations degenerate into a single case, one may consider the (unstable) setup of Z = −η 0 ; the input impedance seen at the plane of Z = Z1 equals to (−η0 ) at the frequencies where kd/π = m − 1/2 for m ∈ N∗ . This means that plane waves traveling in the normal direction, along the axis z, are eigensolutions of this system; see (1) and (2) which become identical for θ = 0. The reflection for normal incidence at these frequencies increases without bound verifying its instability predicted by the derived conditions; such a system can, in principle, act as a laser, radiating a plane wave in the negative direction of axis z. The simplifying assumption of dispersion-free active and passive sheets has allowed us to solve the problem fully analytically. This is because, under this constraint, reactive energy in the system can be stored only in the resonator formed by two parallel and partially reflecting sheets. The resonant condition for the 5 New J. Phys. 23 (2021) 113045 C Valagiannopoulos and S A Tretyakov frequencies of possible instabilities is, in this case, given by the simple formulas (3) and (4), that guarantee the regime of forming a partially standing wave between the two sheets. If we included in the analysis dispersive properties of the sheets, it would be not possible to make any explicit and generic analysis, as only numerical approaches would be feasible. However, our general conclusions would not change. Indeed, the only difference would be that, in addition to the reactive energy stored in the space between the sheets, some extra energy would be stored inside the sheets (modeled by their surface reactances). The frequencies of possible instabilities would correspond to the resonances of the whole system, and the presence of this additional reactance would lead to a shift with respect to the poles values given by (3) and (4). But this does not change the fact that there are infinitely many solutions for these frequencies, since the possibility for formation of a standing-wave mode between the two sheets exists also for dispersive sheets. Thus, we believe that our simplified analysis properly captures the fundamental physical properties of the considered system of two parallel active/passive metasurfaces also for frequency dispersive sheets. To demonstrate the connection of the device behavior with positions of the poles in figure 2, we show the magnitude of the transfer function for two systems and both polarizations on the complex frequency ω plane. Note that in all cases we keep Re[ω] > 0, otherwise the selection of time dependence exp(−iωt) loses its meaning. As predicted from the equations (3) and (4), the poles have the same imaginary part and appear periodically parallel to the real axis with the period π/(k0 d cos θ). In figure 2(a), we examine TM waves for two coupled metasurfaces with |Z1 | < |Z2 | where the first one is active and the second one passive; we directly observe that the system is unstable due to the presence of poles across the upper half plane of complex ω map. In figure 2(b), we investigate the TE case and realize that if the setup is excited only by this type of waves, it will support a stable operation. In figures 2(c) and (d), we keep the same magnitude of impedances but the activity/passivity role is exchanged; it is noteworthy that we record stability for both families of fields since Z2 possesses large negative values (|Z2 | > η 0 ). Let us stress again that the considered stability dynamics concern homogenizable metasurfaces, structured at the subwavelength scale. Since the sheet impedance is assumed to be dispersionless in the frequency range of interest, possible self-oscillation frequencies are determined by the macroscopic resonances of the partially open cavity between the two parallel sheets and not by the resonant frequencies of the individual meta-atoms. 4. Stable system response Since the conditions for stable operation have been thoroughly identified, let us excite the structure by an obliquely incident plane wave under angle θ expressed by Finc = x̂e+ik(z cos θ+y sin θ) producing a reflective wave Fref = x̂Re+ik(−z cos θ+y sin θ) , where F = x̂F is the sole electric field component E for TM waves and the sole magnetic field component H for TE waves. For this excitation scenario one can define the reflectivity ρ = |R|2 for each polarization. Given the reciprocal nature of the considered metasurfaces, the transmissivity will be equal regardless of the side that the system is excited; thus, the reflectivity ρ will be our main metric in the present study. As remarked in sections 2 and 3, the device operation is stable as long as the derived conditions are respected for all β ∈ R and for both polarizations. However, in the following results, we will take these parameters fixed, namely, investigate the conditional stability of the system [51], in the same way that the angular momentum has been taken fixed in a cylindrical analogue [52]. In figure 3, we represent the reflectivity ρTM in dB under TM illumination on the map of the real surface impedances (Z1 /η 0 , Z2 /η0 ) for various incidence angles θ. The blank regions correspond to unstable operation and, thus, there is no meaning in depicting the variation of any response. In figure 3(a), the excitation direction is close to normal, and we notice an infeasible region around the lines 2Z1 /η 0 = 2Z2 /η0 = −sec θ. Importantly, the reflectivity magnitude can be larger than unity due to the presence of an active metasurface and reaches huge levels across the upper left quarter of the map, where the incoming field meets first the active part (Z1 < 0) complimented by the passive part next to it (Z2 > 0). In figure 2(b), we select a more oblique incidence and see that the unstable region stretches slightly but the reflection in the upper left quarter becomes larger. We also observe that the response gets enhanced when Z1 or Z2 is mildly negative combined with a large positive Z2 or Z1 , respectively. In figure 3(c), we further increase the angle θ and thus the conditional instability domain enlarges. It should be pointed out that reflection is low when both components are active (the lower left quarter) and, most importantly, that instability can occur at small values of reflectivity. Such a finding might seem peculiar since instability ignites huge (unbounded) response from the device. However, what we observe in figure 3 is the variation of the reflection magnitude in frequency domain ω, under the assumption that the system converges to a steady-state regime. On the contrary, the decision on the stability is based on the location of the source-free transfer function of the complex frequency plane leading to exponentially increasing or decreasing oscillations with time t. To put it alternatively, the colored areas of figure 3 give the 6 New J. Phys. 23 (2021) 113045 C Valagiannopoulos and S A Tretyakov Figure 2. Magnitude (in dB) of the transfer function (transmission coefficient) of the system depicted in figure 1(a) represented on the complex frequency plane ω/ω 0 , like in the sketch of figure 1(b), where ω 0 > 0 is determined by the application for: (a) TM waves, Z1 /η 0 = −0.5, Z2 /η 0 = +1.5, (b) TE waves, Z1 /η 0 = −0.5, Z2 /η 0 = +1.5, (c) TM waves, Z1 /η 0 = +0.5, Z2 /η 0 = −1.5, (d) TE waves, Z1 /η 0 = +0.5, Z2 /η 0 = −1.5. Plot parameter(s): k0 d = 2, θ = 45◦ . magnitude of time oscillations at a specific time point (say t = 0), not the way they will evolve as the time goes by. In fact, in the presence of white noise, a configuration represented by a non-blank point (Z1 , Z2 ) on √ the maps of figure 3, will produce reflected waves possessing a magnitude ρ that vanish with time. On the other hand, if the point (Z1 , Z2 ) is blank, the reflection will explode with time but at t = 0 will be finite too; however, since the system is unstable, we avoid to show these finite values. Finally, in figure 3(d), we consider a grazing direction (θ = 70◦ ), and a very extensive unstable area is formed while the response from the structure is rather weak, even when the incoming beam illuminates directly the active component. Note that for θ → 90◦ , every single combination of (Z1 , Z2 ) leads to unstable operation and thus rules out unconditional stability for the system; indeed, both polarizations and all excitation directions will be present in thermal noise. Therefore, as also pointed out earlier, the results of figure 3 are presented under the assumption that an efficient filter is employed to block any other, deterministic or stochastic, incoming signal except for that of the specific angle and TM polarization; similar postulations hold for the following graphs. In figure 4, we repeat the same calculations as in figure 3 but for TE-polarized incident plane waves. A basic difference is the size of the infeasible parametric domain which gets shrunk with an increasing angle θ and located, this time, around the cross defined by the singular lines 2Z1 /η 0 = 2Z2 /η 0 = −cos θ. As becomes obvious by inspection of (1) and (2), the stability region for θ = 0 is the same for both polarizations; simply, the worst case of TE waves regarding the stability coincides with the best case for TM waves. The reflectance ρTE takes substantial values for an intermediate angle (θ = 30◦ ) and with similar combinations (Z1 , Z2 ) as in figure 3. Apparently, in the previous (TM) case, the size of the instability region becomes larger and is shaped around a less central point Z1 /η 0 = Z2 /η 0 = −sec θ as θ increases. Note the suppressed reflectance for the oblique incidence scenario of figure 4(d) once both Z1 , Z2 possess a 7 New J. Phys. 23 (2021) 113045 C Valagiannopoulos and S A Tretyakov Figure 3. Reflectivity of TM waves ρTM with respect to the two surface impedances (Z1 /η 0 , Z2 /η 0 ) for various incident angles: (a) θ = 10◦ , (b) θ = 30◦ , (c) θ = 50◦ , (d) θ = 70◦ . Plot parameter(s): k0 d = 2, ω/ω 0 = 1. Blank areas correspond to parametric combinations leading to instability. non-vanishing value, regardless of the sign, a feature that gives an almost symmetric pattern on the (Z1 /η 0 , Z2 /η0 ) map. In figure 5(a), we show the variation of the reflectance ρTM as a function of Z1 /η 0 , while keeping active the second impedance (Z2 /η0 = −2) for various incidence directions θ. If, for a considered parametric combination, the stability constraint (5) is not respected, we disconnect the corresponding data line. One directly notices the peaks for all the curves regardless of θ, when the first metasurface gets slightly lossy Z1 → 0+ , while the response gets weaker for increasing Z1 > 0. Furthermore, as indicated by figure 3, when Z1 becomes more negative beyond the instability limit Z1 /η 0 < −sec θ/2, the reflectivity also decreases on the average. In figure 5(b), we investigate TE waves and the recorded reflection is much stronger compared to figure 5(a). The increasing trend of the device response when approaching the unstable range is again clearly demonstrated. Note that an asymmetric line shape is formed along Z1 /η 0 axis, indicating the ability of the structure to pump more energy when a passive piece is paired with an active one. In figure 5(c), we consider the gain to be fixed (Z1 /η 0 = −2) at the first metasurface and examine the reflectivity ρTM along a continuous range of the other impedance Z2 for several directions of incoming illumination. The magnitude of reflection is similar to that in figure 5(a) but it is shown that the instability interval and the response with respect to Z2 get more sizable with θ. In figure 5(d), we examine TE waves for the same setup of figure 5(c); once again, a larger variation of ρ is noted when the incidence becomes more oblique. Importantly, reflectivity ρTE becomes, on the average, smaller for increasing θ, while the opposite behavior is exhibited by ρTM . To understand further the response of the paired metasurface of figure 1(a), we represent the reflectance at a stable operation regime as a function of the operational frequency ω/ω0 and the incidence angle θ. In this ‘spectra map’ of figure 6(a) we consider TM waves where the first metasurface, directly illuminated by the incident wave, is passive (Z1 > 0). One readily observes a zone of a very weak reflectance that for increasing frequency concerns more oblique excitation; remarkably, this zone is present in all the considered 8 New J. Phys. 23 (2021) 113045 C Valagiannopoulos and S A Tretyakov Figure 4. Reflectivity of TE waves ρTE with respect to the two surface impedances (Z1 /η 0 , Z2 /η 0 ) for various incident angles: (a) θ = 10◦ , (b) θ = 30◦ , (c) θ = 50◦ , (d) θ = 70◦ . Plot parameter(s): k0 d = 2, ω/ω 0 = 1. Blank areas correspond to parametric combinations leading to instability. cases of figure 6, no matter what is the field type or the relative placement of the metasurfaces. Such a feature is related to a perfect matching regime (ρ = 0) under both polarizations [49] calling for the exponential exp(2ikd cos θ) to be real. Thus, given the fact that Z1 , Z2 are also real, parametric loci defined θ by ω/ω0 = mπk0sec for m ∈ N∗ , like the dark zones of figure 6, may lead to zero reflectivity. In other words, d the resonances in figure 6 correspond to destructive Fabry–Perot interference between the two waves into the cavity, for specific distances d that send nothing back to the plane-wave source. Importantly, a resonance appears for ω/ω 0 ∼ = 1.5 close to θ = 60◦ , beyond which the system shifts to instability. In figure 6(b), we examine TE waves, and the reflectance is much weaker as compared to figure 6(a). Despite the presence of active parts, the reflected power is lower than unity (ρTE < 1); note that we use the same magnitude scale for every single scenario of figure 6 for comparison purposes. In figure 6(c), we regard TM fields but this time the primary excitation comes from the other side and meets first the active component (Z1 < 0). Obviously, the reflectivity is enhanced since pumping of energy is encouraged by the incoming wave; however, the position of the resonance is the same as in figure 6(a). Similarly, when the polarization of the wave changes to TE (figure 6(d)), the response is more substantial compared to figure 6(b), but weaker than the TM reflectance of figure 6(c). The reported results are applicable to a wealth of metasurface-based structures involving gain media since paired active–passive sheets are modules used for an extensive range of operations. Examples include invisible sensors [50], wave teleportation devices [38, 39] or active cloaking [41]. Indicatively, the concept of a self-adaptive cloak driven by deep learning with instant response to an ever-changing incident wave and the surrounding environment has been experimentally demonstrated [53]. Furthermore, claddings employing arrays of elementary sources to cancel scattered fields from a hidden object have been fabricated and tested [41]. Stability analyses as those executed in the study at hand can be crucial toward practical realizations of these and similar systems. In addition, measurements on devices achieving an illusion for the 9 New J. Phys. 23 (2021) 113045 C Valagiannopoulos and S A Tretyakov Figure 5. The reflectivity ρ of the metasurface pair of figure 1(a) at several incidence directions θ as a function of: (a) surface impedance Z1 /η 0 , TM waves (Z2 /η 0 = −2), (b) surface impedance Z1 /η 0 , TE waves (Z2 /η 0 = −2), (c) surface impedance Z2 /η 0 , TM waves (Z1 /η 0 = −2) and (d) surface impedance Z2 /η 0 , TE waves (Z1 /η 0 = −2). Plot parameter(s): k0 d = 2, ω = ω 0 . The curves are disconnected once instability occurs. Laplace equation, where the cloaked region is surrounded by controlled sources to protect it from outside detection, is investigated in [54]; this cluster of active elements can be homogenized to formulate a metasurface of the type considered here. 5. The case of a single slab In section 4, we have investigated the response of a model of sheet impedance, assuming that the layer thickness is negligibly small. However, it would be interesting to examine the stability of a similar planar cavity of thickness d, like that of figure 1(a), when it is filled with a homogeneous dielectric of relative complex permittivity ε = n2 , without impedance metasurfaces; that will be a model of a single layer of a significant thickness. As far as the real part of its refractive index is concerned, it will be always positive in the absence of magnetic effects (Re[n] > 0) while we do not consider the case of lossless plasma where Re[n] can be equal to zero. When it comes to the imaginary part, the medium is passive if Im[n] > 0 and active if Im[n] < 0. For such a configuration, the conditions analogous to (1) and (2), involving TM (E x̂) and TE (H x̂) waves respectively, are derived as follows: eikdu cos θ = ± eikdu cos θ = ± where u =  n2 − sin2 θ. 10 1 + u sec θ ≡ ±gTM , 1 − u sec θ n2 + u sec θ ≡ ±gTE , n2 − u sec θ (6) (7) New J. Phys. 23 (2021) 113045 C Valagiannopoulos and S A Tretyakov Figure 6. The reflectivity ρ of the metasurface pair of figure 1(a) with respect to the normalized operational frequency ω/ω 0 and the incidence direction θ for (a) Z1 /η 0 = +1.2 and Z2 /η 0 = −1.5, TM waves, (b) Z1 /η 0 = +1.2 and Z2 /η 0 = −1.5, TE waves, (c) Z1 /η 0 = −1.5 and Z2 /η 0 = +1.2, TM waves and (d) Z1 /η 0 = −1.5 and Z2 /η 0 = +1.2, TE waves. Plot parameter(s): k0 d = 2. As in sections 3 and 4, we again assume that the propagation constant β = k sin θ for −90◦ < θ < 90◦ , considering possible instabilities in form of waves propagating in free space. It is obvious that guided modes along a parallel-plate dielectric slab are always unstable if the slab material is active and, thus, this study does not account for them. Note that the signs of the real and the imaginary parts of u are identical to those of n, since sin θ is real (Re[u] > 0 always, Im[u] > 0 for passive and Im[u] < 0 for active slabs). The transcendental equations (6) and (7) can be solved with respect to the complex frequency ω as follows:   mπ + Arg[gTM/TE ] Re[u] −4 Im[u] ln |gTM/TE | ω Re , = ω0 k0 d|u|2   mπ + Arg[gTM/TE ] Im[u]   +4 Re[u] ln |gTM/TE | ω , Im =− ω0 k0 d|u|2   (8) (9) for m ∈ Z. We again consider only the frequencies with positive real part, thus mπ + Arg[gTM/TE ] > 4 Im[u] ln |gTM/TE |, Re[u] (10) due to the adopted time dependence exp(−iωt). The system is stable once the poles respect Im[ω] < 0 for all permissible m, namely:   mπ + Arg[gTM/TE ] Im[u] > −4 Re[u] ln |gTM/TE |. 11 (11) New J. Phys. 23 (2021) 113045 C Valagiannopoulos and S A Tretyakov Figure 7. Magnitude (in dB) of the transfer function (transmission coefficient) of a single dielectric slab represented on the complex frequency plane ω/ω0 , where ω 0 > 0 is determined by the application for: (a) TM waves with passive medium of ε = 5 + 2i, (b) TE waves with passive medium of ε = 5 + 2i, (c) TM waves with active medium of ε = 3 − i, (d) TE waves with active medium of ε = 3 − i. Plot parameter(s): k0 d = 2, θ = 45◦ . By taking into account that n2 = u2 + sin2 θ, it is straightforward to show that |gTM |2 = 1 + |gTE |2 = 1 + 4 Re[u] cos θ , Im [u] + {Re[u] − cos θ}2 2 4 Re[u] cos θ(|u|2 + sin2 θ) . Im2 [u]cos2 (2θ)   2 + Re[u] − cos θ(|u|2 + sin2 θ) (12) (13) It is, therefore, apparent that |gTM/TE | > 1, meaning that the right-hand side (rhs) of (11) is always negative. To achieve stability, the left-hand side (lhs) of (11) should be larger that its rhs for all the permissible m, as dictated by (10). If the slab is passive (Im[u] > 0), m can take unbounded positive values according to (10), and the minimum value of lhs in (11) is positive; thus, the system is always stable, as expected. Such a feature is demonstrated by figures 7(a) and (b) where the transfer function magnitude of a passive slab (ε = 5 + 2i) is shown across the complex frequency ω plane for TM and TE waves respectively. As well predicted by (8) and (9), the poles (all with Re[ω] < 0) lie across a downward sloping semi-infinite straight line on ω ∈ C map. If the slab is active (Im[u] < 0), m takes integer values below a threshold according to (11). Therefore, there will be always (infinite) poles as m → −∞ for which the stability constraint (11) cannot be respected; as a result, the system is always unstable. It is also indicated by the numerical example of corresponding figures 7(c) and (d), where an active (ε = 3 − i) layer is investigated. This time the poles are located across an upward sloping semi-infinite straight line and thus some (more accurately, infinite) of them are inevitably unstable (Re[ω] > 0). 12 New J. Phys. 23 (2021) 113045 C Valagiannopoulos and S A Tretyakov This result is easy to comprehend by remembering that the solution for fields inside the slab can be found as a sum of plane waves that are multiple times reflected from the two interfaces. If the waves inside the slab partially reflect from the dielectric-air interface at very oblique angles, the reflection coefficient can always become so close to unity so that the sum of the exponentially growing partial waves diverges. One may ask why this instability was not seen in the analysis of two parallel sheets? Indeed, a sheet is a model of a thin material layer; however, the reason is that we have considered purely resistive sheets. Neglecting the reactive nature of the metasurface is tantamount to assuming that there are no guided modes across the sheet. Another significant reason that the active slab is always unstable is that the waves are actually propagating into a gain medium in the absence of any lossy substance in its vicinity. It is also worth mentioning that in the analysis of the slab above, we regarded only the solutions with real |β| < k, even though dielectric films support guided modes with larger tangential propagation constants; if the structure is active, these modes are also unstable yielding exponentially growing fields across the xy plane. 6. Conclusions An often neglected aspect of studies examining photonic designs that incorporate active media, is the stability analysis; thus, there is always the danger of detecting interesting regimes that are unapproachable due to the unbounded increase of the fields, under the inevitable presence of noise. In this work, we have analytically found the conditional stability constraints for a generic setup comprising two parallel planar, potentially active impedance metasurfaces and studied the reflectivity for variable surfaces impedances, oscillation frequencies, and excitation angles. Our findings can be used to check the stability of numerous electromagnetic and optical structures incorporating this generic arrangement of paired sheets designed for controlling the amplitude, direction, and polarization of reflected and transmitted waves. A meaningful expansion of the work at hand would be to take into account the dispersion of impedances and understand how the positions of the poles of the transfer functions on the complex frequency plane are affected. Interestingly, one may also introduce nonlinearities in the response of metasurfaces and observe how the dynamics are modified or if multistability can be achieved, with significant applicability in photonic signal processing and optical memory. Acknowledgments Both authors are grateful to Professor Silvio Hrabar (University of Zagreb, Croatia) for discussions of stability of systems containing transmission lines. The corresponding author (CV) would like also to thank Professor Francesco Monticone (Cornell University, NY, USA) for useful discussions on the stability of planar photonic structures. This work was partially supported by Nazarbayev University Faculty Development Competitive Research Grant No. 021220FD4051 entitled: ‘Optimal design of photonic and quantum metamaterials’. Data availability statement The data that support the findings of this study are available upon reasonable request from the authors. ORCID iDs Constantinos Valagiannopoulos https://orcid.org/0000-0003-1560-2576 Sergei A Tretyakov https://orcid.org/0000-0002-4738-9987 References [1] Boardman A D et al 2011 Active and tunable metamaterials Laser Photon. Rev. 5 287 [2] Hess O, Pendry J B, Maier S A, Oulton R F, Hamm J M and Tsakmakidis K L 2012 Active nanoplasmonic metamaterials Nat. Mater. 11 573 [3] Ota Y, Takata K, Ozawa T, Amo A, Jia Z, Kante B, Notomi M, Arakawa Y and Iwamoto S 2020 Active topological photonics Nanophotonics 9 547 [4] Pavesi L, Dal Negro L, Mazzoleni C, Franzò G and Priolo F 2000 Optical gain in silicon nanocrystals Nature 408 440 [5] Bergman D J and Stockman M I 2003 Surface plasmon amplification by stimulated emission of radiation: quantum generation of coherent surface plasmons in nanosystems Phys. Rev. Lett. 90 027402 [6] Xiao S, Drachev V P, Kildishev A V, Ni X, Chettiar U K, Yuan H-K and Shalaev V M 2010 Loss-free and active optical negative-index metamaterials Nature 466 735 13 New J. Phys. 23 (2021) 113045 C Valagiannopoulos and S A Tretyakov [7] Chen H-T, Padilla W J, Zide J M O, Gossard A C, Taylor A J and Averitt R D 2006 Active terahertz metamaterial devices Nature 444 597 [8] Lee S H et al 2012 Switching terahertz waves with gate-controlled active graphene metamaterials Nat. Mater. 11 936 [9] Brabec T and Krausz F 2000 Intense few-cycle laser fields: frontiers of nonlinear optics Rev. Mod. Phys. 72 545 [10] Duan X, Huang Y, Agarwal R and Lieber C M 2003 Single-nanowire electrically driven lasers Nature 421 241 [11] Oulton R F, Sorger V J, Zentgraf T, Ma R-M, Gladden C, Dai L, Bartal G and Zhang X 2009 Plasmon lasers at deep subwavelength scale Nature 461 629 [12] Zhou W, Dridi M, Suh J Y, Kim C H, Co D T, Wasielewski M R, Schatz G C and Odom T W 2013 Lasing action in strongly coupled plasmonic nanocavity arrays Nat. Nanotechnol. 8 506 [13] Noginov M A et al 2009 Demonstration of a spaser-based nanolaser Nature 460 1110 [14] Rüter C E, Makris K G, El-Ganainy R, Christodoulides D N, Segev M and Kip D 2010 Observation of parity-time symmetry in optics Nat. Phys. 6 192 [15] Feng L, Xu Y-L, Fegadolli W S, Lu M-H, Oliveira J E B, Almeida V R, Chen Y-F and Scherer A 2013 Experimental demonstration of a unidirectional reflectionless parity-time metamaterial at optical frequencies Nat. Mater. 12 108 [16] Peng B et al 2014 Parity-time-symmetric whispering-gallery microcavities Nat. Phys. 10 394 [17] Alexeeva N V, Barashenkov I V, Sukhorukov A A and Kivshar Y S 2012 Optical solitons in PT-symmetric nonlinear couplers with gain and loss Phys. Rev. A 85 063837 [18] Krasnok A, Baranov D, Li H, Miri M-A, Monticone F and Alú A 2019 Anomalies in light scattering Adv. Opt. Photon. 11 892 [19] Lalanne P, Yan W, Vynck K, Sauvan C and Hugonin J-P 2018 Light interaction with photonic and plasmonic resonances Laser Photon. Rev. 12 1700113 [20] Sounas D L, Fleury R and Alù A 2015 Unidirectional cloaking based on metasurfaces with balanced loss and gain Phys. Rev. Appl. 4 014005 [21] Monticone F, Valagiannopoulos C A and Alù A 2016 Parity-time symmetric nonlocal metasurfaces: all-angle negative refraction and volumetric imaging Phys. Rev. X 6 041018 [22] Chen A and Monticone F 2020 Active scattering-cancellation cloaking: broadband invisibility and stability constraints IEEE Trans. Antennas Propag. 68 1655 [23] Zhiyenbayev Y, Kominis Y, Valagiannopoulos C, Kovanis V and Bountis A 2019 Enhanced stability, bistability, and exceptional points in saturable active photonic couplers Phys. Rev. A 100 043834 [24] Savoia S, Valagiannopoulos C A, Monticone F, Castaldi G, Galdi V and Alù A 2017 Magnified imaging based on non-Hermitian nonlocal cylindrical metasurfaces Phys. Rev. B 95 115114 [25] Kovanis V, Gavrielides A, Simpson T B and Liu J M 1995 Instabilities and chaos in optically injected semiconductor lasers Appl. Phys. Lett. 67 2780 [26] Valagiannopoulos C and Kovanis V 2021 Injection-locked photonic oscillators: legacy results and future applications IEEE Antennas Propag. Mag. 63 51 [27] Ugarte-Muñoz E, Hrabar S, Segovia-Vargas D and Kiricenko A 2012 Stability of non-foster reactive elements for use in active metamaterials and antennas IEEE Trans. Antennas Propag. 60 3490 [28] Loncar J, Hrabar S and Muha D 2017 Stability of simple lumped-distributed networks with negative capacitors IEEE Trans. Antennas Propag. 65 390 [29] Lončar J, Vuković J, Krois I and Hrabar S 2021 Stability constraints on practical implementation of parity-time-symmetric electromagnetic systems Photonics 8 56 [30] Valagiannopoulos C A 2016 Optical PT-symmetric counterparts of ordinary metals IEEE J. Sel. Top. Quantum Electron. 22 5000409 [31] Valagiannopoulos C A, Monticone F and Alù A 2016 PT-symmetric planar devices for field transformation and imaging J. Opt. 18 044028 [32] Fleury R, Sounas D L and Alù A 2014 Negative refraction and planar focusing based on parity-time symmetric metasurfaces Phys. Rev. Lett. 113 023903 [33] Anantha Ramakrishna S and Pendry J B 2003 Removal of absorption and increase in resolution in a near-field lens via optical gain Phys. Rev. B 67 201101 [34] Chen H-T, Taylor A J and Yu N 2016 A review of metasurfaces: physics and applications Rep. Prog. Phys. 79 076401 [35] Fong B H, Colburn J S, Ottusch J J, Visher J L and Sievenpiper D F 2010 Scalar and tensor holographic artificial impedance surfaces IEEE Trans. Antennas Propag. 58 3212 [36] Asadchy V S, Díaz-Rubio A and Tretyakov S A 2018 Bianisotropic metasurfaces: physics and applications Nanophotonics 7 1069 [37] Yang Y, Wang W, Moitra P, Kravchenko I I, Briggs D P and Valentine J 2014 Dielectric meta-reflectarray for broadband linear polarization conversion and optical vortex generation Nano Lett. 14 1394 [38] Xiao Z, Ra’di Y, Tretyakov S and Alù A 2019 Microwave tunneling and robust information transfer based on parity-time-symmetric absorber-emitter pairs Research 2019 7108494 [39] Ra’di Y, Sounas D L, Alù A and Tretyakov S A 2016 Parity-time-symmetric teleportation Phys. Rev. B 93 235427 [40] Ma X, Mirmoosa M S and Tretyakov S A 2020 Parallel-plate waveguides formed by penetrable metasurfaces IEEE Trans. Antennas Propag. 68 1773 [41] Selvanayagam M and Eleftheriades G V 2013 Experimental demonstration of active electromagnetic cloaking Phys. Rev. X 3 041011 [42] Ferreira A, Viana-Gomes J, Nilsson J, Mucciolo E R, Peres N M R and Castro Neto A H 2011 Unified description of the dc conductivity of monolayer and bilayer graphene at finite densities based on resonant scatterers Phys. Rev. B 83 165402 [43] Sarsen A and Valagiannopoulos C 2019 Robust polarization twist by pairs of multilayers with tilted optical axes Phys. Rev. B 99 115304 [44] Qiao X et al 2021 Higher-dimensional supersymmetric microlaser arrays Science 372 403 [45] Midya B, Zhao H, Qiao X, Miao P, Walasik W, Zhang Z, Litchinitser N M and Feng L 2019 Supersymmetric microring laser arrays Photon. Res. 7 363 [46] Tretyakov S 2003 Analytical Modeling in Applied Electromagnetics (Boston, MA: Artech House Publishers) [47] Landau L D and Lifshitz E M 1984 Electrodynamics of Continuous Media 2nd edn (Oxford: Pergamon) [48] Lončar J, Muha D and Hrabar S 2015 Influence of transmission line on stability of networks containing ideal negative capacitors IEEE Int. Symp. on Antennas and Propagation USNC/URSI National Radio Science Meeting (Vancouver, Canada) pp 73–4 14 New J. Phys. 23 (2021) 113045 C Valagiannopoulos and S A Tretyakov [49] Valagiannopoulos C 2020 Analytical inverse design of polarization-insensitive photonic filters by tailoring Fabry–Perot interference IEEE Access 8 141641 [50] Fleury R, Sounas D and Alù A 2015 An invisible acoustic sensor based on parity-time symmetry Nat. Commun. 6 5905 [51] Oppenheim A V, Willsky A S and Nawad S H 1997 Signals and Systems (Englewood Cliffs, NJ: Prentice Hall) [52] Yerezhep B and Valagiannopoulos C 2021 Approximate stability dynamics of concentric cylindrical metasurfaces IEEE Trans. Antennas Propag. 69 5716 [53] Qian C, Zheng B, Shen Y, Jing L, Li E, Shen L and Chen H 2020 Deep-learning-enabled self-adaptive microwave cloak without human intervention Nat. Photonics 14 383 [54] Ma Q, Mei Z L, Zhu S K, Jin T Y and Cui T J 2013 Experiments on active cloaking and illusion for Laplace equation Phys. Rev. Lett. 111 173901 15