Intersubband polaritonic metasurfaces for high-contrast ultra-fast power limiting and optical switching

Nonlinear intersubband polaritonic metasurfaces support one of the strongest known ultrafast nonlinear responses in the mid-infrared frequency range across all condensed matter systems. Beyond harmonic generation and frequency mixing, these nonlinearities can be leveraged for ultrafast optical switching and power limiting, based on tailored transitions from strong to weak polaritonic coupling. Here, we demonstrate synergistic optimization of materials and photonic nanostructures to achieve large reflection contrast in ultrafast polaritonic metasurface limiters. The devices are based on optimized semiconductor heterostructure materials that minimize the intersubband transition linewidth and reduce absorption in optically saturated nanoresonators, achieving a record-high reflection contrast of 54% experimentally. We also discuss opportunities to further boost the metrics of performance of this class of ultrafast limiters, showing that reflection contrast as high as 94% may be realistically achieved using all-dielectric intersubband polaritonic metasurfaces.


Introduction
Nonlinear optical materials with large and fast Kerr ( (3) ) nonlinearity are highly desirable for a wide range of applications, including frequency combs and short laser pulse generation, protection of sensitive equipment via optical power limiting, photonic information processing, and all-optical computing [1].To address these needs, a number of approaches have been pursued with the goal of achieving large values of  (3) in metasurfaces, thin-film and two-dimensions materials, such as using intrinsic metal or dielectric nonlinearities [2][3][4][5][6], nonlinearities of semiconductors and conduction-to-valence band transitions in semiconductor heterostructures [7,8], two-dimensional materials [9,10], epsilon-near-zero materials [11][12][13][14][15] and J-aggregates or other molecules [16][17][18].Large nonlinearities and subwavelength optical thicknesses enable operation with very tight light focusing in the focal plane, while avoiding the phase-matching constraints of bulk nonlinear materials.
Sub-picosecond Kerr nonlinearities have been recently reported in intersubband polaritonic metasurfaces operating in in the mid-infrared (3-30 m) spectral range [19].The operating principle of these structures (Figs.1a-b) is based on a fast change from strong to weak coupling between a cavity mode and an intersubband (ISB) transition in a multi-quantum-well (MQW) semiconductor heterostructure [20][21][22][23].These metasurfaces yieldto datethe largest picosecond third-order nonlinear response in condensed matter systems in the infrared spectral range, with effective  (3) values on the order of 3.410 -13 m 2 /V 2 , and a response time shorter than 2 ps for operating wavelengths around 7.5 m [19].
While large nonlinearities are beneficial in reducing the required intensity levels, for practical applications such as optical power limiting, saturable absorber mirrors for laser mode-locking and frequency comb generation in the mid-infrared spectral range, and the development of nonlinear optical elements for optical computing and all-optical control, achieving large absorption/reflection contrast is equally important.Here and in the following, we define the reflection contrast of a metasurface limiter as the difference between the minimum and maximum reflection levels within a certain input power range.The metasurfaces reported in Ref. [19] and similar structures reported recently in [24] displayed moderately low power-dependent reflection contrast, with a maximum reflectivity change of about 25-30% reported in [19].
In this work, we present Kerr nonlinear metasurfaces based on intersubband polaritons which achieve a nonlinear absolute reflection contrast of 54%, while simultaneously displaying recordhigh values of effective  (3) .Such enhancement is due to the combination of an optimized MQW system based on Ga0.53In0.47As/GaAs0.51Sb0.49heterostructures, supporting ISB transitions with extremely narrow linewidths, more than a factor of 4 narrower than those of the Al0.48In0.52As/Ga0.53In0.47Asheterostructures used in Ref. [19], optimized nanophotonic structures leveraging these material resonances, and vertical, rather than diagonal, intersubband transitions [19], which maximize the transition dipole moment.The measured reflection contrast of 54% is, in fact, limited by the range of input powers available in our setup, and we estimate that reflection contrasts above 60% can be observed within the same realized device.Besides largely increasing the reflection contrast, the narrow linewidth of our ISB transitions also unveils exotic absorption features associated to the bending of the semiconductor bands close to the Schottky metal contacts, which may be leveraged to further improve metasurface performance.The main factor limiting the reflection contrast performance in our devices is the residual absorption in the metallic nanoresonators.We theoretically show that, by combining our optimized MQW system with all-dielectric metasurfaces, reflection contrasts as high as 94% may be realistically achieved.

