# Cavity modification of magnetoplasmon mode through coupling with intersubband polaritons

Lucy L. Hale, Daniele de Bernardis, Stephan Lempereur, Lianhe Li, A G Davies, Edmund Linfield, T. Blaikie, Chris Deimert, Zbigniew Wasilewski, Iacopo Carusotto, et al.

► To cite this version:

Lucy L. Hale, Daniele de Bernardis, Stephan Lempereur, Lianhe Li, A G Davies, et al. Cavity modification of magnetoplasmon mode through coupling with intersubband polaritons. 2025. hal-05327842

HAL Id: hal-05327842

https://hal.science/hal-05327842v1

Preprint submitted on 23 Oct 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Cavity modification of magnetoplasmon mode through coupling with intersubband polaritons

Lucy L. Hale¹,* Daniele De Bernardis², Stephan Lempereur¹, Lianhe H. Li³, A. Giles Davies³,
Edmund H. Linfield³, Trevor Blaikie⁴, Chris Deimert⁴, Zbigniew R. Wasilewski⁴, Iacopo Carusotto⁵,
Jean-Michel Manceau⁶, Mathieu Jeannin⁶, Raffaele Colombelli⁶, Jérôme Faist¹, and Giacomo Calari¹

¹ ETH Zürich, Institute of Quantum Electronics,
Auguste-Piccard-Hof 1, 8093 Zürich, Switzerland

²National Institute of Optics (CNR-INO), c/o LENS via Nello Carrara 1, Sesto F.no 50019, Italy

³School of Electronic and Electrical Engineering,
University of Leeds, Leeds LS2 9JT, United Kingdom

⁴Department of Electrical and Computer Engineering, University of Waterloo,
200 University Ave W, Waterloo, Ontario N2L 3G1, Canada

⁵ Pitaevskii BEC Center, CNR-INO and Dipartimento di Fisica,
Università di Trento, I-38123 Trento, Italy and

⁶Centre de Nanosciences et de Nanotechnologies (C2N),
CNRS UMR 9001, University of Paris-Saclay, Palaiseau 91120, France

(Dated: October 23, 2025)

We investigate the coupling of a multi-mode metal-insulator-metal cavity to a two-dimensional electron gas (2DEG) in a quantum well in the presence of a strong magnetic field. The TM cavity mode is strongly hybridized with an intersubband transition of the 2DEG, forming a polaritonic mode in the ultrastrong coupling regime, while the TE mode remains an almost purely cavity mode. The magnetoplasmon excitation emerging from the presence of the magnetic field couples with both TM and TE modes, exhibiting different coupling strengths and levels of spatial field inhomogeneity. While the strong homogeneity of the bare TE mode gives rise to the standard anticrossing of strong coupling, the inhomogeneous polaritonic TM mode is shown to activate an observable Coulombic effect in the spectral response, often referred to as non-locality. This experiment demonstrates a cavity-induced modification of the 2DEG response and offers a new route to probing the effect of Coulomb interactions in ultrastrongly coupled systems via reshaping of their cavity mode profiles.

# I. INTRODUCTION

Quantum wells (QWs) embedded in terahertz (THz) metal-insulator-metal (MIM) cavities are a powerful platform for exploring strong light-matter interactions, with high coupling strengths enabled by large dipole moments, collective enhancement, and deeply subwavelength field confinement [1–5]. Typically, the cavity modes are coupled to QW intersubband transitions (ISBTs), where the coupling strength can reach the ultrastrong coupling regime [6–13] characterized by a light-matter interaction strength comparable to the bare resonance frequencies. This so-called ultrastrong coupling regime holds the promise of bringing quantum effects and technologies into the far-infrared range of the spectrum. Examples include the predicted emission of correlated photon pairs through abrupt modulation of the system [4], or the electrical probing of quantum phase transitions [14].

Along a different line, significant efforts have focused on exploiting the bosonic nature of ISB polaritons to develop inversionless lasers [15]. Notable demonstrations include ISB polariton-LO phonon scattering [16, 17], and more recently, final-state stimulation in polariton-polariton scattering schemes [18]. With an eye towards applications, there has also been a surge of activity

aimed at toggling the system in and out of the strong coupling regime via electrical or optical modulation. Such dynamic control is particularly valuable for implementing amplitude or phase modulation on a continuous-wave (CW) carrier at GHz frequencies, as demonstrated in [19–22], and has been recently proposed as a means to realize ultrafast saturable absorber mirrors [23, 24].

While the ISB-MIM cavity approach relies on fixed cavity and material configurations to engineer the light-matter coupling, magnetic field-induced cyclotron resonances offer dynamic, in-situ tunability [3, 25]. Landau polaritons are formed by strongly coupling light to the collective cyclotron resonance – known as a magnetoplasmon (MP). These systems have enabled the highest recorded coupling strengths, as well as allowing direct tuning of the frequency the magnetoplasmon with applied magnetic field [26, 27].

In Landau polariton systems, the description of the coupling relies on Kohn's theorem, which states that the cyclotron resonance frequency in a translationally-invariant two-dimensional gas is unaffected by electron-electron interactions [28]. Despite its simplicity, this result turns out to be extremely powerful, reducing the system to a non-interacting one, and allowing for the development of a simple, intuitive picture of light-matter interactions in terms of a Hopfield-like polaritonic description [29, 30]. At the same time, Kohn's theorem is also a

<footnote>* luhale@phys.ethz.ch</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. Metal-Insulator-Metal Cavities: a) Schematic showing the sample, indicating multi quantum well layer (red), magnetic (B) field direction, incident THz field and intersubband ($P_{ISB}$) and MP ($P_{MP}$) resonances in the QW. Below: side profile of cavity showing dimensions in µm. b). Finite-element simulation of the reflectance spectra of cavity in the absence of QW, indicating the TM and TE modes. In the following of the work, both modes couple to the MP resonance; in addition, the TM mode is strongly coupled to the ISB resonance. c) Simulated in-plane electric field in a cavity cross-section at frequencies corresponding to the TM and TE modes. Color map shows electric field magnitude in the y-z plane, white arrows show the in-plane field orientation.</figcaption></figure>

limitation, excluding non-linearities and possibly many-body physics from the optical properties of these systems. Going beyond its range of validity is thus becoming an important topic of research [31], which turns out to be surprisingly challenging: besides the use of intense THz sources as in Ref. [32], the fundamental strategy to break Kohn's theorem and observe the effect of Coulomb interactions in the cyclotron resonance is to design the system to explicitly break translational invariance. This can be naturally achieved in resonant cavities with highly sub-wavelength field confinement and inhomogeneous fields [33], where the violation of Kohn's theorem is manifested as a non-local polaritonic response [25, 34, 35].

In this work, we demonstrate the breakdown of Kohn's theorem in a strongly hybridized cavity-electron system, comprising of quantum wells strongly coupled to a metal-insulator-metal (MIM) cavity. In contrast to conventional MIM polariton systems, where the cavity couples only to the intersubband transition (ISBT), we apply a magnetic field, activating also the magnetoplasmon (MP) resonance of the two-dimensional electron gas, in a tripartite MIM-ISBT-MP interplay.

Using THz time-domain spectroscopy, we observe a pronounced renormalization of the magnetoplasmon resonance frequency through electron-electron Coulomb interactions that is activated by the spatial dependence of the polaritonic field arising from the hybridization of the TM cavity mode with the ISBT. Unlike Ref. [25], where non-locality resulted from a deeply sub-wavelength cavity miniaturization design (λ/1000), here the violation of Kohn's theorem occurs with moderate field confinement (λ/60), provided that the resonant mode is sufficiently inhomogeneous. Furthermore, we demonstrate

that by tuning the magnetoplasmon with the magnetic field, it is possible to strongly couple it with a different TE mode, which is mostly homogeneous and in-plane polarized. As a consequence, the TE mode is not hybridized with the ISBT, and thus exhibits purely cavity-like behavior. Moreover, the non-locality disappears and Kohn's theorem is re-established due to a high level of spatial homogeneity.

This device thus allows sampling of the MP degree of non-locality by sweeping the magnetic field and by resonantly selecting cavity/polaritonic modes with varying homogeneity, largely increasing the flexibility of previous standard sub-wavelength cavity miniaturization designs [25, 33]. Interestingly, the resonant control on the non-local response of strongly hybridized cavity-electron systems is purely cavity-induced, effectively realizing a semi-classical counterpart of the modification of materials induced by the cavity quantum vacuum [36, 37].

## II. EXPERIMENTAL SET-UP

### A. Cavity composition

The MIM cavity design is shown in Figure 1a. The QWs (red) are sandwiched between a gold back plane and a gold grating with a periodicity of 40 µm and a duty cycle of 75%. The cavity supports a variety of modes when excited at normal incidence with THz light polarized perpendicular to the grating (shown in Fig. 1a). The first order transverse-magnetic (TM) and transverse-electric (TE) modes are shown in Fig. 1b,c. The TM mode
<figure><img src="image_3.png" /><figcaption>FIG. 2. a) Schematic of experimental reflection THz-TDS set up. THz transmitter and receiver indicated with a 'T' and 'R' respectively. b) Measured reflectance spectra of the QW-loaded cavity as a function of temperature (with no applied magnetic field). The lower polariton branch (LP) is indicated. Below schematic indicating the frequency mismatch between the cavity and ISB resonances. Inset - single spectra slice at 10K. c) Measured THz-TDS transmission spectra of semiconductor heterostructure in the absence of a cavity at 3K. The black dashed line indicates the unperturbed MP frequency given by $\omega_B = eB/m^*$, where $m^* = 0.067m_e$. Inset - single spectra slice at 3T.</figcaption></figure>

