Free-Space Nonreciprocal Transmission Based on Nonlinear Coupled Fano Metasurfaces

: Optical nonlinearities can enable unusual light–matter interactions, with functionalities that would be otherwise inaccessible relying only on linear phenomena. Recently, several studies have harnessed the role of optical nonlinearities to implement nonreciprocal optical devices that do not require an external bias breaking time-reversal symmetry. In this work, we explore the design of a metasurface embedding Kerr nonlinearities to break reciprocity for free-space propagation, requiring limited power levels. After deriving the general design principles, we demonstrate an all-dielectric ﬂat metasurface made of coupled nonlinear Fano silicon resonant layers realizing large asymmetry in optical transmission at telecommunication frequencies. We show that the metrics of our design can go beyond the fundamental limitations on nonreciprocity for nonlinear optical devices based on a single resonance, as dictated by time-reversal symmetry considerations. Our work may shed light on the design of ﬂat subwavelength free-space nonreciprocal metasurface switches for pulsed operation which are easy to fabricate, fully passive, and require low operation power. Our simulated devices demonstrate a transmission ratio >50 dB for oppositely propagating waves, an operational bandwidth exceeding 600 GHz, and an insertion loss of <0.04 dB.


Introduction
Optical nonreciprocity allows asymmetric light transmission from the opposite sides of a two-port device [1,2]. A common approach to break reciprocity is to utilize a magnetic bias [3,4], or in general a bias with odd symmetry under time-reversal, such that the propagating signals from opposite sides see entirely different devices. While a DC magnetic bias is bulky and difficult to integrate along with other components [3,4], the use of artificial angular or linear momentum bias achieved by periodically modulating the material parameters has enabled a paradigm shift to realize nonreciprocal components [5]. Time modulation has been successfully implemented on different platforms, from acoustics [6] to RF [7,8], THz [9], and optical frequencies [5,10],. Free-space nonreciprocal metasurfaces based on time modulation have been explored spanning the RF frequency range [11][12][13] and in the visible frequency range [14]. However, modulating the material parameters comes at the cost of power consumption, and it generates spurious harmonics that may interfere with nearby components [15], unless the devices are judiciously designed [16]. In addition, time modulation becomes harder as the operating frequency increases because the modulation frequency should scale accordingly [15]. Breaking the passivity of the system through the use of unidirectional gain elements (amplifiers) has also been explored to break reciprocity, leading to small size, ease of fabrication, and integrable CMOS technology for nonreciprocal responses [17][18][19][20]. However, the main disadvantage of such approaches is the large power consumption and noise, because amplifiers are active elements that need to be properly biased in order to operate [20]. Additionally, the use of a CMOS amplifier limits the dynamic range of operation, so these metasurfaces are limited to work in the low power regime not exceeding a few dBm in order to avoid triggering nonlinearities in the amplifiers.
An alternative to these techniques, with the appealing property of not requiring any form of external bias or power, is the use of nonlinearities [21] combined with geometrical asymmetries. Light waves traveling through these devices can modulate the material refractive index while passing through it [22], and the geometrical asymmetry ensures that the induced variations are different for opposite propagation directions, breaking reciprocity. Due to the typically weak optical nonlinearity of nonlinear materials, these selfbiased nonreciprocal devices have to rely on strong optical resonances and large contrast in the induced field distributions [23][24][25].
An ideal two-port optical isolator supports unitary transmission in one direction and full absorption (zero transmission) in the opposite direction, where the deviation from unitary transmission in the forward direction defines the insertion loss and the finite transmission in the reverse direction defines the isolation of the device [26]. In reality, nonlinearity-based nonreciprocal devices do not work like conventional isolators, in that they do not absorb the backward power, but instead reflect it [23,25,27]. While it is strictly forbidden by thermodynamic considerations to realize a passive two-port isolator using a lossless material, as we consider throughout our work [28,29], we emphasize that nonlinearity-based nonreciprocal devices are not true isolators, as they cannot isolate a weak backward-propagating signal in the presence of a strong forward signal [29], as dictated by dynamic reciprocity [30]. Therefore, our nonreciprocal device only ensures asymmetric transmission features in a switched mode operation, based on which it is separately excited at the two ports. Thermodynamic considerations [31][32][33] also imply that asymmetric transmission cannot arise for all power levels, as is consistent with the fact that nonreciprocity is inherently rooted in the nonlinear response of the system. Different metrics can be used to characterize the performance of the devices under consideration. We can define a Nonreciprocal Intensity Range (NRIR) [5,27] as the range of input intensities for which we expect to find a largely different transmission level for opposite excitations. Interestingly, a nonlinearity-based nonreciprocal device based on a single nonlinear resonator, with an arbitrary Fano lineshape, suffers from an important trade-off between the minimum insertion loss in the system and its NRIR [27]. These limitations, along with the high power requirements and the associated signal distortions, hinder several applications of these concepts, despite their appealing bias-free, fully passive nature. Only recently, a few proposals to realize free-space nonreciprocal devices based on these concepts have been presented, based on coupled dielectric spheres in multilayers [23] or narrow slits in a dielectric grating [24,27], in addition to proposals based on guided wave implementations in optical [1] and RF frequencies [34]. In this work, we show theoretically and with full-wave simulations that, by using nonlinear metasurface bilayers separated by a subwavelength spacer, we can obtain a nonreciprocal compact device that overcomes the trade-off between insertion loss and NRIR in a free-space implementation. Moreover, we discuss the design guidelines and the response of the device at different wavelengths in order to assess its spectral bandwidth. Figure 1a schematically depicts the functionality of the nonlinear metasurface bilayer described in this work; a few designs will be explored in the following sections to demonstrate both high and low power operations. We use silicon (Si) as the material of choice, given its wide availability, compatibility with CMOS fabrication processes, and relatively large third-order nonlinear susceptibility χ (3) ≈ 2.8 × 10 −18 m 2 /V 2 [22]. The choice of Si is also motivated by the possibility of comparing our results with recent reports using silicon metasurfaces for similar purposes [24,27]. We use glass as the spacer in our design due to similar reasons, apart from the fact that it has no appreciable optical nonlinearity in the frequency and intensity ranges of interest. The device operation in Figure 1a can be explained in analogy with the microwave circuits introduced in [34], but with the Fano resonance realized using dielectric metasurfaces instead of electrical circuits. Fano resonant metasurfaces have been widely studied in nanophotonics. More recently, there has been significant work devoted to realize high-quality factor Fano resonances for optical and visible frequencies based on bound states in the continuum (BICs) [35][36][37][38][39][40][41]. Due to their enhanced light-matter interactions, these BIC metasurfaces can be efficiently employed in nonlinear optics applications [35]. (a) Schematic depicts the nonlinearity-based nonreciprocal metasurface bilayer presented in this work. The arrows represent high power plane waves pointing to the direction of propagation that are blocked and reflected when propagating from the top to the bottom (red arrow) and transmitted when propagating from the bottom to the top (blue arrows). The zoomed in box shows one unit cell of the metasurface top layer made of Si (the bottom layer is similar, but with different geometrical parameters). Upon the plane wave excitation, the transmission from each layer takes the Fano line shape in the dashed box, and can be explained by the interference between: (1) the dark narrow band mode similar to the waveguide mode inside the unperturbed slab (i.e., with filled Si in the dashed box) and (2) the bright wideband mode that resembles a Fabry-Pérot mode or background reflection, as shown in the top and bottom of the shadowed rectangle, respectively. (ω 1 , ω 2 ) are the resonance frequencies of the modes, while their complex amplitudes are (a 1 , a 2 ) and their coupling rates to the free space are (κ 1 , κ 2 ) for the dark and bright mode, respectively. (b) Transmission in the forward direction vs the nonreciprocal intensity range for the various nonlinear resonator designs presented in this work (red asterisks). The blue shaded region corresponds to the bound in Equation (5).