Results
Physical model.We begin by briefly outlining the physics behind the nonlinear response of a polaritonic metasurface induced by the saturation of ISB transitions in MQWs [19].We consider a metasurface with a unit cell schematically shown in Fig. 1a, formed by a stack of MQW layers interacting with a metal patch antennaa platform that has been investigated in several recent works [19,23,24].This system is analyzed by a combined coupled-mode theory/Maxwell-Bloch (CMT/MB) model [19], modified to account for the absorption from higher-energy electron subbands, and described by the equations = 4Im( * ) +  1 ( + 1), The patch antenna is modeled as a low-loss optical resonator with resonant frequency   = 2  , complex modal amplitude a (normalized so that || 2 is the number of stored photons), and radiative and nonradiative decay rates   and   , respectively.The optical resonance of the antenna can be directly excited by an incident plane wave with time-dependent amplitude  + (normalized so that | + | 2 is the incident power per unit cell area ) and carrier frequency  .The backpropagating outgoing field  − (Eq.4) is given by the superposition of the direct reflection of the incident field and the field leaked by the optical resonance.The MQW is modeled by N identical two-level quantum oscillators, each with transition frequency   = 2  , total decoherence rate  2 and relaxation rate  1 , coupled to the optical resonance with rate  = , where r is the relative permittivity of the MQW system at zero doping,  12 is the transition dipole moment from state 1 to state 2, and V is the mode volume.This effective two-level system describes the ISB electronic transition between the ground state and the first excited state of the conduction band (inset of Fig. 1a), assuming that N electrons are in the ground state at thermal equilibrium.The state of each two-level emitter is described by the off-diagonal element of the reduced density matrix, denoted by  (Eq.2), which describes the polarization of the transition, and by its inversion  =   −   , where   and   are the electron populations in the ground and first excited electron levels in the MQW system, normalized to the total electron population.Additionally, Eq.
(1) contains an effective loss term   ( + 1) [19], which accounts for light absorption due to transitions of the excited electrons in higher energy states.This term becomes relevant at high pump intensities, depending on how significant this additional absorption is compared to the nonradiative decay rate   of the cavity and the Rabi splitting (see Eq. ( 5)) at low intensity.The inversion  is the source of nonlinearity in our system.At low intensities (|| ≈ 0), the inversion approaches  → −1, and Eqs. ( 1)-( 2) describe two classical coupled oscillators.In this linear regime, when the system is degenerate (  0 =   =   ) and in strong coupling, the reflection spectrum () ≡ | − ()| 2 /| + ()| 2 (blue line in Fig. 1b) features a characteristic Rabi splitting, displaying two distinct polaritonic resonant dips with equal linewidths and depth, and with resonant frequencies symmetrically detuned from the central frequency  0 by ±Ω , where Ω ≡ √.As a result of this splitting, the reflection is large at  =  0 , approaching the limit ( 0 ) → 1 for  ≫ { 1 ,  2 ,   ,   ,   }.As the impinging intensity increases, the saturation of the ISB transition causes the inversion to approach zero ( → 0), and the ISB transitions (Eq.2) are completely decoupled from the optical resonance.In this regime, the reflection spectrum of the system approaches the one of the bare patch antenna (red line in Fig. 1b), consisting of a single reflection dip centered at the frequency  =  0 =   .Because of this saturation effect, the energy stored in the optical resonance, and thus the outgoing field (Eq.4), depends nonlinearly on the incident intensity | + | 2 .In particular, at a fixed excitation frequency  =  0 (yellow arrow in Fig. 1b) the reflection strongly decreases as the input intensity (and thus the emitter inversion) increases.Instead, for excitation frequencies close to the frequencies of the polaritonic dips ( 0 ± Ω, green arrow in Fig. 1b), the reflection level evolves from low (at low intensities) to high values (at high intensities).The evolution of the reflection versus saturation, defined as  + 1 , is displayed in Fig. 1c for excitation frequencies equal to  0 (yellow line) and  0 + Ω (green line).
The relaxation time 1/ 1 is typically the slowest time constant in Eqs.(1)-(3).Thus, it dictates the timescales over which the inversion and saturation evolve in time, and it determines the speed of the metasurface nonlinear response.Values of 1/ 1  1-2 ps have been determined experimentally [19], confirming that this system can be used as a fast power limiter or a fast saturable mirror.Another important figure of merit towards the practical application of these devices for power limiting is the reflection contrast Δ, defined as the difference between the reflection of the system in the fully linear case ( = −1) and the fully saturated case ( = 0), i.e., Δ ≡ ( = −1) − R( = 0) (see also vertical black arrow in Fig. 1b).For an ideal power limiter Δ = 1, while an ideal saturable absorber features Δ = −1.

