Stability of active photonic metasurface pairs

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.


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, 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][31][32][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.

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 Z 1 and Z 2 (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.
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 Z 2 = 0 while Z 1 = −η 0 , where η 0 is the free-space impedance; for that selection of parameters, the input impedance seen at the plane of Z 1 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.

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 = k yŷ + k zẑ 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 impedance of a layer with thickness h can be approximated by Z = −i η 0 (ε−1)kh , where η 0 , k = ω/c and c are 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 where ω p is the plasma frequency), the sheet is inductive. In the following, we regard passive or active metasurfaces without reactances (Im[Z 1 ] = Im[Z 2 ] = 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 (Z 1 , Z 2 ∈ R).
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 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: 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 Z 1 and Z 2 , 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).
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: Re Im where ω 0 > 0 is a reference frequency used for normalization purposes, k 0 = ω 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:  (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 Z 1 = −η 0 /2, setting β = 0 (assuming normal incidence, θ = 0) for brevity. We can use an equivalent circuit where the impedance Z 2 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 Z 2 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 Z 1 = −η 0 /2. At this special value of Z 1 , 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 ), Z 2 , and η 0 ; thus, the system allows lasing for any value of Z 2 . 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 Z 1 , Z 2 < 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 Z 1 = Z above a PEC ground plane, we assign Z 2 → 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 = Z 1 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 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 π/(k 0 d cos θ). In figure 2(a), we examine TM waves for two coupled metasurfaces with |Z 1 | < |Z 2 | 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 Z 2 possesses large negative values (|Z 2 | > η 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.

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 F inc =xe +ik(z cos θ+y sin θ) producing a reflective wave F ref =xRe +ik(−z cos θ+y sin θ) , where F =xF 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 (Z 1 /η 0 , Z 2 /η 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 2Z 1 /η 0 = 2Z 2 /η 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 (Z 1 < 0) complimented by the passive part next to it (Z 2 > 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 Z 1 or Z 2 is mildly negative combined with a large positive Z 2 or Z 1 , 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 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 (Z 1 , Z 2 ) on the maps of figure 3, will produce reflected waves possessing a magnitude √ ρ that vanish with time. On the other hand, if the point (Z 1 , Z 2 ) 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 (Z 1 , Z 2 ) 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 2Z 1 /η 0 = 2Z 2 /η 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 (Z 1 , Z 2 ) 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 Z 1 /η 0 = Z 2 /η 0 = −sec θ as θ increases. Note the suppressed reflectance for the oblique incidence scenario of figure 4(d) once both Z 1 , Z 2 possess a non-vanishing value, regardless of the sign, a feature that gives an almost symmetric pattern on the (Z 1 /η 0 , Z 2 /η 0 ) map.
In figure 5(a), we show the variation of the reflectance ρ TM as a function of Z 1 /η 0 , while keeping active the second impedance (Z 2 /η 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 Z 1 → 0 + , while the response gets weaker for increasing Z 1 > 0. Furthermore, as indicated by figure 3, when Z 1 becomes more negative beyond the instability limit Z 1 /η 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 Z 1 /η 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 (Z 1 /η 0 = −2) at the first metasurface and examine the reflectivity ρ TM along a continuous range of the other impedance Z 2 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 Z 2 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 (Z 1 > 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 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 Z 1 , Z 2 are also real, parametric loci defined by ω/ω 0 = mπ sec θ k 0 d for m ∈ N * , like the dark zones of figure 6, may lead to zero reflectivity. In other words, 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 (Z 1 < 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 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.

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 ε = n 2 , 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: e ikdu cos θ = ± n 2 + u sec θ n 2 − u sec θ ≡ ±g TE , where u = n 2 − sin 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: Re Im for m ∈ Z. We again consider only the frequencies with positive real part, thus due to the adopted time dependence exp(−iωt). The system is stable once the poles respect Im[ω] < 0 for all permissible m, namely: By taking into account that n 2 = u 2 + sin 2 θ, it is straightforward to show that It is, therefore, apparent that |g TM/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). 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.

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.