Principle of Operation of Nonlinearity-Based Nonreciprocal Devices Based on Coupled Fano Metasurfaces
To start, we designed the top and bottom silicon surfaces in Figure 1a such that each layer supports a Fano resonance in the frequency range of interest based on the guided mode resonances concept (GMR) [24]. The zoomed-in box in Figure 1a shows one unit cell of the top metasurface layer. The Fano resonance of each layer can be explained similarly; here, we focus only on the top layer, of which the resonance can be explained by considering the interference between the broadband (bright) and narrowband (dark) modes. The (dark, bright) mode has a complex amplitude (a 1 , a 2 ), resonance frequency (ω 1 , ω 2 ) and linewidth (κ 1 , κ 2 ) defining its coupling to free space. Here, these two modes correspond to (1) the waveguide mode that can be coupled to free-space radiation through the periodic perturbation along the slab interface, and (2) the Fabry-Pérot mode of the Si slab, or the broadband reflection from this slab, as shown in the top right and bottom right panels of Figure 1a, respectively. In general, κ 2 κ 1 because the bright mode has a much larger coupling rate to the radiation continuum than the dark mode has. Therefore, when excited by plane waves, these modes interfere, creating a Fano line shape in the transmission curve due to their different coupling rates, as shown in the dashed box in Figure 1a. We define ω 0 as the resonance frequency of the Fano resonator, at which the linear transmission coefficient goes to zero. This frequency also equals the resonance frequency of the dark mode, ω 0 = ω 1 , as shown in the dashed box in Figure 1a. In particular, when we excite the Fano resonator with an incident plane wave with frequency ω, the transmission curve takes the form (see Equation (8) in Section 5.1 (Coupled Mode Theory)) where δ 1 = ω − ω 1 and δ 2 = ω − ω 2 . It is readily seen that, when δ 1 or δ 2 = 0, we have zero transmission; therefore, ω 1 = ω 0 in a linear case, and when κ 1 δ 2 = −κ 2 δ 1 we have unitary transmission. Additionally, because we assume κ 2 κ 1 because it describes the bright mode, it follows that the mode frequency ω 2 ω 1 to obtain unitary transmission is very close to the zero transmission and we have appreciable Fano resonance, as shown in the dashed box in Figure 1a.
Next, we now include the Kerr nonlinearity of Si, so that the refractive index n changes as a function of the incident field intensity according to the relation n = n 0 + n 2 I, where n 0 is the refractive index at low incident power and n 2 characterizes the Kerr nonlinearity according to the relation n 2 = 3 µ 0 / 0 χ (3) / 4n 2 0 > 0. When we have high incident power, the resonance frequency ω 1 is no longer constant, i.e., ω 1 = ω 0 ; instead, it depends on the frequency and power of the incident wave, such that (see Equation (11) in Section 5.2 (Nonlinear Bistability)) where δ 01 = ω − ω 0 is the frequency shift where we recall that ω 0 is the resonance frequency of mode 1 before including the nonlinearity, and P is a normalized term proportional to the incident power. Because we assume that the mode frequency ω 2 is very far away from the frequency range of interest, its properties are weakly affected by the high power, so we assume its parameters ω 2 and κ 2 to be fixed in the nonlinear analysis. In the following, we study optical switching in nonlinear Fano resonators near their resonance frequencies as a function of the incident power P and input frequency ω. First, we assume an incident monochromatic wave with fixed normalized power (P = 0.5), and we study the response of the nonlinear Fano resonance as a function of the input frequency ω. Inspecting Equation (2), we can immediately recognize that a regime can arise for which the response of the nonlinear Fano resonator becomes bistable, a scenario of particular importance for the findings described in this paper. In particular, we expect a bistable response at some threshold input power if (see Section 5.3 Bistability Condition)) where In Figure 2a, we show the resonance frequency shift δ 1 as a function of the applied incident frequency for P = 0.5. Generally, the resonance frequency shift increases with the incident frequency; however, in a range of frequencies slightly below the resonance frequency for low intensities ω 0 , it abruptly jumps to the upper branch of the bistability curve, as indicated by the black arrow at the critical point c 1 . The further increase of the incident frequency leads to a monotonic increase of the resonance shift. If we now consider a decreasing input frequency, we find an abrupt jump of the resonance shift from the upper branch to the lower branch at the critical point c 2 , implying a hysteresis in the bistable response. Interestingly, for Fano resonators, we can choose the parameters (κ 1 , κ 2 , ω 1 , ω 2 ) such that this discontinuous jump enables zero to unitary switching in transmission, by requiring δ 1 = 0 at the lower branch, while κ 2 δ 2 = −κ 1 δ 1 at the upper branch. This particular scenario is identified by the horizontal dashed line in Figure 2a, where we choose the parameters to ensure δ 1 = 0 at the point c 1 on the lower branch, and δ 1 κ 2 = −δ 2 κ 1 on the upper branch. This interesting response may find potential applications for frequency filters [22]. In order to demonstrate the abrupt zero to the unitary transmission feature, we plot the transmission of this Fano resonator as a function of the incident frequency in Figure 2b, where it indeed shows that the point c 1 offers very low transmission and, as we increase the frequency beyond a threshold value, the transmission abruptly switches to unity (0 dB). In general, we can design the bistability region as desired, e.g., to have the abrupt change in transmission at some desired frequency; however, we emphasize that, because we assumed n 2 > 0, the bistability can only happen for input frequencies smaller than the resonance frequency at low intensities. We can prove this important property by solving Equation (2) exactly such that δ 1 has a single real solution, or we can graphically compare the bistable transmission curve to the linear transmission curve, i.e., the solid blue line in Figure 2b, showing that we can find a bistability only when the incident frequency is lower than the linear resonance frequency, as in the case of nonlinear Lorentzian resonators [42]. This behavior can be observed for any desired power level P, as shown in Figure 2c. A larger incident power generally results in a lower switch-up frequency threshold and a wider hysteresis of the bistable response, and in the limit of very low power levels the transmission converges to the linear scenario, confirming again that the bistability can only occur for incident frequencies lower than the resonance frequency. The opposite remarks apply for n 2 < 0. Figure 2c confirms that, for P = 0.5, we find an abrupt change in transmission from 0 to 1, corresponding to the dashed line in Figure 2a. Next, we study the nonlinear Fano resonator properties for fixed incident frequency while changing the input power P. Similar conclusions can be drawn by observing the results in Figure 2d-f. For instance, Figure 2d shows the nonlinear shift δ 1 as a function of the incident power P while keeping the frequency fixed (ω = 0.995ω 0 ). Similar to Figure 2a, we select the parameters that enable an abrupt jump in transmission from 0 (−40 dB) to 1 (0 dB) as we increase the incident power. This can be confirmed in the transmission plot in Figure 2e, in which at the critical point c 1 the lower branch shows 0 transmission, and a slight increase in power causes an abrupt jump to unitary transmission. This behavior has potential in optical applications like optical limiters, nanoswitches [22,43] and nonreciprocal devices [23]. Similar to Figure 2c, we plot in Figure 2f the transmission as the input power increases for different frequencies. We observe again that the bistability only arises when the input frequency is lower than the resonance frequency of the linear resonator. In addition, as the incident frequency becomes smaller, the power needed to trigger a bistable transition becomes larger. Because this bistability regime can be used to induce a nonreciprocal response, by making sure that the system operates in different stable states for opposite propagation directions, we deduce that a low-power nonreciprocal device can be achieved if the operating frequency is close to, and less than, the resonance frequency of the linear Fano resonator ( Figure 2c).
Having analyzed the response of a single Si layer operated as a Fano resonator, we can now consider nonlinear optical devices based on coupled nonlinear Fano resonators, as shown in Figure 1a, where each Si layer supports a Fano resonance with its own resonance frequency and decay rate, and they are separated from each other by an electrical length θ. A functional working principle of a nonreciprocal device based on this geometry consists in the following scheme: for a given excitation frequency and power, one of the resonators (resonator 2) is designed to enter its bistable region, while the other one (resonator 1) remains outside this region. In order to distinguish the dark mode resonance frequency of each Fano resonator, we call the dark mode of resonator 1 ω 01 , while the dark mode of resonator 2 is ω 02 , and we recall that ω 0i is the frequency at which the linear transmission goes to zero. Therefore, we can design resonator 2 with the linear resonance frequency ω 02 > ω in > ω 01 . The linear spectral transmission curves for resonator 1 and resonator 2 are shown in Figure 3a for the parameters in Table 3 in the first and second rows in Section 5.5 (Fano Resonator Parameters) for resonator 1 and resonator 2, respectively. We now consider an incident excitation with increasing power P and of frequency ω in indicated by the vertical dashed line in Figure 3a, such that ω 01 < ω in < ω 02 . As we increase the power of the incident wave, resonator 2 exhibits a sharp transition in the transmission curve at the critical point P = 0.7, confirming that it operates in the bistable region, while resonator 1 maintains a smooth variation of transmission, as shown in Figure 3b. For operation as an ideal nonreciprocal device, it is particularly important to design the two resonators to have unitary transmission at the same power level [31], as shown in Figure 3b for the intensity around P = 1. By placing the resonators separated by an electrical distance θ, as shown in the inset of Figure 3c, we can test the working principle of this device by exploring the response for excitation from opposite sides. First, we excite the structure with frequency ω in from the side of resonator 2. Initially, with the very small incident power P, the reflection is large and, as we increase the power, we reach the critical point at which the transmission for resonator 2 goes from 0 to 1; therefore, the signals are transmitted from resonator 2 to resonator 1. Because we designed resonator 1 to support unitary transmission at the same critical power as resonator 2, the incident signal is transmitted to the left port. Overall, the transmission as a function of intensity for excitation from right to left T RL is similar to the transmission of the individual resonator 2, as shown by the dashed line in Figure 3c. Next, we consider the excitation from the left side. Again, we start with a small input power, for which the incident wave is partially transmitted to resonator 2. Because at this power level resonator 2 is highly reflective, the wave experiences multiple reflections between the resonators, and it reaches a steady-state transmission level T a . The incident power on resonator 2, P eff , can be calculated as (see Equation (18) in Section 5.4) where T a is the nonlinear transmission for excitation from the side of resonator 1, and φ a and φ b are the nonlinear reflection phases from resonator 1 and resonator 2, respectively, as shown in Figure 3e, which are not to be confused with the reflection phase in the linear regime shown in Figure 3d. While the former shows an abrupt change in phase at the bistability transition, the latter obviously does not. We used the parameters for resonator 1 and resonator 2 as given in Table 3 (row 1 and row 2) in Section 5.5 (Fano Resonator Parameters).
The solution of Equation (4) considers an effective value of T a stemming from the multiple reflections between the two resonators, which may be calculated based on numerical simulations, but it can also be evaluated with good approximation assuming that T a is only affected by the incident power from the left, ignoring the reflected power from resonator 2. In other words, when the wave experiences the first round trip, resonator 1 will see two incoming waves: the incident wave from left, and the reflected wave from resonator 2 from the right. Therefore, T a should be updated based on the total input power from left and right. This process should repeat; in each iteration we update T a based on the total incoming power from left and right until we reach a steady state; however, in Equation (4) we assumed that T a is updated only once, neglecting that the reflection from resonator 2 will update the transmission T a . This approximation makes it easier to develop physical insights into the response of the system; however, as shown in the numerical results in the next section, this approximation is quite accurate. We calculate the value of P eff for different electrical lengths θ as a function of the incident power P, as shown in Figure 3f. It is evident that P eff can be larger or smaller than P, thereby shifting the transmission curve of resonator 2 to higher or lower input power, respectively. In particular, for 2θ = 2.1π, P should be 0.4 < 0.7 in order to obtain a bistable transition in resonator 2 which is lower than the transition power when excited from the right; in contrast, 2θ = 1.5π shows that P should be 0.95 > 0.7 in order to enable the bistable transition. These particular values are indicated by the vertical dashed lines in Figure 3f. For the two values of θ considered, we plot the transmission from left to right T LR , as shown in Figure 3c. Thus, we obtain a nonreciprocal response with a large transmission from left to right, and zero transmission from right to left for a specific value of θ. We stress that this is not equivalent to a conventional isolator, as the energy from right to left is not absorbed in the device, but rather it is reflected. Interestingly, as shown in Figure 3c, we can obtain an arbitrarily large forward transmission~100% and complete reflection~100% for a given range of incident power. We define the ratio between the power levels at the edges of the isolation region as the nonreciprocal intensity range (NRIR). Time-reversal symmetry implies a fundamental limit to NRIR for devices based on a single nonlinear Fano resonator with asymmetric decay rates to the two ports [27].
where T is the forward transmission. This is due to the fact that time-reversal symmetry requires the transmission to be identical, even in largely asymmetric resonators, if they fully transmit the input signal from one port [44]. This trade-off is plotted in Figure 1b, where the shaded blue region denotes the range of transmission and nonreciprocal intensity ranges available for any nonlinearity-based nonreciprocal device based on a single resonance satisfying the inequality (4). Interestingly, in our case, similar to previously demonstrated nonlinearity-based nonreciprocal devices based on coupled resonators [34], we do not need to comply with this limitation, and we can have an arbitrarily large forward transmission over a wide range of incident power levels. In particular, we show in our numerical results that some of the proposed designs (red asterisks in Figure 1b) can overcome the bound in Equation (5) (blue dashed region). Our approach is consistent with the method used in [34], with the difference that, in our case, we use two coupled Fano resonators, instead of a coupled Fano and Lorentzian resonator. This approach makes it easier to design the two resonators so that they have the same geometric topology with only slight deviations in the geometrical parameters for the different resonators, a property that is important for optical components where fabrication errors can easily occur and drift the parameters of the two resonators in a correlated manner. In addition, we found that it is difficult to design Fano and Lorentzian resonators similar to the approach in [34] with similar quality factors based on metasurfaces. For instance, in one of the designs we will present, the quality factor of the Fano resonance based on guided mode resonance is in the order of~10 3 , attained using only a thin Si layer of a few nanometers thickness~λ/6; if we try to design a Lorentzian resonant metasurface with the same quality factor based on, e.g., a Fabry-Pérot resonance, we need a homogeneous slab of Si with thickness ∼ 10 3 n × λ, where λ is the operation wavelength, meaning that the device will be extremely thick. Finally, we notice that the power factor P scales with the nonlinear coefficient, P ∝ 1/ χ (3) 2 ; therefore, in order to obtain both a lower operating power and larger bandwidth, we can use materials with larger χ (3) , for example multiple quantum wells, which have been shown to support extremely large nonlinear coefficients and low incident power requirements to trigger nonlinear responses [45]. In addition, using high quality-factor Lorentzian [46] or Fano [47] resonant metasurfaces, it is possible to further boost the nonlinear response.
Based on these design requirements, we present two implementations of nonlinearitybased nonreciprocal optical metasurface devices. The first design is based on low quality-factor metasurfaces, which are therefore less prone to fabrication errors, but requiring large intensities in the range of GW/cm 2 . The second design is based on high qualityfactor resonators, which are more sensitive to fabrication errors but require much lower input power levels. The results for the proposed designs are plotted with red asterisks in Figure 1b. Some of them lie in the red shaded region, thereby overcoming the fundamental limitation Equation (5).