Reflection contrast limits.
In [19] we showed that the achievable reflection contrast Δ is fundamentally limited by the different decay channels of the system  2 ,   ,   , and   .In particular, if low reflection at high powers is desired, i.e., for a power limiter operation, the system needs to be close to critical coupling when fully saturated, i.e.,   ≈   +   .In this limit, we find that where the upper limit of inequality is achieved when   = 0. Equation (5) shows that  is limited by the largest decay channel of the system.Since the ratio  2 / is present in all the 'damping' terms in Eq. ( 5), large values of || can be achieved with ISB systems featuring small decoherence rates  2 and large vacuum Rabi frequencies , while keeping   and   as small as possible [19].As an example, Fig. 1d shows the calculated  as a function of  2 / and for a fixed value of   = 0 and   /Ω = 0.1.
In our earlier work [19], values of   / ≈ 0.04 ,  2 / ≈ 0.87 and   / ≈ 0.37 were demonstrated.Based on Eq. ( 5), for these parameters the maximum theoretical value of the reflection contrast is || ≈ 0.23 or 23%, which matches well the response reported in [19].As follows from Eq. ( 5), we can improve the reflection contrast by increasing the value of , reducing the product of   and  2 , and minimizing   .We analyze the practical limits of each of these approaches in the following.For a given operating frequency , the vacuum Rabi frequency can be written as [25] where e is the electron charge,   * is the effective electron mass,  12 = 2  *    12 2 /ℏ is the oscillator strength of the intersubband transition coupled to the cavity mode, N1 and N2 are the 2D electron densities in the ground and excited states in a single MQW period, and ℎ  is the thickness of a single MQW period.The maximum oscillator strength f12 is achieved for a vertical transition in a quantum well, with the value f12  0.96 for an infinite quantum well [26].Assuming a highly degenerate 2D electron gas, the maximum value of concentration difference is (N1-N2)= , which is achieved when the Fermi energy is   = ℏ  .For larger doping levels, both ground and excited states are being filled with electrons at the same rate.The minimum MQW period L is set by the requirement that a single quantum well supports a transition at the desired energy ℏ  , which determines the minimum quantum well width to be ℎ  = √ 3 2 ℏ 2  *   .In practice, the value of L is always larger than ℎ  because a quantum barrier is also required to separate individual quantum wells (in the limit of high conduction band offset, the barrier width can be very small, but in realistic systems it is often comparable to LQW).By including these equalities into Eq.( 6), we find that the maximum value of vacuum Rabi frequency is This value shows a weak dependence on the effective electron mass, and thus the choice of effective mass in the MQW system is not of particular importance [27].The relative dielectric constant can be assumed to be r 10 for most suitable ISB systems [27].Assuming  12 = 0.96,   * ≈ 0.04  , where   is the mass of a free electron, and ℏ  ≈ 150 meV we obtain Ω max ≈ 0.26   .
In the system reported in Ref. [19] we obtained  = 0.13 ω q , which is lower than Ω max by a factor of ~2.This can be attributed to a combination of several factors: first, the doping was ~2 times lower than the one required to achieve the optimal value of N1-N2; second, the use of a diagonal transition in the structure in Ref. [18] led to smaller than optimal transition oscillator strength.;finally, the barriers in the MQW structure had roughly the same thickness as the quantum wells (i.e.,  ≈ 2  ).While the first two issues could be addressed and improved, this would lead to only to a modest increase in .Moreover, significantly reducing the barrier thickness is not a viable route, as thinner barriers result in an energy level anticrossing between adjacent quantum wells, which increases the intersubband transition linewidth.We further note that all damping constants in Eq. ( 5) scale, in first approximation, with the operating frequency, and thus changing the transition frequency   is not expected to improve the metasurface performance.
Based on this discussion, the Rabi frequency Ω cannot be increased significantly without compromising other performance metrics.Hence, we turn our attention to the possibility of increasing the reflection contrast Δ by instead reducing the damping constants  2 ,   and   in Eq. ( 5).The value of   is related to the optical loss of the bare optical cavity, and it is due to Ohmic loss in the metal and scattering on the nanoresonator walls.In the improved metasurface design reported in the next section, we reduce the roughness of the nanoresonator sidewall via optimization of the fabrication process, while keeping the same nanoresonator design as in Ref. [19].We note that the value of   can be also reduced by switching to dielectric nanoresonators [28,29], albeit at the expense of a somewhat reduced electromagnetic field confinement.
The value of   can be reduced by limiting the impact of the ISB transitions from the excited state to higher energy states.This approach is ultimately limited by the conduction band offsets available in mature MQW systems.The value of   for the maximum Rabi frequency can be estimated by assuming that all electrons are excited to a continuum of conduction band states, where they can be described as a free electron gas characterized by the Drude permittivity We can use perturbation theory to find an approximate expression for the loss rate [30], leading to   ≈   Im((  ))/(2 r ) .Assuming the optimal bulk doping concentration of electrons  =   *   ℏℎ  = 1.8 × 10 18 cm −3 derived for Eq. ( 7), a typical electron scattering rate for heavily doped semiconductor materials  ≈ 0.1 ps [31], and operating frequency ℏ = ℏ  ≈ 150 meV, we obtain an upper limit of   = 0.01   .In our earlier work [19], we estimated   ≈ 0.05   based on the nonlinear behavior of the polaritonic metasurface.However, we note that in Ref. [19] the quantum barriers in the MQW structure had a large conduction band offset (520 meV) and most of the electrons were not excited into the continuum for the range of intensities (up to 700 kW/cm 2 ) employed in the experiment.As discussed in [19], the value of   in this case instead likely captured intersubband absorption between higher bound states in the MQW system.In contrast, in the MQW system described in the next section the conduction band offset for the barriers is only 360 meV, facilitating electron excitation to the continuum via two-photon absorption processes, as in Fig. 1(a).Given a factor of 10 larger excitation intensities (up to 10 MW/cm 2 ) employed in the experiments reported here, Eq. ( 8) provides a good estimate for the upper limit of   for ideal doping levels.In practice, however, we operate in an intermediate intensity regime where the intersubband transition is nearly saturated, but the electrons are not fully ionized.
Finally, the ISB transition dephasing rate  2 is present in all the factors in Eq. ( 5).Thus, reducing the value of  2 can lead to significant improvements in the metasurface performance.The value of  2 is typically more than a factor of 10 larger than  1 , and it is due to roughness-induced scattering at the quantum well/quantum barrier interfaces and electron-electron and phonon scattering [32].For higher-doped materials, non-parabolicity [33] also contributes to increased  2 .
Record-contrast  (3) intersubband polaritonic metasurface.It was recently reported that Ga0.53In0.47As/GaAs0.51Sb0.49material systems [34] can produce ISB transitions with extremely narrow transition linewidths.The mechanism behind this material response has not yet been fully explained, but reduced interface roughness scattering and relatively smaller non-parabolicity likely plays a major role.Guided by the analysis in the previous section, we optimized the parameters of molecular beam epitaxy, such as growth speed and temperature, to significantly reduce the interface roughness scattering.As a result, we were able to produce Ga0.53In0.47As/GaAs0.51Sb0.49heterostructures with similar doping levels as in the Ga0.53In0.47As/Al0.48In0.52Assystem used in [19], but with a more than 4-fold reduction of the ISB transition linewidth.
Following the analysis presented above, this new MQW system was designed around a vertical ISB transition in a single quantum well, as shown in the inset of Fig. 1(a).Similar to the MQW systems used in previous metasurfaces, modulation doping was used to avoid scattering of electrons by impurities.The MQW stack was constructed by 15 repetitions, each consisting of an 8.2 nm wide Ga0.53In0.47Aswell separated by a 15 nm wide Ga0.53In0.47Asbarrier.Charge carriers within each well were introduced by two delta-like doping regions, with a concentration of  = 6.5 ⋅ 10 11 cm −2 , positioned 2.3 nm from each side of the well.To minimize band bending effects at the metal contacts, we have incorporated Ga0.53In0.47Asbuffer layers with an n-doping of 1 ⋅ 10 19 cm −3 on both sides of the MQW stack.The bottom buffer layer has a 5 nm thickness, while the top layer also acts as an etch stop layer during fabrication and measures approximately 15 nm in thickness.This approach closely aligns with the design detailed in the previous work from Ref. [19].Additional details of the sample wafer structure are given in Methods.The extremely narrow ISB absorption linewidth is confirmed via 45°-angled absorption measurements [33] shown in Fig. 1(e).The absorption measurement yields an ISB transition linewidth of 2 2 = 7.4 meV, which is over a factor of 4 smaller than that of the MQW structures used in [19].The measured spectrum yields an ISB transition energy of ℏ  ≈ 155 meV, close to the simulated value.The height of the absorption peak can be used to deduce the actual doping level in our structure, which was found to be   = 5.4 × 10 17 cm -3 .The corresponding sheet electron concentration density  2 = 1.3 × 10 12 cm −2 is 2.2 times smaller than the optimal value of 2.5 × 10 12 cm −2 derived in the previous section.The smaller-than optimal value of electron concentration originates from the nonlinear activation rates of silicon dopants experimentally observed in our material system.In particular, we observe that, in order to obtain the desired further two-fold increase of electron concentration under optimized growth conditions, it would be necessary to increase the concentration of silicon dopants by a factor of approximately 9.This much larger concentration of donor impurities increases the electron-impurity scattering rate, leading to a significant increase of the ISB transition linewidth,  2 .
To demonstrate the improved performance unlocked by this MQW system, we used this material to fabricate and test polaritonic metasurfaces.Figures 2(a-b) show SEM pictures of a fabricated sample.The unit cell of the metasurface consists of a single patch antenna, obtained by sandwiching a stack of InGaAs/GaAsSb MQWs of thickness H≈350 nm between a flat layer of metal and a 50-nm-thick rectangular gold antenna with in-plane dimensions Lx and Ly (Fig. 2b).The patch antennas are arranged in periodic rectangular arrays with total lateral dimensions of approximately 200 µm (details on the fabrications are provided in the Supplementary Material).Because of the strong field confinement and the subwavelength unit cells, the optical response of the metasurfaces is governed by the resonance supported by each unit cell, which in turn is largely controlled by the antenna dimensions.In order to systematically tune the interaction between the ISB transition of the MWQ and the antennas, we fabricated devices with different antenna lengths in the range Ly = 1 μm − 2.3 μm, while keeping lattice pitch (  = 2 μm and   = 4 μm) and the short side of the antenna (Lx ≈ 330 nm) fixed.
The color-coded plot in Fig. 2c shows y-polarized reflection spectra of devices with different values of Ly, measured with a commercial Fourier-transform infrared spectroscopy (FTIR) system (see Methods).The measurements show the typical signature of strong coupling: two distinct reflection dips, corresponding to the MQW and antenna resonances, which anti-cross as the antenna resonance frequency (dashed white line) approaches the ISB frequency (horizontal dashed-dotted blue line).When the two bare resonance frequencies are approximately equal (  =   ≈ 38 THz , corresponding to the vertical dashed green line in Fig. 2c,   ≈ 1845 nm), the measured reflection spectrum (shown in Fig. 2d as solid blue line) features two reflection dips located at frequencies larger and smaller than the common frequency  0 =   =   , and a reflection maximum at a frequency close to 0  .While the presence of these two polaritonic dips is expected in the strong coupling regime, the measured spectrum (solid blue lines in Fig. 2d) shows clear deviations from the conventional two resonator model (Eqs.1-4 and Fig. 1b).First, the two polaritonic dips are not symmetrically detuned from the common resonant frequency 0  .Second, the two dips have different linewidths and depths, with the high-frequency dip being sharper and deeper than the low-frequency one.Such asymmetry is not possible based on the strong coupling of two oscillators (see also Fig. 1b), independently of their mutual coupling strength, decay rates, or other parameters, and it strongly contrasts with previous experimental results on similar devices [19].The black dashed line in Fig. 2d shows the best fit of the experimental data by using Eqs.1-4, emphasizing the incompatibility between the two-resonator model and the measured data.
These deviations from our theoretical model are due to the occurrence of an additional ISB transition localized at the top metal-semiconductor interface, where the Schottky effect, in combination with the n-doped top buffer layer, strongly affects the electronic band structure.We note that the effect of this additional transition becomes apparent only in the measurement of the patch antennas, and it does not appear in the absorption measurements of the bare MQW materials (Fig. 1e), since in that case no metal interface is located next to the buffer layer.
To corroborate this statement, we numerically calculated the band structure of the full MQW stack by using a self-consistent 8-band  •  model (Fig. 2f).In our model, we assumed that the top titanium/gold layer has a Schottky potential of 0.23 meV [35] and we included the effect of the band bending due to charge carrier separation at the Schottky interface.The band structure calculations display a clear bending of the conduction band potential energy (solid black line in Fig. 2f) near the metal interfaces and within the top buffer layer region (see the sample structure description in the Supplementary Material).As a consequence, the ISB transition near the interface has a transition energy (130 meV, blue arrow in Fig. 2f) lower than the one of the electronic states existing closer to the center of the MQW stack (153 meV, green arrow in Fig. 2f).The calculated charge carrier density for the barrier ISB transition is   = 8 × 10 17  −3 , which is 1.5 times larger than the charge carrier density inside the MQW stack (  = 5.4 × 10 17  −3 ).We note that another ISB transition is expected to exist also at the bottom buffer layer.However, the bottom layer has a different geometry and smaller thickness (due to its different functionality), leading to a much larger ISB transition energy, exceeding 300 meV according to our calculations.Due to the much larger detuning, this additional ISB transition does not play a relevant role in our model.
In the presence of an ISB transition at approximately 130 meV in the top buffer layer, our metasurface spectra cannot be described by a standard two resonator model for the MQW-cavity interaction [19] [Eqs.(1)-( 3)].We emphasize that, while the additional ISB transitions are expected to exist in any polaritonic metasurface with MQWs close to metal layers [19,24], their effect is often hard to observe because of the large ISB transition linewidths.The effect instead is very evident here due to the very small linewidth of ISB transitions in our MQW material, which allows us to spectrally resolve these additional states.In order to validate and incorporate this effect into our model of the polaritonic metasurface, we expand our CMT/MB model by describing the dynamics of the MQW stack with two separate resonators, one describing the ISB transition at the center of the MQW (with electron density   ) and the other one describing the ISB transition at the top buffer layer [see Fig. 2(f)] near the metal contact (with electron density   ).The system is described by an expanded set of coupled equations: where now   and   (i = c,b) are the polarizations and the inversions of the ISB transitions in the center (c) and in the buffer (b) layer, respectively.Each MQW transition has a different resonant frequency   and total decoherence rate  2 = 1/ 2 , while the relaxation rate 1 1 1/T  = is assumed to be the same for simplicity.Moreover, each ISB transition couples to the cavity mode with a different Rabi frequency Ω  = √    .
We begin by considering the low-power regime (|| 2 ≈ 0) where the inversions are   = −1, and we use this corrected model [Eqs.( 9)-( 13)] to fit the measured reflection spectra (Fig. 2c) of devices with different antenna lengths.In the fitting procedure we assume that   ,  2 ,  1 and Ω  are fixed across all devices (since all devices were made of the same MQW material) while the antenna parameters (  ,   , and   ) were allowed to vary.For the MQW material parameters, we obtained the values   /2 = 36.53THz ,   /2 = 33.37THz,  2 /2 = 0.88 THz,  2 /2 = 3.15 THz, Ω  /2 = 2.83 THz, Ω  /2 = 2.66 THz.We note that the value of  2 is in excellent agreement with the value independently extracted from the fit of the bare MQW absorption spectra in Fig. 1e,  2 /2 = 0.89 THz.For the device corresponding to the vertical dashed green line in Fig. 2c, the retrieved antenna parameters are   /2 = 36.82THz,   /2 = 0.87 THz,   /2 = 1.46 THz.The value of  1 was computed to be 0.9 THz from the opticalphonon-limited electron lifetime in the first excited electron subband in In0.53Ga0.47Asquantum wells, similar to the one computed and used for the MQW system in Ref. [19].With the parameters obtained from the fit, and assuming that the antenna resonance wavelength scales linearly with its length, we can accurately reproduce the linear reflection measurements in Fig. 2(c), as shown by the color plot in Fig. 2e.The accuracy of our model is further revealed in Fig. 2d, where we plot the calculated reflection spectra when the antenna resonance frequencies and the   transition frequencies are equal (solid red line), alongside with the corresponding measured data (blue line).Despite some residual discrepancies, the spectra calculated with the model in Eqs. ( 9)-(13) (red line) correctly reproduces the asymmetric shapes of the two polaritonic peaks observed in the experiment (blue line).Moreover, the value of the transition energy ℏ  = 138 meV obtained from the fit agrees reasonably well with the one expected from the band structure calculations near a metal contact (130 meV, Fig 2f).
After having experimentally characterized and analytically interpreted the linear spectra of our polaritonic metasurfaces, we now show how their sharp spectral response, enabled by the narrow ISB linewidth, can be leveraged to achieve MQW-based saturable absorbers and power limiters with large reflection contrast.We focus on the device with spectrum shown in Fig. 2d (i.e.antennas with length   ≈ 1845 nm), corresponding to the case where the ISB transition at the center of the MQW (described by   in Eqs.9-13) and the patch antenna are resonant.The corresponding reflection spectrum is reproduced in Fig. 3a.We experimentally measured the reflection of this metasurface versus impinging intensity for different frequencies across the metasurface spectrum.The measurements were performed with a custom-built setup described in the Supplementary Material.The excitation is provided by a pulsed laser with tunable carrier frequency, a repetition rate of 80 MHz and a pulse duration of 2 ps. Figure 3b shows the measured nonlinear reflection for 8 different frequencies, corresponding to the color-coded vertical lines in Fig. 3a.As expected from our theory [Figs.1(b-c)], when the excitation frequency is close to the center of the gap between the two dips ( = 38.36THz, dark green curve), the reflection strongly decreases with increasing intensity, starting from a large value of about 0.64 for intensities smaller than 100 KW/cm 2 and reaching values smaller than 0.1 for intensities of about 10 MW/cm 2 , corresponding to a reflection contrast |Δ| ≈ 0.54 or 54%.
For frequencies closer to the center of the dips the reflection evolves from low to high values as the intensity increases.In particular, when the excitation frequency is resonant with the narrowest dip ( = 40.67THz, dark orange line in Figs.3a and 3b), the reflection evolves from 0.1 to almost 0.5.Different nonlinear behaviors are obtained for other frequencies, including scenarios where the reflection level features a local minimum for intermediate intensity values ( = 39.63 THz, light green lines).We also note that, for most of the curves in Figs.3b, the measured reflection contrast Δ is limited by the intensity range available in our experiment.In particular, at the lowest intensity considered in our experiment (~ 50 KW/cm 2 ) the saturation level of the ISB transition appears to be already non-negligible, as manifested by the fact that several curves in Fig. 3b have a non-zero slope already at low intensities.Thus, as discussed later, the reflection contrast achievable with this metasurface is expected to be even larger than what observed in Fig. 3b.The color-coded plot in Fig. 3d shows the measured nonlinear reflection for a larger and more refined set of excitation frequencies.As clear from this plot, the reflection spectrum evolves from an asymmetric two-dip lineshape at low intensities to a single-dip lineshape at high intensities.The high-intensity spectrum corresponds to the case of a fully saturated ISB transition, and thus the single dip is attributed to the response of the bare antenna.The low reflection minimum achieved by the high-intensity spectra (min [ high ()] ≈ 0.1) confirms that the antenna is close to the critical coupling condition   ≈   .From the measurements in Fig. 3d we can extract the frequency-dependent reflection contrast Δ() ≡  low () −  high (), which is plotted in Fig. 3c (blue symbols).The experimental reflection contrast varies between Δ ≈ −0.4 for the saturable absorber operation and Δ ≈ 0.54 for the power limiter operation, which is approximately a factor of 2 larger than the previous state-of-the-art reported in [19].
To model the nonlinear measurements in Fig. 3d, we numerically solved Eqs. ( 9)-( 13) by assuming that the input  + () is a Gaussian pulse with duration of 2 ps, consistent with the experimental setting, and we calculated the resulting output field  − () for different input peak intensities.Importantly, in these nonlinear calculations all system parameters are fixed by the fit done on the linear measurements, i.e., Figs.2c and 2e, except for the electron density of each ISB transitions,   and   , which determine the corresponding saturation intensities.These values, together with the estimated volume of the unaffected and bent quantum wells [see Fig. 2(f)], yield   = 5.4 × 10 17 cm −3 and   = 8 × 10 17 cm −3 .The numerical calculations (Fig. 3e) show excellent agreement with the experimental measurements, reproducing both the asymmetric lineshape at low intensities and the intensity-dependent evolution of the spectrum.The excellent agreement between simulations and measurements is highlighted in Fig. 3f, where we show the measured and calculated intensity-dependent reflection for two representative frequencies.Moreover, the frequency-dependent reflection contrast Δ() extracted from the calculations (orange solid line in Fig. 3c) agrees well with the one extracted from measurements (blue circles in Fig. 3c), and it reproduces both the absolute values of Δ() and the asymmetric tails.Overall, the excellent match between the experimental measurements and the theoretical predictions for both linear (Figs.2c and 2e) and nonlinear (Figs.3(c-f)) scenarios confirms that the dynamics of these polaritonic metasurfaces can be quantitatively described by expanding the simple two-level model [19] to account for the presence of additional ISB transitions, induced by the Schottky barrier formed near the interface between the MQW and the metal layers.As mentioned above, the measured reflection contrast is limited by the range of input power levels used in the experiment, and it could be increased by starting from lower input powers in the measurements in Fig. 3b.Indeed, the calculation shown in Fig. 3c (orange line) confirms that the reflection contrast of our device could be as high as ΔR≈0.64 if a larger range of powers is considered.
In the context of saturable materials and power limiting applications, an important figure of merit is the change of reflection (or absorption) upon small variations of the impinging intensity I, i.e., the derivative / .In Fig. 4a we show the experimental value of |/| versus impinging frequency (blue symbols), extracted from the data in Fig. 3d at an impinging intensity of about 200 KW/cm 2 .A slope |/| ≈ 0.8 ⋅ 10 −6 cm 2 /W is obtained for frequencies close to the center of the two-dip linear reflection spectra.The measured |/| agrees well with values extracted from numerical simulations (dashed blue line in Fig. 4a).Since our system does not transmit any power due to the metal backing, the absorption is  = 1 − , and thus / = −/.The slope / can be used to estimate an effective nonlinear third-order imaginary susceptibility Im[ eff (3) ], via [19] Im[ eff (3) ] = 2 0  2  2 3  eff (11) where is the effective nonlinear absorption coefficient, and d is the thickness of the metasurface.By assuming n = 1 and d = 400 nm , we obtain Im[ eff (3) ] = 2.2 ⋅ 10 −13 m 2 /V 2 for the central frequency  0 ≈ 38 THz.As mentioned above, the saturation of the ISB transitions is nonnegligible already at the lowest available intensities in our experiment.This suggests that even larger values of |/| may be observed at lower intensities.In order to numerically verify this effect, we repeated the calculations shown in Fig. 3e for lower impinging intensities (while keeping all the other parameters the same as in Fig. 3e), and we report the corresponding values of |/| versus intensity and frequency in Fig. 4b.This calculation shows that derivatives as large as |/| ≈ 4.5 ⋅ 10 −6 cm 2 /W are expected in our metasurface for peak intensities lower than 10 KW/cm 2 .In turn, this value corresponds to an effective nonlinear susceptibility Im[ eff (3) ] > 10 −12 m 2 /V 2 , i.e., at least a factor of 3 times larger than the record nonlinearity reported in [19].
Further improvementsdielectric intersubband polaritonic metasurfaces.As discussed in Section 3, a key limiting factor in the polaritonic metasurfaces described above is the absorption rate   , mostly due to the presence of metals.In principle, these absorption losses can be largely reduced by using all-dielectric metasurface designs [28,29].Using the expressions derived for   and , assuming   ≫ 1/, and setting the absorption loss   to zero, we can rewrite Eq. 5 as Interestingly, this expression implies that the performance of a dielectric polaritonic metasurface in absence of scattering or absorption losses does not depend on the actual doping density, but only on the quality factor of the intersubband transition  2 /  , the relative scattering time   in the ionized electron gas, and the oscillator strength.For the present MQW, we have  12 ≈ 0.8 and  2 / 0 ≈ 0.024, and we estimate   = 22.The reflection contrast for these parameters yields || ≈ 0.94 or 94%.With the maximum oscillator strength of  12 ≈ 0.96 the contrast increases to 0.95.
We can distinguish between local and nonlocal dielectric polaritonic metasurfaces.In the local approach, each element of the metasurface supports a Mie-like resonance.For example, strong coupling between an intersubband transition and a geometric resonance has been observed cylinders supporting a magnetic dipole [28].Crucially, the magnetic dipole has a field that is largely z-oriented, enabling efficient interaction between the transition and the optical mode.Nonetheless, it should be noted that in Eq. ( 12) we assumed that all electric fields of the resonance are overlapping and aligned with the intersubband transition.In the dielectric metasurfaces, a fraction of the field will be stored outside of the dielectric quantum wells or not polarized in the vertical direction, and as a result  will be reduced, together with the resulting contrast.
Due to the small mode volume, local metasurfaces generally have a relatively low quality factor.Hence, to reach critical coupling and full absorption at saturation, significant absorption losses are also required, which in the present design are provided by the metal.In a full-dielectric design, absorption would be reduced significantly once the intersubband transitions are saturated, which means that achieving critical coupling to the small remaining losses is challenging with a local metasurface.This means that such local metasurfaces work best as saturable absorbers (e.g., as the current metasurface does at 41.5 THz in Fig. 3c), going from strong absorption at low intensity at one of the polaritonic branches to full reflection at saturation.One could also achieve limiting in local metasurfaces by using alternative "loss" mechanisms such as diffracted orders or transmission, but a more direct approach may be to use nonlocal dielectric metasurfaces instead.In these metasurfaces, the mode is delocalized and extends over many unit cells of a (quasi-)periodic structure.This enables much larger quality factors: in fact, such nonlocal metasurfaces may exhibit quasi-bound states in the continuum (Q-BIC) with arbitrarily high quality factors [36].This control over the quality factor enables critical coupling to low absorption loss rates, which for example has already been utilized to achieve unity emissivity and absorptivity for thermal emitters [37].Starting the design from a symmetry protected BIC with out-of-plane electric field ensures efficient coupling between the optical mode and the intersubband transition.Hence, depending on the desired application, using either local or nonlocal dielectric metasurfaces is a promising path towards even larger absorption contrasts than those demonstrated in the present paper.