sits at 2.4 THz and exhibits a strong electric field amplitude perpendicular to the growth plane (z-direction in 1a) underneath the gold grating. In the gaps, however, the TM mode also supports considerable in-plane fields (y-direction in 1a). In addition, the TE mode at 3.8 THz also lies within the spectral bandwidth of our measurement system. This has primarily in-plane fields, which are strongest underneath the gold grating.

## B. Quantum Well

A 3µm stack of 53 GaAs/AlGaAs square quantum wells (red) with an electron density per well of $3 \times 10^{11} \text{cm}^{-2}$ is placed in the cavity below an 8µm GaAs spacer. The QWs exhibit two resonances in the THz region, which we investigate in our study: the ISBT and the MP resonance.

The ISBT is at 2.74 THz and strongly couples to the TM cavity mode due to its strong z-polarized fields. Figure 2b shows the THz-TDS reflection spectra measured as a function of temperature in the absence of a magnetic field (experimental set-up schematic shown in Fig. 2a). At 150 K, a single spectral feature around 2.4 THz is resolved, corresponding to the TM cavity resonance. As the temperature is reduced, the ISBT becomes well-resolved and couples to the TM mode, resulting in the formation of polariton branches. The difference in spectral contrast observed between upper and lower polariton branches is a result of the mismatch in frequency between the cavity TM resonance (2.4 THz) and the ISBT (2.74 THz).

When a magnetic field is applied, the QW also exhibits an MP resonance due to the collective cyclotron resonance of the electrons [38]. Unlike the ISBT, which is fixed in frequency with the material growth, the MP resonance frequency increases linearly with applied magnetic field as $\omega_B = eB/m^*$, where $m^* = 0.067m_e$ is the effective mass of electrons in the GaAs QW. The MP resonance measured in the bare QW heterostructure without the cavity is shown in Fig. 2c. The THz response shows an almost perfect agreement with the expected linear dependence on the magnetic field (shown as a dashed black line). The linewidth can be estimated to have an almost B-independent value of $\kappa/(2\pi) \approx 300 \text{ GHz}$ (see inset), which is the expected order of magnitude considering the number of quantum wells and electron density [39]. We will see later that this clear agreement of the observed MP frequency with the expected $\omega_B$ is in contrast to the cavity case, where the cavity modes play a prominent role in reshaping the MP optical response.

## III. THZ MAGNETOSPECTROSCOPY OF STRONGLY COUPLED CAVITY

The cavity-coupled QW device is then measured by reflection-mode THz-TDS at 3K with an applied DC magnetic field from 0 – 9 T. In this magnetic field range, the MP resonance sweeps across the entire spectral range of the measurement and interacts with the different cavity modes consecutively. The resulting spectra are shown
<figure><img src="image_4.png" /><figcaption>FIG. 3. Interaction of cavity modes with MP resonance through magnetic field tuning: a). Measured reflectance spectra at 3K as a function of magnetic field strength. Dashed line indicates the cyclotron frequency, $\omega_B$. TM (LP) indicates the lower polariton branch of the ISB polariton with the TM mode, TE indicates the TE mode. b) Calculated spectra using the formalism in Section IV A, accounting for Coulombic interaction. For the theory parameters, see Table I in App. C. Both TM upper (UP) and lower (LP) polaritons are visible, as well as TE mode. c) Close-up of non-smoothed spectra of lower polariton branch. In (a,c), the black dashed line indicates the bare MP resonance as illustrated in Fig.2(c). d) Individual spectra of TM LP for three values of the magnetic field corresponding to $\omega_B < \omega_{TM}$ (purple), $\omega_B = \omega_{TM}$ (orange) and $\omega_B > \omega_{TM}$ (yellow).</figcaption></figure>

in Figure 3a. The MP resonance frequency, given by $\omega_B = eB/m^*$, is shown by the dashed line. For both the TE and TM modes, we observe a clear coupling of the cavity mode to the MP; however, the nature of the coupling is different depending on the mode.

The high contrast feature observed in the spectra at 2.4 THz corresponds to the lower ISB polariton, where the TM cavity mode is coupled to the ISBT. The much lower contrast upper polariton is less visible in the spectra, and we therefore restrict our analysis to the lower polariton only. Although the TM mode is primarily polarized orthogonally to the MP, it still supports sufficient in-plane field components to couple to the MP. The observed coupling around 4 - 6 T (shown more clearly in Fig. 3c,d) therefore corresponds to a tripartite coupling between the ISBT and cavity (which form the ISB polariton) with the MP. Given its weak strength, this LP-

MP coupling does not result in the typical anti-crossing shape, but rather in an enhanced broadening of the spectral cavity feature due to its mixing with the broadened MP resonance. Most interestingly, the location of this broadening feature is displaced from the crossing point with the bare MP frequency $\omega_B$, and can be understood as resulting from a significant blue-shift of the MP which is the signature of the non-local effects stemming from the Coulomb interaction. At larger B values, the interaction with the MP is still visible as a weak red-shift of the ISB polariton frequency, quantified by approximately 40 GHz.

On the other hand, the higher-frequency TE mode at 3.5 THz is a purely photonic mode that has no contribution from the ISBT. It contains mostly in-plane fields which directly couple to the MP resonance, resulting in strong coupling and a clear anti-crossing behaviour. In
contrast to the spectral region around the TM mode, here the MP appears to fit well the $\omega_B$ dependency on magnetic field. Fitting the response using a Hopfield model (assuming an MP frequency $\omega_B$) results in a polariton gap of 0.6 THz and therefore normalized coupling of $\eta = 0.04$. We note that the value of the uncoupled cavity mode frequency extracted from the experiment in Fig. 3a is slightly shifted from the value 3.8 THz of ab-initio simulations in Fig. 1b. This is likely due to slightly mismatched cavity dimensions from their nominal values.

The physical understanding of the different coupling of the TM and TE mode to the MP and ISBT resonances in the QW will be discussed in the next following Section IVA.

## IV. DISCUSSION

### A. Cavity-induced Coulomb interactions

In order to fully understand the interactions between the cavity modes and the MP - particularly where this deviates from what is expected - it is necessary to consider the impact of the spatially-resolved inhomogeneous fields inside the cavity.

From the theory of macroscopic dielectrics, it is well known that a polarizable medium is strongly affected by its own internal electric field. This self-interaction is due to the fact that the dipolar constituents of the material interact among each other via the Coulomb force, which is often approximated by the dipole-dipole interaction. When the medium is excited by a radiative source, it develops a macroscopic polarization density vector $\mathbf{P}(\mathbf{r}, t)$, oscillating together with the external drive. This macroscopic motion of charges also gives rise to a macroscopic electric field inside the dielectric [40]:

$$ \begin{aligned} \mathbf{E}_{\text{inside}}(\mathbf{r}, t) &= \mathbf{E}_{\text{bulk}}(\mathbf{r}, t) + \mathbf{E}_{\text{surf}}(\mathbf{r}, t) \\ &= \frac{(\mathbf{r} - \mathbf{r}')}{4\pi\epsilon_0} \left[ \int_V d^3\mathbf{r}' \frac{\vec{\nabla}' \cdot \mathbf{P}(\mathbf{r}', t)}{|\mathbf{r} - \mathbf{r}'|^3} - \int_{\partial V} d\mathbf{S}' \frac{\mathbf{P}(\mathbf{r}', t)}{|\mathbf{r} - \mathbf{r}'|^3} \right], \end{aligned} \quad (1) $$

which splits into bulk and surface contributions.

The energy contribution due to this internal electric field, $-\mathbf{P} \cdot \mathbf{E}_{\text{inside}}$, which is dependent on the dielectric geometry, makes it harder to polarize the medium and thus leads to what is called depolarization shift [6, 41–43]. In most situations, the system is so large that the surface contribution becomes negligible, and the only remaining contribution comes from the inhomogeneous bulk term.

Using a classical coarse-grained approach [44] (as described in detail in the Appendix) the inclusion of such a depolarization effect results in a $\mathbf{k}$-dependent (wavevector dependent) MP frequency in the 2D QW:

$$ \bar{\omega}_B^2(\mathbf{k}) \approx \omega_B^2 + N_{\text{qw}}\omega_P^2\zeta_{\mathbf{k}}. \quad (2) $$

with an explicit dependence on the material’s plasma frequency

$$ \omega_P = \sqrt{\frac{e^2 n_{2D}}{\epsilon m L_c}}. \quad (3) $$

Here $N_{\text{qw}}$ is the number of quantum wells, and the function $\zeta_{\mathbf{k}}$ depends on the specific cavity-QW geometry. Following App. A2, for the simplest cavity geometry, with two infinite metallic parallel plates at distance $L_c$, we obtain the long-wavelength limit:

$$ \zeta_{\mathbf{k}} \approx \frac{L_c k}{2} - \frac{(L_c k)^2}{4}, \quad (4) $$

holding for $L_c k \ll 1$. The first term of Eq. (4) accounts for the free space contribution of Coulomb interactions, while the second term, proportional to $(L_c k)^2$, is the screening correction coming from the metallic boundaries. This phenomenology is also commonly known as non-locality [25].

Without inhomogeneities in a translationally invariant geometry, the MP is excited only on the $\mathbf{k} = 0$ mode so the system does not experience any depolarization shift. In this case the system is in agreement with Kohn’s theorem and is well described by a completely non-interacting system, without any trace of its internal Coulomb forces. However, in our observed configuration, the cavity-ISBT polariton excites the quantum well with a highly inhomogeneous electric field profile (due to its TM mode structure), leading to a non-homogeneous MP oscillating polarization $\vec{\nabla} \cdot \mathbf{P} \neq 0$. The cavity thus induces a visible Coulombic contribution in the MP dynamics.

### B. Modeling the Coupled System

Tuning the magnetic field directly selects the level of inhomogeneity in the MIM cavity through resonant excitation of the different modes coupled to the MP. The choice of cavity mode in light-matter coupled systems is therefore important not only for maximizing the strength of the coupling with the electronic transition, but also for determining the degree of inhomogeneity and, consequently, the breaking of translational invariance of the system.

To demonstrate this, in Figure 3b we calculate the spectra taking into account the effect of Coulomb interaction and resulting wavelength-dependent MP frequency. This calculation is performed with a simple classical linear model, starting from Maxwell equations in the presence of a polarizable material and projecting them on the cavity two-mode subspace (see App. A). While the cavity modes are transverse modes [45], we do not discard the longitudinal component of the electric field, which gives rise to what we call the *Coulomb contribution* mentioned in Eq. (1). As indeed described in the previous section IVA, and detailed in Appendix A-B, it is this
<figure><img src="image_5.png" /><figcaption>FIG. 4. a) Cavity containing parabolic quantum well heterostructure: schematic (above) and simulated cavity fields (below) b) Experimentally measured spectra with magnetic field for cavity in a. Arrows indicate positions of line cuts in c) for lower polariton (bottom) and upper polariton (top). c) Individual spectra of TM LP (left) and UP (right) for three values of the magnetic field corresponding to the purple, orange and yellow arrows in Fig. 4a. d) Calculated transmission of sample in b. For the theory parameters, see Table I in App. C.</figcaption></figure>