Low Quality Factor Resonators
As illustrated in the previous section, we designed two optical Fano resonators with shifted resonance frequencies such that they have the same unitary transmission at the same incident power level; in addition, one exhibits a bistable transmission curve, while the other smoothly varies as we increase the input power for the individual resonators. Next, we will arrange the resonators in order to realize large optical nonreciprocity.
In order to design a Fano resonator, we consider a thin dielectric slab of Si ( Si = 12) of thickness d 1 with an etched groove on the top surface of thickness d 2 patterned periodically with period p, as shown in the inset of Figure 4a. We choose the period p = λ/2 to couple the normally incident wave to the guide mode of the unpatroned slab. This structure supports Fano resonances based on two coupled modes. In order to inspect the two modes in this periodic array, let us consider first the special case d 2 = 0, so that the thick dielectric slab acts like a Fabry-Pérot resonator with a broadband resonance, i.e., a bright mode. When we next structure the slab with periodic grooves, a normally incident plane wave, with a magnetic field polarized along the z axis (TM), can couple to the dielectric waveguide mode, which is essentially a dark mode and cannot be accessible without the presence of the grooves. Upon the TM excitation, the two modes interfere, and the transmission spectrum experiences a Fano line shape as a function of the incident frequency f, as shown in Figure 4a for different parameters of d 1 and d 2 . In order to confirm the role of interference between the two modes, we fit the transmission spectrum with a standard CMT model, which indeed shows excellent agreement with the full wave simulation (COMSOL Multiphysics, Stockholm, Sweden). In order to further confirm that the waveguide mode is excited, we plot the magnetic field component H z , as shown in the inset in Figure 4a, which indeed shows a field distribution very similar to the TM 01 dielectric waveguide mode [42]. The quality factor of the resonance is in the order of 10, suggesting that the design is not very prone to fabrication errors and imperfections, and has a wide bandwidth.
We choose (d 1 and d 2 ) for each resonator, such that they satisfy the requirements illustrated in the previous section. In order to confirm this, we numerically investigate the nonlinear response of the individual resonators by including in the simulations the Kerr nonlinearity of Si, Si (r) = 12 + χ (3) |E(r)| 2 , where |E| is the magnitude of the electric field, and r is the position vector. Figure 4b shows the response at the normalized resonance frequency, f/(0.5 c/p) = 1 (vertical dashed line in Figure 4a), with a high-power incident plane wave with intensity P inc . The transmissions for two different structures with parameters shown in the legend as a function of the input intensity are plotted in Figure 4b, showing that they both exhibit very high transmission T = −0.03dB at the same power level, P inc = 15 GW/cm 2 , as shown in the inset at point ℵ, confirming the design requirements presented in the previous section.  Table  3 in Section 5.5 (Fano Resonator Parameters). The color plot inset shows the z-component of the magnetic field distribution inside the slab, consistent with the guided mode inside the unperturbed dielectric slab of the same thickness, and the color bar is normalized to the maximum field. (b) Transmission coefficient for individual resonators in dB versus the incident power for normal incidence assuming χ (3) = 2.8 × 10 −18 m 2 /V 2 at the resonance frequency.
We can now form a nonreciprocal device by arranging the two periodic resonators back-to-back, separated by an electrical length θ (a top grooved slab with d 1 = 0.18λ 0 and d 2 = d 1 /4, separated by a thickness l = λθ/(2π) and a bottom grooved slab with d 1 = 0.203λ 0 and d 2 = d 1 /3, p = λ 0 /2 where λ 0 = 1.55 µm, and θ = 1.5π), where one unit cell is shown in the inset of Figure 5a. We consider high power excitation (H z ) with increased input intensities (P inc ) from top and bottom, and plot the transmission coefficients T 12 and T 21 at different operating wavelengths. Figure 5a shows the transmission coefficient when excited from different sides at the operation wavelength λ 0 . As expected, at a low incident power P inc < 10.5 GW/cm 2 , the system is almost reciprocal, resulting in T 12 ≈ T 21 . As we increase the incident power, specifically at P inc = 10.5 GW/cm 2 , T 21 experiences a sharp transition from 0 to 0.95, while at the same time T 21 remains very low. As we further increase the input intensity, the transmission T 12 exhibits a similar abrupt transition from 0.07 to 0.98. The intensity range between the two transitions defines our NRIR = 0.27 dB, and the highest transmission between the two transitions is −0.132 dB, plotted in Figure 1b, showing that at wavelength λ 0 = 1.55 µm the device does not overcome the bounds of an optimal single resonance device.
We can also calculate the transmission at larger wavelengths λ 0 = 1.56 µm and at λ 0 = 1.565 µm, as shown in Figure 5b,c, respectively. Interestingly, the response at λ 0 = 1.565 µm shows an NRIR = 0.984 dB and the highest transmission T = −0.037 dB, and it indeed shows that the response breaks the limitation of a single Fano nonlinear device, as denoted by its asterisk point in Figure 1b. However, the design at λ 0 = 1.56 µm still satisfies the bounds for a single nonlinear resonator, as indicated in Figure 1b. We also calculated the transmission in the dB scale, finding that T 12 = −50 dB at the peak transmission for T 21 when λ = 1.565 µm, demonstrating a large contrast in transmission, a sort of pseudo-isolation of 50 dB for this device at an incident power P inc = 18 GW/cm 2 . We stress that we should not consider this quantity as a conventional isolation metric, because it is based on the assumption that the two ports are excited independently, and the device cannot fully isolate in the general case of continuous wave excitation from opposite ports. At this level of incident power and frequency, we show the full-wave simulation results of the field distribution when the device is excited from each side in the inset of Figure 5c, which indeed shows close to unitary transmission in the forward propagation, while the field is largely reflected in backward propagation, forming a standing wave at the opposite port. Notice that this functionality is very different from a conventional isolator, which would instead necessarily absorb the input energy coming from the input port. This reflectivity can be used to our advantage, for instance by routing it to a third port for circulation purposes [1]. Overall, the results in Figure 5 are far superior to those reported in [23,24] for designs using the same material parameters; however, here we use larger input intensities because of the low-quality factor of the involved resonators. As shown in Figure 5a-c for each frequency, the device works for slightly different power levels. Therefore, in order to characterize the response in the frequency domain at a fixed input power, we plot in Figure 5d the transmission versus the frequency at a fixed input power of 16.5 GW/cm 2 , which shows the isolation of the bandwidth greater than 600 GHz around the telecommunication wavelength, confirming the wideband response of the device due to the low quality factor of the constituent resonators, at the cost of increased required power levels.