Discussion
In this work, we have theoretically and experimentally investigated optimized nonlinear metasurfaces based on ISB polaritons in MQW systems.We have provided a comprehensive theoretical analysis showing how the key metrics of these devices for limiting operations, such as the reflection contrast Δ, can be linked to the intrinsic material properties of the underlying MQW system and the metasurface parameters.Our analysis elucidates that the achievable values of Δ can be largely increased by reducing the linewidth of the ISB transition ( 2 ) and by minimizing the impact of higher-order electronic states (captured by the effective decay rate   ).In order to test these predictions, we optimized a MQW system based on Ga0.53In0.47As/GaAs0.51Sb0.49heterostructures to produce ISB transitions with extremely narrow linewidths, more than a factor of 4 narrower than those used in previous works [19].These new MQW materials allowed us to experimentally demonstrate nonlinear polaritonic metasurface limiters with absolute reflection contrast of about 54%, while simultaneously displaying record-high values of effective  (3) .The measured values of reflection contrast are partially limited by the range of input powers available in our setup.Specifically, we estimate that if our measurements would include even lower input powers, a reflection contrast as high as 64% could be demonstrated within the same fabricated device.Additionally, we predict that, by using the same MQW material, the reflection contrast can approach values as high as 94% if the absorption induced by the metasurface could be removedfor example, by working with fully-dielectric designs.
Besides boosting the reflection contrast, the narrow linewidth of the ISB transitions also allowed us to observe fine-structure features in the metasurface absorption spectra that are induced by the bending of the semiconductor bands close to the metal contacts, resulting in additional ISB transitions with different transition frequencies.These additional states are crucial to correctly model the linear and nonlinear dynamics of polaritonic metasurfaces featuring giant nonlinearities.We anticipate that the properties of these additional ISB transitions may be tailored in future studies in order to further improve the metasurface performance as new degrees of freedom.
Overall, the results presented in this work provide an important step forward towards the design of efficient nonlinear polaritonic metasurfaces for application as saturable absorbers, power limiters, and nonlinear elements for all-optical control.