longitudinal contribution which provides a *k*-dependent depolarization shift affecting the MP during the interaction with the TM mode.

We see a good agreement of the measured spectra with the calculated spectra, in Fig. 3a and Fig. 3b, respectively. Comparing more closely in Fig. 3d, the calculated spectra accurately demonstrate the broadened linewidth of the ISB polariton as a result of the interaction with the Coulomb-shifted MP. The TE mode is also well reproduced, showing a clear anti-crossing at the MP frequency, in contrast to the TM mode. This demonstrates that to accurately model light-matter coupled systems with electronic transitions and inhomogeneous electric fields, it is crucial to include the influence of the field inhomogeneity on the Coulomb interactions in the 2DEG.

### C. Generality of Non-Locality

For sufficiently doped quantum wells, the nonlocality behavior only depends on the inhomogeneity of the cavity mode, and can thus be recovered in other types of heterostructure, and for different cavity geometries. To demonstrate this, in Figure 4 we show the case of ISB polaritons in a cavity containing parabolic quantum wells [11]. Here, the ISBT is still coupled to the TM mode of the cavity, but the continuously-graded material composition of the quantum well allows for a lower frequency ISB transition, and therefore larger ISB polariton coupling strength. To match the TM mode frequency to the ISB frequency, the cavity dimensions are slightly different from the square well case, with no GaAs buffer region (Fig. 4a). Despite the altered cavity dimensions, we still observe significant inhomogeneity in the fields, particularly at the interface between the metallic cavity plate and the gap. In addition, the lack of a buffer region means that the electrons in the quantum wells also ex-
experience the fields at the grating edges, which are highly inhomogeneous.

As for the previous samples, Figure 4b shows the measured spectra when a magnetic field is applied, tuning the MP resonance across the ISB polaritons. Here, interaction of the MP with both polariton branches is visible thanks to the closer frequency alignment between the TM cavity mode and ISB, although the spectral features are slightly less clear due to the broadened polariton branches. For both polariton branches, these features include a cavity-induced broadening around the crossing point, a blue-shift of the MP from its bare frequency $\omega_B$ indicated as a black line, and an eventual red-shift of the polariton branches at strong B fields, as seen in the square well case in Figure 3.

By demonstrating the effect of Coulomb interaction in different cavities with two different heterostructures, we highlight that the level of nonlocality can be tailored through cavity design, independently from the quantum well system. To further support this, in Appendix D we show the calculated spectra for both the cavities in Fig. 3 and Fig. 4 without including Coulomb interaction, demonstrating how the experimental spectra is clearly modified in the experiment by nonlocality. Moreover, by further increasing the coupling strength either through different cavity design or increasing the electron density, the Coulomb interaction and thereby nonlocality can be further enhanced.

## V. CONCLUSION

In conclusion, we have experimentally demonstrated how the spatial inhomogeneity of a cavity mode can fundamentally alter the response of a two-dimensional electron gas (2DEG) in the strong light-matter coupling regime in the presence of a strong magnetic field. Using a multi-mode metal-insulator-metal (MIM) cavity containing quantum wells, we access different coupling regimes by tuning the magnetoplasmon resonance as a function of magnetic field and selectively coupling it to distinct cavity or polaritonic modes, the latter resulting from the hybridization of cavity modes with the intersubband transition in the QW.

We find that the MP-cavity coupling in spatially inhomogeneous field configurations results in significant new features, including spectral shifts attributable to Coulomb interaction. This breakdown of translational invariance - induced purely by the spatial profile of the cavity field - demonstrates a method for altering materials through cavity coupling that arises not from vacuum fluctuations but simply from the cavity field spatial profile. Furthermore, the cavity-induced Coulomb interactions can be activated on demand via magnetic field tuning and are controlled through the cavity design.

These findings broaden the understanding of light-matter interactions in the strong coupling regime and provide a versatile route for tailoring electronic

excitations in solid-state systems. Beyond this, the ability to manipulate polariton properties dynamically through cavity mode selection and magnetic field tuning opens new possibilities for tunable polaritonic devices, where the tripartite interplay can be also exploited to study the cavity protection phenomena or the quantum vacuum of ultrastrongly coupled systems [46, 47].

## ACKNOWLEDGMENTS

We thank Gian Marcello Andolina, Bianca Turini, and Alberto Nardin for interesting and insightful discussions. D.D.B. acknowledges funding from the European Union - NextGeneration EU, "Integrated infrastructure initiative in Photonic and Quantum Sciences" - I-PHOQS [IR0000016, ID D2B8D520, CUP B53C22001750006]. IC acknowledges support from the Provincia Autonoma di Trento; from the Q@TN Initiative; from the National Quantum Science and Technology Institute through the PNRR MUR Project PE0000023-NQSTI, co-funded by the European Union - NextGeneration EU. J.F., G.S. and L.H. acknowledge funding through the SNF grant 200021-227521. This work was also supported by the European Union Future and Emerging Technologies (FET) Grant No. 737017 (MIR-BOSE) and partially by the French RENATECH network.

## Appendix A: General theoretical description

Here we develop a fully classical description of the transmittivity/reflectivity of the THz cavity coupled to two-dimensional quantum well (QW). The quantum well is characterized by a polarization density $P(r)$, having the dimensions of dipole moment per volume. The polarization density can be then linearly decomposed in all the contributions due to the independent transitions in the material. While these transitions have fully quantum origin, it will be clear that a classical description fully matches the pure quantum one [4, 30]. Here for simplicity we focus on the case where the polarization is fully given by a MP transition [3], but the generalization to other type of transitions (i.e. intersubband [4]) is straightforward.

While the theory is developed in a fully general way, we will focus here on the simplest case of a planar, translationally invariant cavity filled with a medium having $\epsilon_r = 1$ with no grating.

### 1. Maxwell equation with polarizable matter

We first start by considering Maxwell equations

$$ \vec{\nabla} \cdot \mathbf{E} = \frac{\rho}{\epsilon_0} \qquad (A1) $$
$$ \vec{\nabla} \times \mathbf{B} = \frac{1}{c^2} \left( \frac{\mathbf{J}}{\epsilon_0} + \frac{\partial}{\partial t} \mathbf{E} \right) \quad (A2) $$

$$ \vec{\nabla} \cdot \mathbf{B} = 0 \quad (A3) $$

$$ \vec{\nabla} \times \mathbf{E} = -\frac{\partial}{\partial t} \mathbf{B}. \quad (A4) $$

As a first task in developing our theory we rewrite these equation in a minimal form to treat problems with polarizable matter.

Here, the charge density is fully given by the polarization density of the QW

$$ \rho = -\vec{\nabla} \cdot \mathbf{P}_{3D}, \quad (A5) $$

which is related to the current by the dipole (long wavelength) approximation

$$ \mathbf{J} = \partial_t \mathbf{P}_{3D}. \quad (A6) $$

Here we use the label 3D in $\mathbf{P}_{3D}$ to distinguish from the two-dimensional QW polarization $\mathbf{P}$ which will be introduced later on.

Combining Eq. (A2) with the rotor of Eq. (A4) we obtain

$$ -c^2 \nabla^2 \mathbf{E} + c^2 \vec{\nabla} (\vec{\nabla} \cdot \mathbf{E}) = - \partial_t^2 \mathbf{E} - \frac{1}{\epsilon_0} \partial_t^2 \mathbf{P}_{3D}. \quad (A7) $$

To further reduce the system of equations and make it solvable, we need to take the Gauss law into account

$$ \vec{\nabla} \cdot \mathbf{E} = - \frac{\vec{\nabla} \cdot \mathbf{P}_{3D}}{\epsilon_0}, \quad (A8) $$

so to isolate the so called transverse and longitudinal parts of the electric field

$$ \mathbf{E} = \mathbf{E}^{\parallel} + \mathbf{E}^{\perp}. \quad (A9) $$

here

$$ \mathbf{E}^{\parallel} = \frac{1}{\epsilon_0} \vec{\nabla} (G \star \vec{\nabla} \cdot \mathbf{P}_{3D}) \quad (A10) $$

and

$$ \mathbf{E}^{\perp} = \partial_t \mathbf{A}. \quad (A11) $$