High Quality Factor Resonators
As shown in the previous section, the power levels required for the device to work in its optimal condition are in the order of GW/cm 2 , due to the low-quality factor of the resonator, which may hinder its broad applicability. However, one possible solution is to employ materials with high nonlinearity, e.g., multiple quantum wells have shown extremely high nonlinear susceptibilities both in both the second and third orders [45,48,49].
Higher-Q Fano resonators can be used to reduce the required power, however, and for this reason we aim to maximize the field enhancement levels inside the resonators by designing a high-Q Fano resonator based on guided mode resonance coupling through narrow slits in a dielectric slab, as shown in inset of Figure 6a. The transmission coefficient for a normally incident plane wave, with an electric field in the z-direction (TE), is shown in Figure 6a, exhibiting a typical Fano line shape with zero transmission at λ 0 = 1.5608 µm and unitary transmission at the very nearby wavelength λ 0 = 1.5598 µm. In order to evaluate the quality factor of this resonator, we plot the average electric field enhancement inside the nonlinear material |E| 2 , as shown in Figure 6b, which shows a Lorentzian line shape with quality factor Q = 1700, compared to Q = 10 of the design in the previous section. Additionally, the inset shows that the electric field distribution inside the slab is similar to the TE 01 waveguide mode [42]. The second resonator is formed by adding a glass slab (n glass = 1.5) to another Si grating, as shown in the inset of Figure 6c. Similar to Figure 6a, the transmission exhibits a sharp Fano resonance; however, the average electric field enhancement is smaller than in Figure 5d. This is understood, as we chose the guided mode resonance to be formed inside the glass slab, and not in the Si layer, as confirmed by the electric field distribution in the inset of Figure 6d. We chose the second resonator to be attached to the glass substrate in order to control the field enhancement in the Si layer and hence obtain unitary transmission at the same power levels when the nonlinearity is included (similar to point ℵ in Figure 4b). In addition, we plot the CMT standard model results in Figure 6a, and Figure 6c using the fitting parameters from Table 3 in Section 5.5 (Fano Resonator Parameters), which shows excellent agreement with the full wave simulation.  Table 3 (see Section 5.5 (Fano Resonator Parameters)) for resonator 1 (w/o glass). (b) Average squared value of the electric field inside the Si material shown in (a); the inset shows the electric field distribution at resonance, and the scale bar is normalized to the maximum field. (c) The same as (a), but with an added glass slab. The red markers are the CMT standard model results using the parameters given in Table 3 (see Section 5.5 (Fano Resonator Parameters)) for resonator 2 (w glass). (d) The same as (b) but for the structure in (c). The geometry parameters used are d = 0.1179 µm, p = 1.4763 µm, dg = 0.45875 µm, and slit width = 0.0536 µm.
By exciting each metasurface separately with a normally incident TE plane wave with input intensity P inc , we can calculate the transmission as a function of the intensity, considering the Si nonlinearity using full wave simulation (COMSOL Multiphysics). The nonlinear transmission coefficient is shown in Figure 7a for the two resonators independently. We stress that the excitation for each structure from different sides results in a symmetric transmission curve following identically the same curves in Figure 7a, essentially having NRIR = 0 dB. Following a similar design process as above, we are able to achieve unitary transmission for each individual resonator at the same incident power level P inc = 1.5 MW/cm 2 (point ℵ in the inset of Figure 7a), which is much smaller than the input intensity of 15 GW/cm 2 obtained in the previous section. Additionally, we can plot the CMT results for this nonlinear behavior, which shows excellent agreement with the numerical simulation. We can form an optimal nonreciprocal device by stacking the two gratings of Figure 7a, separated by the distance l = λθ/(2π), where λ = 1.56104 µm and θ = 6.4π, as shown in the inset of Figure 7b. The nonlinear transmission coefficient for excitation from opposite sides for increased input intensities P inc is shown in Figure 7b. We obtain very high transmission T = −0.2 dB and large transmission contrast of over 60 dB at an incident power level of only 1.4 MW/cm 2 . Interestingly, the NRIR = 2.1 dB, and by plotting the point (NRIR, T) in Figure 1b, we see that this design exceeds the limitation of a single Fano nonlinear resonator. We emphasize that this power level is very small, and can be achieved even with continuous wave lasers; however, this comes at the cost of being vulnerable to fabrication errors [50]. Concepts borrowed from topological photonics may be used to enhance the robustness of the response in high-quality factor systems, even in the presence of fabrication errors [51,52].  Table 3 (see Section 5.5 (Fano Resonator Parameters)). (b) Transmission for the coupled metasurface design with θ = 6.4π shown in the inset. T 12 indicates the transmission from top to bottom, while T 12 indicates the transmission from bottom to top. The NRIR is 2.1 dB, and the peak forward transmission is −0.2 dB.
We also studied the response of this optimal device for different electrical lengths θ in Figure 8. Due to the nonlinearity of the system, the response is interestingly not necessarily periodic, even after the evanescent field effects disappear. At very large input powers, past the bistable response region, the response becomes periodic, as shown for the case θ = 10.4π, which is very similar to θ = 6.4π in Figure 7b. Additionally, as illustrated above, the value of θ largely controls the direction of isolation; for example, the values of θ (θ = 3.6, 4.4, 5.2, 10.4π) give forward transmission and backward isolation, and the converse is true for the other values of θ = 3.8π, 4.8π. Additionally, there are some values of θ where the response becomes symmetric and the field enhancement becomes symmetric for different directions of propagation (θ = 4π, 4.2π). We can characterize the operation of the proposed metasurface by comparing it to other free-space nonreciprocal metasurface devices operated under different principles. In Table 2, we present a comparison highlighting the main aspects characterizing the response of these different nonreciprocal metasurfaces, focusing on the main approaches used to break the reciprocity, i.e., breaking time invariance [11], passivity [17] and linearity [23], including this work.
As seen in Table 1, we notice that the bandwidth of all of the proposed approaches is very small, and it does not exceed a 5% fractional bandwidth in the best case [19]. This can be illustrated because all of the approaches employ resonant elements with a high quality factor in order to reduce the power requirements, and as a consequence decrease the bandwidth. We also observe that it is often difficult to characterize the efficiency of a metasurface at RF frequencies because there is a non-trivial number of parameters that need to be taken into account in order to estimate the total power efficiency. For example, in time-modulated metasurfaces, the efficiency should be calculated by normalizing the output power to the total input power, including the modulation power. In fact, it is very valuable that the authors in [13] did a full characterization of the metasurface performance by calculating all of the relevant power quantities. Although the insertion loss in the their work is about −5dB, this does not mean that the efficiency is low compared to other works based on the same concept, because other works [11,14] do not provide a full estimation of the relevant power quantities. However, we emphasize that for nonlinear metasurfaces, it is very straightforward to characterize the overall efficiency because the signal self-modulates the medium through the Kerr nonlinearity, and the metasurface itself is passive.  [16]. £ Based on the datasheet of the amplifier and the DC bias used in the paper and similar transistors used in previous publications [17], it is expected to be 0.1-0.2 W [18]. * From Figure 3c [23]. ** From Figure 3 [22]. *** From Figure 4 [22]. **** From Figure 5 [22]. ¿ Notice that the NRIR of this metasurface does not break the limitation of the single Fano resonator nonlinear isolator, even though 1.52 dB NRIR goes beyond the limit, the isolation is not infinite, so a fair comparison cannot be guaranteed here.