Experimental Setup
The linear reflection spectra in Fig. 2 were acquired with a commercial FTIR spectrometer (Bruker INVENIO) coupled with a Hyperion microscope equipped with a low-NA ZnSe objective lens (Edmund, LFO-5-18).The spectra were normalized with respect to the reflection from the nearby gold substrate.The nonlinear reflection measurements (Fig. 3) were acquired with the custom-built setup shown in Fig. 5.The sample is excited with a pulsed laser (pulse duration  = 2 ps, repetition rate f = 80 MHz) with central frequency tunable within the mid-IR.The mid-IR signal is obtained by difference-frequency generation module (APE, HarmoniXX DFG) which is fed by the signal and idler produced by an optical parametric oscillator (APE, Levante IR ps), which is in turn excited by a Yb pump laser (APE, Emerald Engine).A couple of holographic wire grid polarizers (LP1 and LP2, Thorlabs, WP25H-B) were used to control the power and polarization of the laser.The second polarizer (LP2) is set such that polarization of the laser is parallel to the long side of the antennas.The laser is mechanically chopped at 500 Hz, and then directed to a 50:50 ZnSe beamsplitter (Thorlabs, BSW710).Half of the signal impinges on first HgCdTe detector (MCT1, Thorlabs, PDAVJ10).The other half of the signal is directed towards the sample.The signal is focused on the sample with a geltech aspheric lens (Lexc) with focal length f = 5.95 mm (Thorlabs, C028TME-F).The reflected signal is collected via the same lens and measured via a second (identical) HgCdTe detector (MCT2).Several neutral-density filters were used in front of both detectors to avoid saturation.Each detector is connected to a dedicated lock-in amplifier, and both lock-in amplifiers are referenced by the chopper controller.In order to measure the reflection level of the metasurface versus impinging power, the first polarizer (LP1) was systematically rotated in equal steps via a stepper motor (Thorlabs, K10CR1), and the corresponding signal generated by the two detectors was recorded.This procedure was repeated twice, first with the excitation beam focused on the metasurface and then with the beam focused on the bare substrate.By comparing the signal read by MCT2 in the two cases (and using the signal from MTC1 to monitor the laser power), the reflection level of the metasurface (with respect to the bare substrate) is then obtained.In order to measure the absolute power and intensity levels impinging on the sample, an additional calibration run was performed by measuring the average input power  avg () with a thermal powermeter (Thorlabs, S401C) placed before the lens Lexc.Due to the duty cycle of the chopper (50%), the actual input average power is   = 2 ×  avg () .Following standard formulas, the impinging peak intensity (both in time and space) was then calculated by where f = 80 MHz and  = 2 ps.The spot radius w0 (defined as the distance from the center of the spot at which the intensity drops by a factor  2 ) was determined by both knife-edge measurements and by imaging the laser spot on a bolometric camera and comparing it with features of known sizes on the sample.We found that the spot radius is almost diffraction limited, and it increases linearly with the wavelength following the relation 0 w 0.73  .