The vector $\mathbf{A}$ is the vector potential, having the property $\vec{\nabla} \times \mathbf{A} = \mathbf{B}$ and $\vec{\nabla} \cdot \mathbf{A} = 0$. As a consequence, $\vec{\nabla} \cdot \mathbf{E}^{\perp} = 0$ by construction, coinciding with the standard definition of transverse vector [45]. We also introduced the Green's function of the Poisson equation, $G$, which solves $-\nabla^2 G(\mathbf{r}, \mathbf{r}') = \delta(\mathbf{r} - \mathbf{r}')$ with metallic boundary conditions on the plates and zero electric potential difference between them. The $\star$ denotes the convolution operator in real space. As a consequence also $\vec{\nabla} \times \mathbf{E}^{\parallel} = 0$

by construction. These definitions are general and true in a cavity setup or any other confined geometry.

We finally arrive at the main equation describing the electric field dynamics in the grating cavity

$$ -c^2 \nabla^2 \mathbf{E}^{\perp} + \partial_t^2 \mathbf{E}^{\perp} = -\frac{1}{\epsilon_0} \partial_t^2 [\mathbf{P}_{3D} + \vec{\nabla}(G \star \vec{\nabla} \cdot \mathbf{P}_{3D})]. \quad (A12) $$

The term inside the square brackets on the left-hand side is the transverse projected polarization density

$$ \mathbf{P}_{3D}^{\perp} = \mathbf{P}_{3D} + \vec{\nabla}(G \star \vec{\nabla} \cdot \mathbf{P}_{3D}), \quad (A13) $$

having the property

$$ \vec{\nabla} \cdot \mathbf{P}_{3D}^{\perp} = 0. \quad (A14) $$

Its longitudinal projection is instead

$$ \mathbf{P}_{3D}^{\parallel} = -\vec{\nabla}(G \star \vec{\nabla} \cdot \mathbf{P}_{3D}). \quad (A15) $$

Evidently $\mathbf{P}_{3D}^{\perp} + \mathbf{P}_{3D}^{\parallel} = \mathbf{P}_{3D}$.

We consider the transverse electric field decomposed on the cavity eigenmodes

$$ \mathbf{E}^{\perp}(\mathbf{r}, z) = \sum_{\lambda} \mathbf{w}_{\lambda}(\mathbf{r}, z) \mathcal{E}_{\lambda}(t). \quad (A16) $$

Here we use $\mathbf{r} = (x, y)$ as the in-plane position, keeping explicit $z$, and $V$ is the system's volume. This convention will be especially convenient when we introduce the two-dimensional quantum well. In particular, $\mathbf{w}_{\lambda}(\mathbf{r}, z)$ is an adimensional function that fully characterises the $\lambda$ electromagnetic eigenmode. The index $\lambda$ is intended as an generalized index, comprising the polarization index and eventual wavevector. It has the properties that $\vec{\nabla} \cdot \mathbf{w}_{\lambda}(\mathbf{r}, z) = 0$ and it satisfy the Helmholtz equation

$$ -c^2 \nabla^2 \mathbf{w}_{\lambda}(\mathbf{r}, z) = \omega_{\lambda}^2 \mathbf{w}_{\lambda}(\mathbf{r}, z). \quad (A17) $$

These functions are orthonormal and normalized such that

$$ \int d^3 r \, \mathbf{w}_{\lambda}(\mathbf{r}, z) \cdot \mathbf{w}_{\lambda'}(\mathbf{r}, z) = \delta_{\lambda\lambda'}, \quad (A18) $$

where $V$ is the system's volume. Projecting both sides of Eq. (A12) on these eigenmodes, we finally obtain the main equation for the amplitude of the cavity modes

$$ \partial_t^2 \tilde{\mathcal{E}}_{\lambda} + \omega_{\lambda}^2 \tilde{\mathcal{E}}_{\lambda} = \partial_t^2 \left[ \int d^2 r \, dz \, \mathbf{w}_{\lambda}(\mathbf{r}, z) \cdot \mathbf{P}_{3D}(\mathbf{r}, z) \right]. \quad (A19) $$

Here $\tilde{\mathcal{E}}_{\lambda} = \epsilon \mathcal{E}_{\lambda}$ and has the same dimensional units of $\mathbf{P}_{3D}(\mathbf{r}, z)\sqrt{V}$.

Importantly in Eq. (A19) the longitudinal polarization
has disappeared because

$$
\begin{align*}
& \int d^2r dz \mathbf{w}_\lambda(\mathbf{r}, z) \cdot \vec{\nabla} [G \star \vec{\nabla} \cdot \mathbf{P}_{3D}](\mathbf{r}, z) = \\
& \int d^2r dz \mathbf{w}_\lambda(\mathbf{r}, z) \cdot \vec{\nabla} [\varphi_\mathbf{P}(\mathbf{r}, z)] \\
& = \int d^2r dz \vec{\nabla} \cdot [\mathbf{w}_\lambda(\mathbf{r}, z) \varphi_\mathbf{P}(\mathbf{r}, z)] \tag{A20} \\
& - \int d^2r dz [\vec{\nabla} \cdot \mathbf{w}_\lambda(\mathbf{r}, z)] \varphi_\mathbf{P}(\mathbf{r}, z) \\
& = 0.
\end{align*}
$$

where

$$
\varphi_{\mathbf{P}}(\mathbf{r}, z) = \int d^2 r' dz' G(\mathbf{r}, z, \mathbf{r}', z') \vec{\nabla}' \cdot \mathbf{P}_{3D}(\mathbf{r}', z'). \quad (\text{A21})
$$

The term in the third line of Eq. (A20) is a boundary
term, vanishing in a infinitely large system. The term in
the fourth line is zero as well by definition of the trans-
verse mode functions $\vec{\nabla} \cdot \mathbf{w}_\lambda(\mathbf{r}, z) = 0$.

2. Long wavelength Green's function

The most important object to compute the longitudi-
nal contribution to the dielectric response is the Green’s
function of the Poisson equation

$$
-\nabla^2 G(\mathbf{r}, z, \mathbf{r}', z') = \delta^{(2)}(\mathbf{r} - \mathbf{r}')\delta(z - z'), \quad (\text{A22})
$$

this gives rise to the well-known Coulomb interaction.

In free space it gives the standard Coulomb potential
$G(\mathbf{r}, z, \mathbf{r}', z') = 1/(4\pi|\mathbf{r}_{3D} - \mathbf{r}'_{3D}|)$, where $\mathbf{r}_{3D} = (\mathbf{r}, z)$,
while in a infinite parallel mirrors cavity one must in-
cludes all the image charges extra terms [48].

For our purposes, we don't actually need the real space
representation, but rather its Fourier k-space one. To
have the maximally accurate description, we would need
to consider the solution relative to the grating patch cav-
ity, with appropriate boundary conditions. In practice,
this is quite challenging, so we restrict our discussion to
the simplified case where two parallel, perfectly metallic
mirrors form the cavity. The distance between the two
mirrors is the cavity length, and is called $L_c$. Introducing
$k = \sqrt{k_x^2 + k_y^2}$, the Poisson equation is rewritten as

$$
(-\partial_z + k^2) G_k(z, z') = \delta(z - z'), \quad (\text{A23})
$$

imposing the metallic boundary conditions $G_k(z = 0, z') = G_k(z = L_c, z') = 0$. The solution is given by

$$
G_k(z, z') = g_k(z, z')\Theta(z - z') + g_k(z', z)\Theta(z' - z), \quad (\text{A24})
$$

where

$$
\begin{equation}
\begin{split}
g_k(z, z') = \frac{e^{-k|z-z'|}}{2k} & \\
& - \frac{1}{2k \sinh(kL_c)} \left[ \sinh(kz)e^{k(z'-L_c)} + \sinh(k(L_c-z))e^{-kz'} \right],
\end{split}
\tag{A25}
\end{equation}
$$

and $\Theta(z)$ is the Heaviside step function.

The long-wavelength expansion for $L_c k \to 0$ gives

$$
\lim_{k \to 0} g_k(z, z') = \frac{1}{2k} - \frac{zz'}{L_c}. \quad (A26)
$$

It is worth stressing that this expansion holds only if $k \ll 1/L_c$. Having a box of size $L_x$, for which $k = 2\pi/L_x$, this immediately implies $L_x \gg L_c$, suggesting that the limit $L_c \to \infty$ does not commute with the longwavelength limit.

Taking $z = z' = L_c/2$, we have

$$
k^2 G_k(L_c/2, L_c/2) \approx \frac{k}{2} - \frac{L_c k^2}{4} \quad (\text{A27})
$$

Instead, to calculate the Green’s function in the free
space limit (no cavity), we have to take the limit $L_c \to \infty$
already in (A25), which recovers the standard Coulomb
dispersion,

$$
k^2 G_k(L_c/2, L_c/2) \approx \frac{k}{2}. \qquad (A28)
$$

3. Cyclotron polarization and non-locality

While the cyclotron (or MP) transition in a 2DEG
(two-dimensional electron gas) is a fully quantum phe-
nomenon, due to harmonicity, it can be equivalently de-
scribed by a classical system. In particular, here we
develop a coarse-grained description of the MP polar-
ization, which is also often dubbed in the literature as
macroscopic dielectric theory.

Let's consider a single electron as a classical charged
particle, confined in a two dimensional plane at $z = 0$,
in a strong homogeneous perpendicular magnetic field
**B**<sub>ext</sub> = (0, 0, B<sub>ext</sub>) and in-plane electric field **E**(r, z) =
(E<sub>x</sub>(r, z), E<sub>y</sub>(r, z), 0). As before we take as a convention
**r** = (x, y) as the in-plane position. Calling δ**r** the in-
plane coordinate of the electron, the equations of motion
are

$$
m \ddot{\mathbf{r}} = e \delta \dot{\mathbf{r}} \times \mathbf{B}_{\text{ext}} + e \mathbf{E}(\delta \mathbf{r}, z=0). \quad (\text{A29})
$$

By integrating out the equation for $\delta y$, we obtain the
equation for $\delta x$ (the same applies to $\delta y$)

$$
\delta \ddot{x} + \omega_B^2 \delta x = \frac{e}{m} E_x(\delta \mathbf{r}, 0), \quad (A30)
$$

which describes an harmonic motion coupled to the elec-
tric field.

We then extend the description to many electrons
$n = 1, 2, \dots, \delta N_e$, all localized around the position **r**. Here
we completely neglect the interaction between these elec-
tron. While this assumption is surely wrong in a real
physical system, here is motivated by the fact that we
can interpret $\delta \mathbf{r}_n$ as the coordinate of effective degrees
of freedom, resulting from a proper quantum treatment
where the Fermionic nature of electrons allows to remap
the interacting problem in an independent particle model [49].

We take a long-wavelength approximation approximating the electric field as homogeneous around that position $\mathbf{E}(\delta\mathbf{r}_1, \delta\mathbf{r}_2, \dots, \delta\mathbf{r}_N) \approx \mathbf{E}(\mathbf{r})$. We identify the material total dipole moment at this fixed position as $e \sum_n \delta\mathbf{r}_n$. From Eq. (A30) we obtain the equation of the oscillating total dipole

$$ \partial_t^2 \sum_n \delta x_n + \omega_B^2 \sum_n \delta x_n = \frac{e \delta N_e}{m} E_x(\mathbf{r}). \quad (A31) $$

Here, the displacements $\delta x_n$ are relative to the fixed position $\mathbf{r}$, which can be interpreted as the center of mass position of this electronic ensemble. In general, we now understand that the macroscopic two-dimensional polarization can be defined as

$$ \mathbf{P}(\mathbf{r}) = e \sum_n \frac{\delta \mathbf{r}_n}{\delta S} \quad (A32) $$

where $\delta S$ is a small surface element containing the $\delta N_e$ electrons [44]. A similar treatment applies to the more general three-dimensional polarization $\mathbf{P}_{3D}$ defined in the previous subsection, but here, for simplicity, we focus directly on the specific case of our interest. Since the general electrodynamics equations developed in Sec. A1 require a three-dimensional polarization, it is worth noticing that the two-dimensional one can be casted to a three-dimensional one via delta function

$$ \mathbf{P}_{3D}(\mathbf{r}, z) = \mathbf{P}(\mathbf{r})\delta(z). \quad (A33) $$

We then arrive at the main equation describing the harmonic motion of the in-plane polarization

$$ \partial_t^2 \mathbf{P}(\mathbf{r}) + \omega_B^2 \mathbf{P}(\mathbf{r}) = \omega_P L_c \epsilon \mathbf{E}(\mathbf{r}, 0). \quad (A34) $$

Here we have introduced the plasma frequency as

$$ \omega_P = \sqrt{\frac{e^2 n_{2D}}{\epsilon m L_c}}, \quad (A35) $$

where we have used the three-dimensional electron density as $n_{3D} = \delta N_e / \delta V$, leading to the two-dimensional density $n_{2D} = n_{3D} L_c$. $\epsilon = \epsilon_r \epsilon_0$ is the dielectric permittivity of the considered material while $L_c$ is the cavity height. From here on the system's volume is given by $V = L_c S$. A full quantum treatment leads to a different coupling constant, which depends explicitly on the filling factor $\nu$ [30]. However, our classical description gives very similar quantitative results without the complexity of a full quantum description and with the advantage of having a more transparent interpretation.

In light of the development in the previous subsection, we further manipulate Eq. (A34) to account for the split between transverse and longitudinal electric field components. We use $\mathbf{E}(\mathbf{r}, 0) = \mathbf{E}^{\perp}(\mathbf{r}, 0) + \vec{\nabla}[G \star \vec{\nabla}]$.

$\mathbf{P}_{3D}](\mathbf{r}, 0)/\epsilon = \mathbf{E}^{\perp}(\mathbf{r}, 0) - \mathbf{P}^{\parallel}(\mathbf{r})/\epsilon$, which is a consequence of Eqs. (A10)-(A15). The main equation of motion for the cyclotron polarizability is then

$$ \partial_t^2 \mathbf{P}(\mathbf{r}) + \omega_B^2 \mathbf{P}(\mathbf{r}) = \omega_P^2 (\epsilon L_c \mathbf{E}^{\perp}(\mathbf{r}, 0) - \mathbf{P}^{\parallel}(\mathbf{r})). \quad (A36) $$

The two-dimensional version of the longitudinal polarization appearing in this equation is rewritten explicitly as

$$ \mathbf{P}^{\parallel}(\mathbf{r}) = L_c \int d^2 r' \vec{\nabla} \vec{\nabla}' G(\mathbf{r}, z=0, \mathbf{r}', z'=0) \cdot \mathbf{P}(\mathbf{r}'). \quad (A37) $$

Here we have integrated by parts and discarded the boundary terms.

Eq. (A38) is particularly important because it shows how the Coulomb force enters the polarization dynamics through the longitudinal part of the field. This term is what is often referred to as non-locality [25, 31, 50], and we will see that it is the major factor responsible for depolarization shift phenomena in non-homogeneous MP excitations. To this aim it is worth to express Eq. (A38) in Fourier space, by considering $\tilde{\mathbf{P}}_{\mathbf{k}} = \int d^2 r / S e^{-i\mathbf{k} \cdot \mathbf{r}} \mathbf{P}(\mathbf{r})$ (same for $\mathbf{E}^{\perp}(\mathbf{r}, 0)$). We have that

$$ \partial_t^2 \tilde{\mathbf{P}}_{\mathbf{k}} + \omega_B^2 \tilde{\mathbf{P}}_{\mathbf{k}} = \omega_P^2 (\epsilon L_c \tilde{\mathbf{E}}_{\mathbf{k}} - L_c \mathbf{k} \tilde{\mathbf{G}}_{\mathbf{k}} \mathbf{k} \cdot \tilde{\mathbf{P}}_{\mathbf{k}}). \quad (A38) $$

Here we have assumed that the system is translation invariant in the plane, so that $G(\mathbf{r}, z=0, \mathbf{r}', z'=0) = G(\mathbf{r} - \mathbf{r}', z=0, z'=0)$. Restricting for simplicity the system only the $x$-dimension and considering that $\lim_{k \to 0} \tilde{\mathbf{G}}_{\mathbf{k}} = 1/(2k)$ [43, 51] we have

$$ \partial_t^2 \tilde{P}_k + \tilde{\omega}_B^2(k) \tilde{P}_k = \omega_P^2 \epsilon L_c \tilde{E}_k, \quad (A39) $$

where the cyclotron frequency is now shifted by the $k$-dependent depolarization shift [25]

$$ \tilde{\omega}_B^2(k) = \omega_B^2 + \omega_P^2 \frac{L_c k}{2}. \quad (A40) $$

More generally, we can write

$$ \tilde{\omega}_B^2(k, z) = \omega_B^2 + \omega_P^2 \zeta_k(z), \quad (A41) $$

where

$$ \zeta_k = \lim_{k \to 0} k^2 \tilde{G}_k(z = z'). \quad (A42) $$

### 4. Polaritonic transmission in the presence of non-locality

Combining Eq. (A19) with Eq. (A38), we can finally address the coupled light-matter system, where the cyclotron MP excitation is hybridized to the polaritonic/cavity modes to form new hybridized modes.

In order to account for the polaritonic origin of the TM mode we must also include the ISBT polarization,
splitting the total polarization in Eq. (A19) into $\mathbf{P}_{3D} = \mathbf{P}_{3D}^{\text{mp}} + \mathbf{P}_{3D}^{\text{isb}}$ [47].

Again, restricting to only the x-direction, and taking a Fourier transform with respect to time, we obtain the coupled modes equations

$$ (\omega_{\lambda}^{2} - \omega^{2}) \tilde{\mathcal{E}}_{\lambda} = \omega^{2} \sum_{k} (\tilde{u}_{\lambda,k} \tilde{P}_{k}^{\text{mp}} + \tilde{u}_{\lambda,k} \tilde{P}_{k}^{\text{isb}}), \quad (\text{A43}) $$

$$ (\bar{\omega}_{B}^{2}(k) - \omega^{2}) \tilde{P}_{k}^{\text{mp}} = \omega_{P}^{2} \sum_{\lambda} \tilde{w}_{\lambda,k} \tilde{\mathcal{E}}_{\lambda}, \quad (\text{A44}) $$

$$ (\omega_{\text{isb}}^{2} - \omega^{2}) \tilde{P}_{k}^{\text{isb}} = \omega_{P}^{2} \sum_{\lambda} \tilde{u}_{\lambda,k} \tilde{\mathcal{E}}_{\lambda}, \quad (\text{A45}) $$

Here we have shifted the electric field in Eq. (A16) to

$$ \tilde{\mathcal{E}}_{\lambda} \mapsto \tilde{\mathcal{E}}_{\lambda} \sqrt{V}/L_c. \quad (\text{A46}) $$

In this way, the coupling matrices are given by the dimensionless quantities (by restoring the full two-dimensionality dependence)

$$ \tilde{w}_{\lambda, \mathbf{k}} = \sqrt{V} \int \frac{d^2 r}{S} e^{-i \mathbf{k} \cdot \mathbf{r}} \mathbf{w}_{\lambda}(\mathbf{r}, 0) \cdot \mathbf{u}_{x}, \quad (\text{A47}) $$

$$ \tilde{u}_{\lambda, \mathbf{k}} = \sqrt{V} \int \frac{d^2 r}{S} e^{-i \mathbf{k} \cdot \mathbf{r}} \mathbf{w}_{\lambda}(\mathbf{r}, 0) \cdot \mathbf{u}_{z}, \quad (\text{A48}) $$

where $\mathbf{u}_x = (1, 0, 0)$, $\mathbf{u}_z = (0, 0, 1)$ are the x or z unit pointers. Moreover, we have introduced the ISBT frequency, which is assumed to be k-independent and not affected by non-locality.

The transmission/reflection can then be obtained through the formula [46, 47]

$$ \mathcal{T} = \sum_{\lambda} \gamma_{\lambda} \omega_{\lambda} \sum_{\lambda'} \mathcal{M}_{\lambda\lambda'}^{-1} A_{\lambda'}. \quad (\text{A49}) $$

Here we introduced the photon losses of each mode $\gamma_{\lambda}$, which is incorporated in Eqs. (A43)-(A44) by shifting the poles $\omega_{\lambda}^{2} - \omega^{2} \mapsto \omega_{\lambda}^{2} - \omega^{2} - i\gamma_{\lambda}\omega$ for each mode. The dynamical matrix is defined as

$$ \mathcal{M}_{\lambda\lambda'} = (\omega_{\lambda}^{2} - \omega^{2}) \delta_{\lambda\lambda'} - \omega^{2} \omega_{P}^{2} \sum_{k} \frac{\tilde{u}_{\lambda,k} \tilde{u}_{\lambda',k}}{\omega_{\text{isb}}^{2} - \omega^{2}} - \omega^{2} \omega_{P}^{2} \sum_{k} \frac{\tilde{w}_{\lambda,k} \tilde{w}_{\lambda',k}}{\bar{\omega}_{B}^{2}(k) - \omega^{2}}. \quad (\text{A50}) $$

Here we also introduce losses: for the cyclotron excitation, $\kappa$, by substituting $\bar{\omega}_{B}^{2}(k) - \omega^{2} \mapsto \bar{\omega}_{B}^{2}(k) - \omega^{2} - i\kappa\omega$, for the cavity, $\gamma$, by substituting $\omega_{\lambda}^{2} - \omega^{2} \mapsto \omega_{\lambda}^{2} - \omega^{2} - i\gamma\omega$, for the ISBT $\kappa_{\text{isb}}$, by substituting $\omega_{\text{isb}}^{2} - \omega^{2} \mapsto \omega_{\text{isb}}^{2} - \omega^{2} - i\kappa_{\text{isb}}\omega$. The adimensional amplitudes $A_{\lambda}$ are instead given by the projections of the in-plane external drive amplitude $\mathbf{I}(\mathbf{r})$ on the cavity eigenmodes

$$ A_{\lambda} \sim \frac{\int d^2 r \mathbf{w}_{\lambda}(\mathbf{r}, z = L_c) \cdot \mathbf{I}(\mathbf{r})}{\sqrt{\int d^2 r |\mathbf{I}(\mathbf{r})|^2}} \quad (\text{A51}) $$

<figure><img src="image_6.png" /><figcaption>FIG. 5. Schematic of the TM-modes spatial field polarization.</figcaption></figure>

## Appendix B: Three-modes cavity theory

Here we specialize the transmissivity/reflectivity theory developed in Sec. A to describe the experimental setup described in the main text, with the parameters derived in Sec. C.

We restrict the description to only the main two TM, TE cavity modes described before, and we can neglect the coupling between different cavity modes and essentially only keep the diagonal terms of $\mathcal{M}$ in Eq. (A50). Notice that the TM mode is a polaritonic mode, split into lower-polariton (LP) and upper-polariton (UP), so in the end we will see three modes in the transmission, TM$_{LP}$, TM$_{UP}$, TE, all hybridized with the MP. The transmissivity can then be written as

$$ \mathcal{T} \approx \sum_{\lambda} \frac{\gamma_{\lambda} \omega_{\lambda} A_{\lambda}}{\omega_{\lambda}^{2} - \omega^{2} \epsilon_{\lambda}(\omega) - i \gamma_{\lambda} \omega} \quad (\text{B1}) $$

where we introduced the ISBT and MP relative permitivities per mode $\epsilon_{\lambda}(\omega) = 1 + \chi_{\lambda}^{\text{isb}}(\omega) + \chi_{\lambda}^{\text{mp}}(\omega)$, determined by the ISBT and MP susceptibility

$$ \chi_{\lambda}^{\text{isb}}(\omega) = \omega_{P}^{2} \sum_{k} \frac{|\tilde{u}_{\lambda,k}|^{2}}{\omega_{\text{isb}}^{2} - \omega^{2} - i\kappa_{\text{isb}}\omega}, \quad (\text{B2}) $$

$$ \chi_{\lambda}^{\text{mp}}(\omega) = \omega_{P}^{2} \sum_{k} \frac{|\tilde{w}_{\lambda,k}|^{2}}{\bar{\omega}_{B}^{2}(k) - \omega^{2} - i\kappa\omega}. \quad (\text{B3}) $$

For the TM mode, we assume that the eigenmode $x$-projection is a box function, which is constant corresponding to the grating gap between the metallic patches, and zero in correspondence to the patch

$$ \mathbf{w}_{\lambda=\text{TM}}(x, y, 0) \cdot \mathbf{u}_x \approx \sqrt{\frac{f}{2\Delta x}} \Pi(x) \cos(k_1 x). \quad (\text{B4}) $$

Here $\Pi(x) = [1 - \Theta(x - \Delta x/2)]\Theta(x + \Delta x/2)$ is the box function and $\Theta(x)$ is the Heaviside step function. The
adimensional factor $f < 1$ accounts for the non-optimal overlap between the QW and the cavity modes. Moreover, the wavevectors $k_1$ are given by the periodicity of the structure, and we have that $k_1 \approx 2\pi/L_x$. This approximation is schematically represented in Fig. 5. Here, $L_x$ is the x-size of the quantum well. Notice that the normalization factor $\sqrt{1/\Delta x}$ follows from the normalization condition of the eigenmodes in Eq. (A18), taking the y-projection vanishing and the z-projection complementary to the x-one. Their Fourier transform gives

$$ \tilde{w}_{\lambda=\text{TM},k} = \sqrt{f \frac{\Delta x}{8L_x}} \left[ \begin{aligned} & \text{sinc}\left((k - k_1)\frac{\Delta x}{2}\right) \\ & + \text{sinc}\left((k + k_1)\frac{\Delta x}{2}\right) \end{aligned} \right] \quad (B5) $$

Here $\text{sinc}(kx) = \sin(kx)/(kx)$. For the z-projection instead, we only keep the $k=0$ term, since the ISBT is not substantially affected by non-locality. We then have $\tilde{u}_{\lambda=\text{TM},k} \approx \sqrt{f_{\text{isb}}\delta_{0k}}$, where $f_{\text{isb}} < 1$.

For the TE mode, instead, we assume that it is completely constant across the plane so to have

$$ \mathbf{w}_{\lambda=\text{TE}}(x, y, 0) \cdot \mathbf{u}_x \approx \sqrt{\frac{f}{L_x}}, \quad (B6) $$

and $\mathbf{w}_{\lambda=\text{TE}}(x, y, 0) \cdot \mathbf{u}_z \approx 0$. For the MP, its normalized Fourier transform gives

$$ \tilde{w}_{\lambda=\text{TE},k} = \delta_{k0} \sqrt{f}, \quad (B7) $$

and $\tilde{u}_{\lambda=\text{TE},k} = 0$ for the ISB.

We also assume that all the $N_{\text{qw}} = 57$ QWs couple in the same way to the cavity, and contribute equally to the Coulomb non-locality. We can thus include the factor $\sqrt{N_{\text{qw}}}$ as a collective enhancement of the plasma frequency. We have that the two modes MP susceptibility is now simplified to

$$ \begin{aligned} \chi_{\lambda=\text{TM}}^{\text{mp}}(\omega) &\approx f \frac{N_{\text{qw}} \omega_P^2}{2} \frac{\Delta x}{L_x} \sum_{k>0} \frac{\text{sinc}^2((k - k_1)\Delta x/2)}{\tilde{\omega}_B^2(k) - \omega^2 - i\kappa\omega} \\ \chi_{\lambda=\text{TE}}^{\text{mp}}(\omega) &\approx f \frac{N_{\text{qw}} \omega_P^2}{\omega_B^2 - \omega^2 - i\kappa\omega}, \end{aligned} \quad (B8) $$

while the ISBT gets

$$ \begin{aligned} \chi_{\lambda=\text{TM}}^{\text{isb}}(\omega) &\approx f_{\text{isb}} \frac{N_{\text{qw}} \omega_P^2}{2} \frac{1}{\omega_{\text{isb}}^2 - \omega^2 - i\kappa_{\text{isb}}\omega} \\ \chi_{\lambda=\text{TE}}^{\text{isb}}(\omega) &\approx 0. \end{aligned} \quad (B9) $$

It is worth noticing that in Eq. (B8) we neglected the cross contribution proportional to $\text{sinc}((k - k_1)\Delta x/2)\text{sinc}((k + k_1)\Delta x/2)$, which is a valid approximation if $\Delta x$ is sufficiently large, so that the overlap between these terms is negligible. Then, we can resum the $\pm k_1$ contributions, rewriting the formula limited to only the $k > 0$ modes.

In Eq. (A41) of $\tilde{\omega}_B^2(k)$, we must also use $N_{\text{qw}}\omega_P^2$ instead of $\omega_P^2$. This accounts for dipole-dipole near-field coupling between the slabs [47], for which the MP mode is actually a superradiant collective mode emerging from coupled motion of all the quantum wells together. This effective description holds well under the long-wavelength condition, for which $k\Delta z \ll 1$. Here $k$ is the typical MP wavelength and $\Delta z$ is the separation between the quantum wells along the $z$ axis. Since $k \sim 0 - 1/\Delta x$ we have that the inter-well separation must be $\Delta z \ll 10\mu$m, which is surely satisfied.

## Appendix C: Parameter estimation

Here we focus on the estimation for the parameters needed in the theory. Since the QW consists in a slab of GaAs, we take $\epsilon \approx 13\epsilon_0$ and $m \approx 0.067m_e$, where $\epsilon_0 \approx 8.85 \text{ pF/m} \approx 53.1 \text{ e/(}\mu\text{V m)}$ is the vacuum permittivity, $e$ and $m_e$ are the electron charge and mass. The three most important quantities in this work are

$$ \begin{aligned} \frac{e^2}{\epsilon_0} &\approx 1.8 \times 10^4 \text{ meV nm}, \\ \frac{\hbar^2}{m} &\approx 1136 \text{ meV nm}^2, \quad (C1) \\ \mu_B &= \frac{e\hbar^2}{2m} \approx 0.86 \text{ meV T}. \end{aligned} $$

Assuming a 2D electron density $n_{2D} \approx 10^{11} \text{cm}^{-2} = 10^{-3} \text{ nm}^2$ we have that $\hbar\omega_P \approx 0.4$ meV, corresponding to

$$ \frac{\omega_P}{2\pi} \approx 100 \text{ GHz}. \quad (C2) $$

(we remind the energy-to-frequency conversion is given by $2\pi\hbar \approx 4.13$ meV/THz).

Considering the collective enhancement due to the presence of many quantum wells, for instance $N_{\text{qw}} = 57$ and $n_{2D} \approx 3 \times 10^{-11}$, as reported in the main text, we have

$$ \frac{\sqrt{N_{\text{qw}}\omega_P}}{2\pi} \approx 1.3 \text{ THz}. \quad (C3) $$

This value is quite large and it would give a much bigger Rabi splitting than what is observed on the TE mode in the main text $\Omega_R/(2\pi) \approx 600$ GHz. However, accounting for the smaller overlap between the QW and the cavity field, $f < 1$ returns the correct value $\Omega_R/(2\pi) \approx 2f\sqrt{N_{\text{qw}}\omega_P/(2\pi)}$.

From the experimental data we observe the MP to cross the TE mode at around $B \approx 8.5$ T. Considering that $\omega_B = 2\mu_B B$, from this value we obtain

$$ \frac{\omega_B}{2\pi} \approx 3.5 \text{ THz}, \quad (C4) $$

indicating that the electron's mass is correctly $m* = 0.067m_e$.
From the independent observation of the MP resonance
linewidth, we obtain

$$
\frac{\kappa}{2\pi} \approx 300 \text{ GHz.} \qquad (C5)
$$

The two main cavity modes are given by $\lambda = \{\text{TM, TE}\}$, called transverse magnetic, and transverse electric, based on their polarization properties [44]. Both the square-well and parabolic well devices illustrated in the main text support these modes, due to their MIM, metallic-grating patch-cavity design.

<table><caption>TABLE I. Parameters used for the numerical simulations in the main text</caption>
  <thead>
    <tr>
      <th></th>
      <th>Square Well</th>
      <th>Parabolic Well</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>ω<sub>TM</sub>/(2π)</td>
      <td>2.408 THz</td>
      <td>1.5 THz</td>
    </tr>
    <tr>
      <td>ω<sub>TE</sub>/(2π)</td>
      <td>3.508 THz</td>
      <td>3.5 THz</td>
    </tr>
    <tr>
      <td>ω<sub>isb</sub>/(2π)</td>
      <td>2.73 THz</td>
      <td>1.4 THz</td>
    </tr>
    <tr>
      <td>√N<sub>qw</sub>ω<sub>P</sub>/(2π)</td>
      <td>1 THz</td>
      <td>1.4 THz</td>
    </tr>
    <tr>
      <td>f</td>
      <td>0.3</td>
      <td>0.3</td>
    </tr>
    <tr>
      <td>f<sub>isb</sub></td>
      <td>0.3</td>
      <td>0.4</td>
    </tr>
    <tr>
      <td>Δx</td>
      <td>10 μm</td>
      <td>6 μm</td>
    </tr>
    <tr>
      <td>L<sub>c</sub></td>
      <td>11.5 μm</td>
      <td>4 μm</td>
    </tr>
    <tr>
      <td>L<sub>x</sub></td>
      <td>30μm</td>
      <td>30 μm</td>
    </tr>
    <tr>
      <td>γ/(2π)</td>
      <td>50 GHz</td>
      <td>50 GHz</td>
    </tr>
    <tr>
      <td>κ/(2π)</td>
      <td>300 GHz</td>
      <td>300 GHz</td>
    </tr>
    <tr>
      <td>κ<sub>isb</sub>/(2π)</td>
      <td>200 GHz</td>
      <td>300 GHz</td>
    </tr>
  </tbody>
</table>

All the other parameters, such as cavity linewidth,
ISBT frequencies etc... are taken phenomenologically from
the experimental data, and are consistently benchmarked
with the standard literature. The parameters used in the
numerical simulations reported in the main text are sum-
marized in Table I. Notice that the cavity linewidth is
take equal for all the modes $\gamma_{\lambda} = \gamma$.

Appendix D: Comparison to simulation without the
longitudinal Coulomb contribution

In this last section, we present a simulation example
with the same parameters as described above and used
in the main text, but where we artificially remove the
shift induced by Coulomb effects, setting $\bar{\omega}_B(k) = \omega_B$.
The result is in Fig. 6. Without the inclusion of the
k-dependent shift due to the longitudinal Coulomb inter-
action, we can see much better the Rabi splitting between
the MP and the polariton lines, which is perfectly sym-
metric. This type of plot does not correctly match the
experimental result, as one can directly check by com-
parison with Fig. 3 and 4.

<figure><img src="image_7.png" /><figcaption>FIG. 6. (a-c) Same simulation as in Fig. 3. (b-d) Same simulation as in Fig. 4. In both cases $\zeta_k = 0$.</figcaption></figure>

[1] A. F. Kockum, A. Miranowicz, S. D. Liberato, S. Savasta, and F. Nori, Ultrastrong coupling between light and matter, Nature Reviews Physics 1, 19 (2019).

[2] P. Forn-Díaz, L. Lamata, E. Rico, J. Kono, and E. Solano, Ultrastrong coupling regimes of light-matter interaction, Rev. Mod. Phys. 91, 025005 (2019).

[3] G. Scalari, C. Maissen, D. Turčinková, D. Hagenmüller, S. D. Liberato, C. Ciuti, C. Reichl, D. Schuh, W. Wegscheider, M. Beck, and J. Faist, Ultrastrong coupling of the cyclotron transition of a 2D electron gas to a THz metamaterial, Science 335, 1323 (2012).

[4] C. Ciuti, G. Bastard, and I. Carusotto, Quantum vacuum properties of the intersubband cavity polariton field, Phys. Rev. B - Condensed Matter and Materials Physics 72, 1 (2005).

[5] D. D. Bernardis, A. Mercurio, and S. D. Liberato, Tutarial on nonperturbative cavity quantum electrodynamics: is the Jaynes Cummings model still relevant?, J. Opt. Soc. Am. B 41, C206 (2024).

[6] Y. Todorov, A. M. Andrews, R. Colombelli, S. D. Liberato, C. Ciuti, P. Klang, G. Strasser, and C. Sirtori,

Ultrastrong light-matter coupling regime with polariton dots, Phys. Rev. Lett. 105, 1 (2010).

[7] M. Geiser, F. Castellano, G. Scalari, M. Beck, L. Nevou, and J. Faist, Ultrastrong Coupling Regime and Plasmon Polaritons in Parabolic Semiconductor Quantum Wells, Phys. Rev. Lett. 108, 106402 (2012).

[8] D. Dietze, A. M. Andrews, P. Klang, G. Strasser, K. Unterainer, and J. Darmo, Ultrastrong coupling of inter-subband plasmons and terahertz metamaterials, Applied Physics Letters 103, 201106 (2013).

[9] B. Askenazi, A. Vasanelli, A. Delteil, Y. Todorov, L. C. Andreani, G. Beaudoin, I. Sagnes, and C. Sirtori, Ultrastrong light-matter coupling for designer reststrahlen band, New Journal of Physics 16, 043029 (2014).

[10] M. Jeannin, G. Mariotti Nesurini, S. Suffit, D. Gacemi, A. Vasanelli, L. Li, A. G. Davies, E. Linfield, C. Sirtori, and Y. Todorov, Ultrastrong Light-Matter Coupling in Deeply Subwavelength THz LC Resonators, ACS Photonics 6, 1207 (2019).

[11] P. Goulain, C. Deimert, M. Jeannin, S. Pirotta, W. J. Pasek, Z. Wasilewski, R. Colombelli, and J. Manceau,
THz Ultra-Strong Light-Matter Coupling up to 200K with Continuously-Graded Parabolic Quantum Wells, *Advanced Optical Materials*, **22020724** (2023).

[12] Y. Todorov, A. M. Andrews, I. Sagnes, R. Colombelli, P. Klang, G. Strasser, and C. Sirtori, Strong light-matter coupling in subwavelength metal-dielectric microcavities at terahertz frequencies, *Phys. Rev. Lett.* **102**, 1 (2009).

[13] F. Berkmann, T. Venanzi, L. Baldassarre, E. Campagna, E. T. Simola, L. D. Gaspare, C. Corley-Wiciak, G. Capellini, G. Nicotra, G. Sfuncia, A. Notargiacomo, E. Giovine, S. Cibella, M. Virgilio, G. Scalari, M. D. Seta, and M. Ortolani, Ultrastrong Coupling of $Si_{1-x}Ge_x$ Parabolic Quantum Wells to Terahertz Microcavities, *ACS Photonics* **11**, 2776 (2024).

[14] U. Iqbal, Y. Todorov, and C. Mora, Dynamical Coulomb blockade: An all electrical probe of the ultrastrong light-matter coupling regime, *Phys. Rev. Res.* **6**, 033097 (2024).

[15] S. De Liberato and C. Ciuti, Stimulated Scattering and Lasing of Intersubband Cavity Polaritons, *Phys. Rev. Lett.* **102**, 136403 (2009).

[16] A. Delteil, A. Vasanelli, P. Jouy, D. Barate, J. C. Moreno, R. Teissier, A. N. Baranov, and C. Sirtori, Optical phonon scattering of cavity polaritons in an electromagnetic device, *Phys. Rev. B* **83**, 081404 (2011).

[17] J.-M. Manceau, N.-L. Tran, G. Biasiol, T. Laurent, I. Sagnes, G. Beaudoin, S. De Liberato, I. Carusotto, and R. Colombelli, Resonant intersubband polariton-LO phonon scattering in an optically pumped polaritonic device, *Applied Physics Letters* **112**, 191106 (2018).

[18] M. Knorr, J. M. Manceau, J. Mornhinweg, J. Nespolo, G. Biasiol, N. L. Tran, M. Malerba, P. Goulain, X. Lafosse, M. Jeannin, M. Stefinger, I. Carusotto, C. Lange, R. Colombelli, and R. Huber, Intersubband polariton-polariton scattering in a dispersive microcavity, *Phys. Rev. Lett.* **128**, 10.1103/PhysRevLett.128.247401 (2022).

[19] A. A. Anappara, A. Tredicucci, G. Biasiol, and L. Sorba, Electrical control of polariton coupling in intersubband microcavities, *Applied Physics Letters* **87**, 051105 (2005).

[20] S. Pirotta, N.-L. Tran, A. Jollivet, G. Biasiol, P. Crozat, J.-M. Manceau, A. Bousseksou, and R. Colombelli, Fast amplitude modulation up to 1.5 GHz of mid-IR free-space beams at room-temperature, *Nature Communications* **12**, 799 (2021).

[21] H. Chung, I. Hwang, J. Yu, G. Boehm, M. A. Belkin, and J. Lee, Electrical Phase Modulation Based on Mid-Infrared Intersubband Polaritonic Metasurfaces, *Advanced Science* **10**, 2207520 (2023).

[22] M. Malerba, S. Pirotta, G. Aubin, L. Lucia, M. Jeannin, J.-M. Manceau, A. Bousseksou, Q. Lin, J.-F. Lampin, E. Peytavit, S. Barbieri, L. H. Li, A. G. Davies, E. H. Linfield, and R. Colombelli, Ultrafast (≈10 GHz) mid-IR modulator based on ultrafast electrical switching of the light-matter coupling, *Applied Physics Letters* **125**, 041101 (2024).

[23] M. Jeannin, J. M. Manceau, and R. Colombelli, Unified Description of Saturation and Bistability of Intersubband Transitions in the Weak and Strong Light-Matter Coupling Regimes, *Phys. Rev. Lett.* **127**, 187401 (2021).

[24] M. Jeannin, E. Cosentino, S. Pirotta, M. Malerba, G. Biasiol, J. M. Manceau, and R. Colombelli, Low intensity saturation of an ISB transition by a mid-IR quantum cascade laser, *Applied Physics Letters* **122** (2023).

[25] S. Rajabali, E. Cortese, M. Beck, S. D. Liberato, J. Faist, and G. Scalari, Polaritonic nonlocality in light-matter interaction, *Nature Photonics* **15**, 690 (2021).

[26] A. Bayer, M. Pozimski, S. Schambeck, D. Schuh, R. Huber, D. Bougeard, and C. Lange, Terahertz light-matter interaction beyond unity coupling strength, *Nano Letters* **17**, 6340 (2017).

[27] J. Mornhinweg, L. K. Diebel, M. Halbhuber, M. Prager, J. Riepl, T. Inzenhofer, D. Bougeard, R. Huber, and C. Lange, Mode-multiplexing deep-strong light-matter coupling, *Nature Communications* **15**, 1847 (2024).

[28] W. Kohn, Cyclotron resonance and de haas-van alphen oscillations of an interacting electron gas, *Phys. Rev.* **123**, 1242 (1961).

[29] J. J. Hopfield, Theory of the contribution of excitons to the complex dielectric constant of crystals, *Phys. Rev.* **112**, 1555 (1958).

[30] D. Hagenmüller, S. De Liberato, and C. Ciuti, Ultrastrong coupling between a cavity resonator and the cyclotron transition of a two-dimensional electron gas in the case of an integer filling factor, *Phys. Rev. B* **81**, 235303 (2010).

[31] F. Monticone, A. Mortensen, A. Fernández-Domínguez, Y. Luo, C. Tserkezis, J. Khurgin, T. Shahbazyan, A. Chaves, N. Peres, G. Wegner, K. Busch, H. HU, F. Della Sala, P. Zhang, C. Ciraci, J. Aizpurua, A. Babaze, A. G. Borisov, X. Chen, T. Christensen, W. Yan, Y. Yang, U. Hohenester, L. Huber, M. Wubs, S. De Liberato, P. A. Gonçalves, J. García de Abajo, O. Hess, I. Tarasenko, J. Cox, L. Jelver, E. Dias, M. Sanchez Sanchez, D. Margetis, G. Gomez-Santos, T. Stauber, S. Tretyakov, C. Simovski, S. Pakniyat, J. S. Gomez-Diaz, I. Bondarev, S.-A. Biehs, A. Boltasseva, V. Shalaev, A. Krasavin, A. Zayats, A. Alu, J.-H. Song, M. Brongersma, U. Levy, O. Long, C. Guo, S. Fan, S. Bozhevolnyi, A. Overvig, f. prudencio, M. Silveirinha, S. A. Hassani Gangaraj, C. Argyropoulos, P. Huidobro, E. Galiffi, F. Yang, J. Pendry, D. A. B. Miller, and X. Zheng, Roadmap on nonlocality in photonic materials and metamaterials, *Optical Materials Express* **15**, 1544 (2025).

[32] T. Maag, A. Bayer, S. Baierl, M. Hohenleutner, T. Korn, C. Schüller, D. Schuh, D. Bougeard, C. Lange, R. Huber, M. Mootz, J. E. Sipe, S. W. Koch, and M. Kira, Coherent cyclotron motion beyond Kohn's theorem, *Nature Physics* **12**, 119 (2016).

[33] J. Keller, G. Scalari, S. Cibella, C. Maissen, F. Appugliese, E. Giovine, R. Leoni, M. Beck, and J. Faist, Few-Electron Ultrastrong Light-Matter Coupling at 300 GHz with Nanogap Hybrid LC Microcavities, *Nano Letters* **17**, 7410 (2017).

[34] S. Rajabali, J. Enkner, E. Cortese, M. Beck, S. D. Liberato, J. Faist, and G. Scalari, Engineered planar plasmonic reflector for polaritonic mode confinement, *Opt. Mater. Express* **13**, 2944 (2023).

[35] S. R. Endo, D. Kim, S. Liang, G. Lee, S. Kim, A. Covarrubias-Morales, M. Seo, M. J. Manfra, D. Lee, M. Bamba, and J. Kono, Cavity-mediated coupling between local and nonlocal modes in landau polaritons (2025), arXiv:2509.05738 [quant-ph].

[36] F. J. Garcia-Vidal, C. Ciuti, and T. W. Ebbesen, Manipulating matter by strong coupling to vacuum fields, *Science* **373**, eabd0336 (2021).
[37] J. Enkner, L. Graziotto, D. Boriçi, F. Appugliese, C. Reichl, G. Scalari, N. Regnault, W. Wegscheider, C. Ciuti, and J. Faist, Tunable vacuum-field control of fractional and integer quantum hall phases, *Nature* **641**, 884 (2025).

[38] K. W. Chiu and J. J. Quinn, Plasma oscillations of a two-dimensional electron gas in a strong magnetic field, *Phys. Rev. B* **9**, 4724 (1974).

[39] Q. Zhang, T. Arikawa, E. Kato, J. L. Reno, W. Pan, J. D. Watson, M. J. Manfra, M. A. Zudov, M. Tokman, M. Erukhimova, A. Belyanin, and J. Kono, Superradiant decay of cyclotron resonance of two-dimensional electron gases, *Phys. Rev. Lett.* **113**, 047601 (2014).

[40] D. J. Griffiths, *Introduction to electrodynamics*, 4th ed. (Pearson, 2013).

[41] Y. Todorov and C. Sirtori, Intersubband polaritons in the electrical dipole gauge, *Phys. Rev. B* **85**, 045304 (2012).

[42] Y. Todorov, Dipolar quantum electrodynamics of the two-dimensional electron gas, *Phys. Rev. B* **91**, 125409 (2015).

[43] D. De Bernardis, T. Jaako, and P. Rabl, Cavity quantum electrodynamics in the nonperturbative regime, *Phys. Rev. A* **97**, 043820 (2018).

[44] J. D. Jackson, *Classical electrodynamics*, 3rd ed. (Wiley, New York, NY, 1999).

[45] C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg, *Photons and Atoms: Introduction to Quantum Electrodynamics* (Wiley, 1989).

[46] D. De Bernardis, M. Jeannin, J.-M. Manceau, R. Colombelli, A. Tredicucci, and I. Carusotto, Magnetic-field-induced cavity protection for intersubband polaritons, *Phys. Rev. B* **106**, 224206 (2022).

[47] D. De Bernardis, G. M. Andolina, and I. Carusotto, Light-matter interactions in the vacuum of ultrastrongly coupled systems, *Phys. Rev. A* **110**, 053713 (2024).

[48] V. Vladimirov, *Equations of Mathematical Physics*, Monographs and textbooks in pure and applied mathematics (M. Dekker, 1971).

[49] E. Lipparini, *Modern Many-Particle Physics*, 2nd ed. (World Scientific, 2008).

[50] J. J. Hopfield and D. G. Thomas, Theoretical and experimental effects of spatial dispersion on the optical properties of crystals, *Phys. Rev.* **132**, 563 (1963).

[51] K. Takae and A. Onuki, Applying electric field to charged and polar particles between metallic plates: Extension of the ewald method, *The Journal of Chemical Physics* **139**, 124108 (2013).