Conclusions
We studied and designed nonlinearity-based nonreciprocal metasurface devices based on coupled Fano resonances. Their response has been shown to overcome the fundamental limitations of a single resonant element associated with time-reversal symmetry [27,44], and our designs show an insertion loss of −0.04 dB for an NRIR of 1 dB with isolation over 50 dB. Our results are presented in normalized units, such that they can be scaled to other frequency ranges, and even applied to other materials, such as multiple quantum wells supporting larger nonlinearities [45]. We focused on two designs, with one supporting a moderate resonance and therefore requiring large power levels, but with robust response to fabrication errors. The second design is more prone to fabrication errors, but it requires much smaller power levels to operate. It is, of course, possible to explore designs that work in between these two extremes, as a function of the desired trade-off between power levels and the robustness of the design. It is interesting that all of the considered designs, while surpassing the bound for single resonant elements in Equation (5), lie very close to it, as seen in Figure 1b. It appears that metasurface designs do not provide the same flexibility of circuit elements, which provided a superior performance in terms of NRIR [34]. Our work provides useful guidelines to design these types of nonlinearity-based nonreciprocal metasurfaces for free-space radiation. Similar concepts can also be extended to photonic integrated circuits for which the control over the structure parameters is easier, and it provides a useful platform for passive, bias-free, magnetic-free CMOS-compatible nonreciprocal components for optical communications, sensing, imaging, and computing.