Sample Design and Fabrication
To fabricate the metasurfaces, a semiconductor heterostructure containing the MQWs and buffer layers is grown on top of a semi-insulating InP substrate using molecular beam epitaxy (MBE).The grown wafer is then coated with metal, thermocompressively bonded to the host wafer, stripped of the original InP substrate, and finally patterned into nanoantennas with the metal coating evaporated on top of the nanoantennas to form the double-metal resonators shown in Fig. 1(a).
The details of the heterostructure, including the substrate, are as follows, starting from the top layer of the original sample produced in the MBE: 1. Bottom GaInAs buffer layer, n-doping 1 ⋅ 10 19 cm −3 , 5 nm thick 2.
MQW structure, the structure is described in the main text, total thickness is 346 nm 3.
300 µm semi-insulating InP substrate A metal layer composed of a titanium, platinum, and gold (layers thicknesses of 10 nm, 80 nm, 300 nm, respectively) is evaporated on top of the wafer.The wafer is then bonded episide-down to a similarly metalcoated n-doped InP host wafer (bonding temperature of 300°C, bonding pressure 1 bar).Subsequently, the substrate (layer 6) is removed using a concentrated HCl solution, followed by the selective removal of the etch-stop layers 5 and 4 using hydrogen peroxide-phosphoric acid solution mix and a diluted HCl solution, respectively.The 200µm x 200µm nanoantenna array patterns for each metasurface variant is lithographically defined using an electron beam system (Raith eLINE 30 kV) and a 100 nm-SiN hard mask.The MQW layer was then etched with Cl-Ar plasma in an inductively-coupled reactive ion etch (ICP RIE) system and the hard mask was removed using a CF4-based RIE plasma.Finally, a 3 nm thick titanium and a 50 nm thick top gold layer was evaporated.