Coupled Mode Theory
In order to illustrate the formation of a Fano resonance, we consider in Figure 9 two generic resonant modes coupled to radiation channels. We assume a complex field inside resonator n, a n such that |a n | 2 is the normalized energy inside the resonator. We derive the CMT equations as Assuming excitation from the left port with the form S 1+ = e jωt , we rewrite the steady-state equations, assuming that all modes and signals have the form A → Ae jωt , as where, δ n = ω − ω n . Then, we solve for a 1 and a 2 , and then substitute in S 3− to obtain the transmission coefficient. (8) There are two zeros in transmission, when δ 1 = 0 or δ 2 = 0. Furhtermore, the unitary transmission happens at ω = κ 2 ω 1 +κ 1 ω 2 κ 1 +κ 2 , or equivalently when κ 1 δ 1 = − κ 2 δ 2 . Throughout the work, we assume that the second resonance frequency is much larger than the first resonance frequency ω 2 ω 1 , with a much higher coupling rate, κ 2 κ 1 .

Nonlinear Bistability
If we assume the resonators to be nonlinear, and ω 2 is far away from the frequency range of interest, the resonance frequency of the mode a 1 changes according to the relation [42], where ω 0 is the resonance frequency of the mode a 1 before including nonlinearity, and |a 0 | 2 quantifies the nonlinearity, specifically, the rate of the resonance frequency shift due to nonlinearity. We can substitute for |a 1 | 2 from (10) into (6) and obtain where, again, δ 01 = ω − ω 0 . This is a cubic equation of the only unknown δ 1 when all other parameters are known, and it has the form where the shorthand parameters (P, , †) are given by By solving (12), we obtain the shift in resonance frequency including nonlinearity, which can later be substituted in (8) to obtain the nonlinear transmission coefficient.

Bistability Condition
The bistability of a single nonlinear Lorentz resonator can be achieved at a certain range of power levels if the resonator is excited with a frequency ω, such that [42] where ω 0L is the resonance frequency of the resonator at which the maximum energy is stored when exited; κ 1L is its total decay rate, and we added the subscript L to refer to Lorentz resonator. We can derive a similar condition for our two-mode Fano resonator. This can be accomplished by writing the energy in a 1 of our Fano resonator in a similar form to the energy of a Lorentz resonator, and then finding the resonance frequency ω mod and coupling rate κ mod that are equivalent to ω 0L and κ 1L . In order to derive these effective parameters, we start from the fact that the Fano resonant modes satisfy the relation |ω 2 − ω 1 | δ 1 , which means that we are interested in a solution in the vicinity of the high-quality mode a 1 . Therefore, a good approximation is Then, we rewrite the mode energy |a 1 | 2 from Equation (6) as It is easily shown from Equation (14) that the energy inside the Fano resonator follows a Lorentzian lineshape, such that we can proceed to derive the equivalent quantities in order to derive the bistability condition. Additionally, we can see that the largest energy is achieved at a frequency that lies between the zero and unitary transmission of the Fano resonator. Thus, this frequency is equivalent to ω 0L , and it is given by Recall that ∆ ω < 0 so ω mod > ω 1 . We can also derive the effective coupling rate (linewidth) for this lineshape as Therefore, the condition for bistabitlity in the proposed Fano resonator is given by