Figure 1 .
Figure 1.(a) Schematic of the polaritonic metasurface.An incoming field ( + ) feeds a nanoantenna, which in turn couples to the ISB transitions of the MWQ stack with a Rabi rate , and determines the outgoing field ( − ).(b) In the strong coupling regime, the low-power reflection spectrum of the system (blue line) displays the characteristic polaritonic splitting, leading to very large reflection at  =  0 .For high impinging intensities (red line) the ISB transitions in all quantum wells are fully saturated, and the reflection spectrum approaches the one of the bare nanoantenna, featuring very low reflection at  =  0 .In these calculations we assumed   =   =  2 = Ω/10.(c) Reflection spectra of the metasurface versus the saturation level of the MQW for two different impinging frequencies:  =  0 (orange line in panel (c), also marked in panel (b) with the orange arrow) and  =  0 + Ω (green line in panel (c), also marked in panel (b) with the green arrow).(d) Reflection contrast Δ defined as the difference in the low-and high-intensity reflection coefficients of the metasurface versus  2 /Ω, for  =  0 .In this plot, the value of  2 is varied while all other parameters are the same as in panel b.(e) Measured absorption spectra (black line) of the bare InGaAs/GaAsSb MQW material.By fitting the main peak with an analytical model (red dashed line) we extract a transition energy of about 153 meV and a linewidth equal to FWHM = 2ℏ 2 = 7.4 meV.