Effective Power of Coupled Nonlinear Fano Resonators
When we connect two Fano resonators (a, b) with shifted resonance frequencies, as shown in Figure 10, we can calculate the incident power to the right Fano resonator when the device is excited from the left with incident power |S 1a+ | 2 . It is worth mentioning that each Fano resonator in Figure 10 can be modeled as a coupled mode resonator, as shown in Figure 9, and that each resonator can be characterized by its scattering matrix, for instance, resonator (a, b) to the (left/right) has the field reflection coefficient (r a , r b ) and transmission coefficient (t a , t b ). We calculate P eff = |S 1b+ | 2 under the condition that we operate the device in power levels that ensure complete reflection from the second resonator, i.e., S 1b− = e jφ b S 1b+ , where φ b is the reflection phase from the second resonator. In addition, we assume that the first resonator reflection coefficient in the form r a e jφ a , which is almost constant with increasing power levels. We can easily write the following equations: We use the relation |r a | = √ R a = √ 1 − T a , and define |S 1a+ | 2 = P, so finally:

Fano Resonator Parameters for Figures 3, 4, 6 and 7
We can interpret the Fano resonances from the structured Si layers used throughout this work based on their definition as the interference of a broadband bright resonance and a narrowband dark resonance. The dark resonance can be coupled to the background continuum when it is coupled evanescently to some bright mode, or through diffraction by perturbing the mode parameters, as in our case. In the case of guided mode resonance, the bright resonance is the Fabry-Pérot resonance supported longitudinally by the unperturbed dielectric slab, while the dark resonance is the guided mode of the dielectric slab, which is perturbed by the thin grooves in the dielectric. These two modes couple to the background with different coupling coefficients, and indeed they have different resonance frequencies.
In order to derive the result in Figure 4, we find the best fitting parameters as given in Table 2. It is interesting to see that the value of the coupling coefficient of the dark mode κ 1 is smaller than the coupling of the bright mode κ 2 . Table 3, on the other hand, gives the fitting parameters for Figures 3, 6 and 7, which also show that the coupling coefficient of the dark mode is much smaller than that of the bright mode.  The characterization of the metasurface was performed using linear and nonlinear simulations in the commercial software COMSOL Multiphysics. They are based on finite element numerical simulations, where the solution domain is divided into nonuniform free triangular elements, and Maxwell's equations are solved inside and on the boundaries of these elements enforcing the boundary conditions. Additionally, the software enables us to define the permittivity values as a function of the local field power. Each of the presented results contains two simulation steps: a linear simulation where the optical parameters of the materials are considered linear or equivalently, assuming very low power operation, and a nonlinear simulation. We first performed linear simulations to obtain the linear reflection and transmission coefficients, and to obtain the resonance of the structure. In this simulation, we assume a unit cell of the metasurface with periodic boundary conditions on the left and right port, and excited the structure with periodic boundaries from top/bottom, while keeping the bottom/top port as receiving ports, as shown in Figure 11. Additionally, we added a nonuniform rectangle mesh to the unit cell, assuming a minimum element size of 38 nm inside the dielectric. A mesh view for the low-quality factor structure is shown in Figure 11 (left panel), and a distribution map of the mesh cell area is shown in Figure 11 (right panel). Notice that a similar mesh is used for the high-quality factor structure; however, it is important in this kind of structure to keep the mesh fixed in both the linear and nonlinear simulation, because any slight change in the mesh would significantly change the resonance location. After obtaining the linear response of the metasurface encoded in the reflection, transmission and energy, and determining the frequency of interest as presented in the main text, we updated the silicon permittivity using a user-defined function, such that the permittivity is written as = linear + χ (3) |E| 2 , where linear is the permittivity used in the linear simulation, and the effect of nonlinearity is encoded in the second term of the permittivity. Notice that |E| 2 is not uniform inside the dielectric; therefore, our simulation takes into account all spatial inhomogeneity in the electric field. Then, we excited the unit cell with a plane wave of fixed frequency, and performed a parametric sweep over its power, in order to obtain the transmission from the top or bottom, respectively. Because the permittivity is a function of the induced field, the convergence of the solver is not always guaranteed. In order to overcome this problem, we increased the number of the iterations in the iterative solver to over 100 iterations in order to make sure that the solution will converge. In some situations, especially for the high Q structure, the convergence was not realized using the parameteric sweep, so we had to revert to manually updating the power and using the solution from the previous step as an initial condition for the next step. This is similar to following one curve on the bistability by either increasing or reducing the power. We repeat this step, but for excitation from the bottom to top in order to obtain the full transmission curves T 12 and T 21 .
Overall, apart from the convergence and manually updating the initial condition, the solution time for both the linear and nonlinear problem takes less than five minutes, assuming 100 points of parametric sweep of power and wavelength. The simulation was performed on a personal laptop.