Figure 2 .
Figure 2. (a) SEM image of the metasurface, scalebar = 3 µm.(b) Zoomed-in SEM image, showing a single unit cell, scalebar = 500 nm.The lattice pitches are  x = 2 μm and  y = 4 μm and the antenna's dimensions are   ≈ 330 nm,   ≈ 1845 nm and the MQW stack thickness is  ≈ 350 nm.(c) Measured y-polarized reflection spectra of several devices as a function of   (all other dimensions are as in panels a-b).The horizontal dashed-dotted blue line indicates the frequency of the bare ISB transition, while the dashed white line indicates the resonant frequency of the bare antenna.(d) Measured reflection spectra of the device with   ≈ 1845 nm (blue solid line), corresponding to vertical dashed line in panel c.The other lines show the low-power reflection spectra calculated either without (black dashed line) or with (red line) the Schottky effect.(e) Calculated low-power reflection spectra of the metasurface, for the same geometries and frequency as in panel c, including the band structure bending caused by the Schottky effect.(f) Numerically calculated band structure of the MQW, including the Schottky effect at the top MQW-titanium/gold interface.

Figure 3 .
Figure 3. Nonlinear response of the metasurface.(a) Measured low-power reflection spectra of the device selected for the nonlinear measurements (L  ≈ 1845 nm), corresponding to the vertical dashed line in Fig. 2b.(b) Measured reflection versus impinging peak intensities at different selected frequencies.The excitation frequency in each measurement is color-coded (see colorbar) and indicated by a corresponding vertical line in panel a.(c) Experimentally measured (blue dots) and numerically calculated (orange line) reflection contrast versus impinging frequency.(d) Measured reflection versus impinging peak intensity and frequency, for a larger number of impinging frequencies.(e) Numerically calculated reflection, corresponding to the experimental data shown in panel d.(f) Comparison between measured (filled circles) and calculated (dashed lines) nonlinear reflection curves for two representative frequencies.

Figure 4 .
Figure 4. (a) Absolute values of |/| extracted from measurements (filled blue circles) and calculations (dashed line) versus the impinging frequency, at a peak intensity of I ≈ 200 KW/cm 2 (horizontal dashed line in panel b).The linear spectrum of the device is also shown for reference (solid gray line, right vertical axis).(b) Calculated absolute values of the derivative |/|, versus impinging frequency and peak intensity.The horizontal dashed line denotes the intensity considered in panel a.

Figure 5 .
Figure 5. Experimental setup used to measure nonlinear reflection.