Mode-multiplexing deep-strong light-matter coupling

Dressing electronic quantum states with virtual photons creates exotic effects ranging from vacuum-field modified transport to polaritonic chemistry, and squeezing or entanglement of modes. The established paradigm of cavity quantum electrodynamics maximizes the light-matter coupling strength \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varOmega }_{{{{{{\rm{R}}}}}}}/{\omega }_{{{{{{\rm{c}}}}}}}$$\end{document}ΩR/ωc, defined as the ratio of the vacuum Rabi frequency and the frequency of light, by resonant interactions. Yet, the finite oscillator strength of a single electronic excitation sets a natural limit to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varOmega }_{{{{{{\rm{R}}}}}}}/{\omega }_{{{{{{\rm{c}}}}}}}$$\end{document}ΩR/ωc. Here, we enter a regime of record-strong light-matter interaction which exploits the cooperative dipole moments of multiple, highly non-resonant magnetoplasmon modes tailored by our metasurface. This creates an ultrabroadband spectrum of 20 polaritons spanning 6 optical octaves, calculated vacuum ground state populations exceeding 1 virtual excitation quantum, and coupling strengths equivalent to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varOmega }_{{{{{{\rm{R}}}}}}}/{\omega }_{{{{{{\rm{c}}}}}}}=3.19$$\end{document}ΩR/ωc=3.19. The extreme interaction drives strongly subcycle energy exchange between multiple bosonic vacuum modes akin to high-order nonlinearities, and entangles previously orthogonal electronic excitations solely via vacuum fluctuations.

Here, we present a conceptually novel approach to DSC which overcomes the limitations of resonant coupling and establishes new record coupling strengths: for sufficiently large spatial overlap of light and matter polarization fields, even electronic excitations that are strongly detuned from the optical mode can substantially boost the vacuum ground state population.We leverage this idea by multiplexing the interaction, exploiting multiple, non-resonant plasmon modes which simultaneously and cooperatively couple with extreme strength to several optical modes of a metallic metasurface (Fig. 1b).Our resonator custom-tailors the plasmon quantum states, while significantly reducing the cavity size as compared to previous approaches, and provides large near-field enhancement.The resulting strongly subcycle transfer of vacuum energy between optical and electronic modes leads to an ultrabroadband spectrum of more than 20 distinct, extremely strongly coupled resonances distributed over up to 6 optical octaves.
The vacuum ground state is populated by up to 1.17 virtual photons, which leads to significant squeezing and corresponds to a record coupling strength of Ω R s /ω s = 3.19 for a hypothetical single pair of light and matter modes, for reference.Surpassing previous values almost twofold, this opens up new possibilities for a range of non-adiabatic c-QED phenomena such as vacuum-field modified transport 13 , control of cavity chemistry [19][20][21] , and Unruh-Hawking radiation 28,29 .Moreover, our cavity mediates extremely strong coupling of plasmons that are orthogonal in their bare state -a unique hallmark of the novel regime of deep-strong multi-mode coupling, which may in the future be exploited to drive phase transitions 32 in quantum materials, solely by the vacuum field.
Our resonator design capitalizes on the principles of subwavelength near-field confinement 33 of THz vacuum fields and current distributions in metasurfaces 34 , exploiting their great flexibility for designing the spectral shape, resonance frequencies, and near-field structure of optical modes.The structures are specifically optimized for multi-mode coupling and address the two central challenges of our c-QED scenario: designing a broadband spectrum of multiple distinct electronic modes to boost the coupling strength as well as reducing the unit cell size multifold over existing designs.The metasurface is inspired by an inverted layout 35 (Fig. 1c).As a key novelty, we discard the thought pattern of individual, spatially isolated resonators which are generously separated to avoid nearest-neighbour couplings.Instead, we perform finite-element frequency-domain (FEFD) calculations of the current distribution within the metal layer to identify and remove areas of low current density, which are irrelevant to the performance of the resonator.In addition, we exploit the symmetries of the structure to merge adjacent current paths of opposite phase, allowing us to eliminate the affected elements entirely (see Supplementary Information).The resulting metasurface (Fig. 1d) is highly compact and, as a result of the changed topology, can no longer be separated into individual resonators.More precisely, its unit cell is approximately only a quarter as large as compared to conventional designs, 4,7 while a similar near-field enhancement is achieved.Its fundamental LC resonance has a centre frequency of  =1 = 0.52 THz, while the second, dipolar mode lies at  =2 = 1.95 THz (Fig. 1e), with j denoting the cavity mode index.
We experimentally evaluate our multi-mode concept by fabricating the gold resonators by electronbeam lithography on top of n-doped GaAs multi-quantum well (QW) heterostructures which host twodimensional electron gases with a nominal charge carrier concentration of  QW = 1.8 × 10 12 cm −2 , and are separated by AlGaAs barrier layers 7 (see Methods).A magnetic bias field B of up to 5.5 T is oriented perpendicularly to the QW plane and Landau-quantizes the electrons, achieving a cyclotron resonance (CR) with a tuneable frequency  c = /2 * , where e is the elementary charge and m * is the electron effective mass.
In the assembled structure, the periodicity of the field localization of the metasurface quantizes the wave vectors of the light field that can couple to the electrons in the QWs (Fig. 2a).As a result, the mode's vacuum field couples to a fan of distinct plasmon resonances 31 .The plasmons feature a frequency of |  | (Refs. 36,37) (Fig. 2a, red graph), where  2D is the total carrier density integrated over all QWs, and  0 and  eff () are the vacuum permittivity and effective dielectric function, respectively.The plasmon wave vector   () = The coupling of the multiple plasmon modes to the electromagnetic modes of our resonator results in a set of multiple light-matter coupled modes.Each of these modes is a characteristic superposition of all MP modes and, since the cavity modes are almost orthogonal 15 , only a single cavity mode.We thus reuse the cavity mode index  and introduce the index  for the resulting 22 magnetoplasmon cavity polaritons (MPPs).These are categorized into a lower polariton, LP ,(=0) , and, as a novelty of multimode coupling, several upper polaritons, UP ,1≤≤ c +1 .The presence of only one LP mode, but several UP modes is the fingerprint of cooperative coupling of all independent MP modes to one common cavity mode.Dark modes ( < 0) are included in our model but not further discussed here (see Supplemental Information).As shown in the simplified example in Fig. 2c, all MPP frequencies rise monotonically when  c is increased from zero.The LP =1 asymptotically emerges from below the CR (dashed line) near  c ≈ 0 and converges towards the cavity frequency for large  c .The UP =1, modes, for  c = 0, start with a finite frequency, which is defined by the plasma frequencies of the constituent MPs as well as the diamagnetic shift caused by the photon field 16 .All UP modes bend upwards with increasing  c and asymptotically converge towards the CR.
In a first experimental campaign, we investigate this new regime of multi-mode coupling for four samples containing 1, 3, 6, and 12 QWs, respectively, by THz time-domain magneto-spectroscopy at cryogenic temperatures as a function of  c (see Methods).For the single-QW structure (Fig. 2d), the single LP1 mode is located below the CR for vanishing  c and approaches the frequency of the bare first cavity mode at 0.5 THz, as  c is increased.For  c = 0, the UP 1, are observed as a dense fan of partially overlapping modes distributed from ~0 THz upwards.The most strongly visible UP mode, highlighted by the uppermost semi-transparent white curve, is highest in frequency with  = 0.78 THz.Increasing  c , the entire UP mode structure curves upwards in frequency as discussed for Fig. 2c and occupies the increasingly smaller spectral bandwidth between the CR and the highest-energy UP mode.Similarly, the LP2 mode related to the second cavity mode (lowermost dotted curve) branches off below the CR as c is increased and reaches  = 1.65 THz at  c = 1.9 THz.The UP 2, ensemble (upper dotted curves) is dominated by a single spectral feature centred near  = 2.0 THz for  c = 0 THz while slightly shifting upwards with increasing  c .
For  QW = 3, the overall electronic dipole moment is boosted and the frequency spacing ω P () of the uncoupled MPs and, as a result, that of the coupled MPPs is increased (Fig. 2e).For  c = 0 THz, the UP 1, extend from ~0 to 1.0 THz and form clearly separated resonances, evidencing the multi-mode character of the interaction.The increased coupling likewise shifts the frequencies of the UP 2, modes further up while lowering the frequencies of the LP  modes.These trends continue for  QW = 6 (Fig. 2f), where for  c = 0 THz, the MPP fan manifests in well-separated UP 1, modes which extend up to The While these values already exceed previous records significantly 7,14 , we push the limits of our approach even further with two additional structures.Employing up to 48 QWs, we occupy almost the entire available cavity mode volume and further boost the electronic oscillator strength.The resulting deepstrongly coupled polariton modes extend over an even wider frequency range and create a highly structured spectrum (Fig. 3a,c).
Our calculations (Fig. 3b,d) confirm yet higher vacuum photon populations and coupling strengths for  QW = 24 (Tab.T1), achieving 〈 1 〉 = 1.00, 〈 2 〉 = 0.17 as well as  1 = 2.83,  2 = 0.88 for  QW = 48.Here, for the first time, the vacuum photon population of a single coupled optical mode reaches unity, while the combined ground state population of both modes, 〈〉 = 〈 1 〉 + 〈 2 〉 = 1.17, comfortably surpasses it.The latter value would be found for an effective coupling strength of Ω R s /ω = 3.19 -exceeding existing structures with  = 1.43 (Ref. 7) and  = 1.83 (Ref. 14) by almost a factor of 2. The strong back-action of the cavity vacuum fields moreover leads to a combined population of the MPs by 1.06 virtual excitations, and results in strong mixing of the formerly orthogonal matter modes.This perspective holds out the prospect of wide-ranging control of transport 13 , chemical reactions [19][20][21] or phase transitions 32 , merely by vacuum fluctuations.
The unconventional nature of this extremely strong multi-mode mixing is especially evident when the subcycle dynamics are investigated directly in the time-domain.To this end, we analyse the polarization dynamics of our structure by the transmitted THz field of the 48-QW sample.At  c = 0.52 THz, the experimental waveform consists of a pronounced initial cycle at a delay time of  = 0 ps (Fig. 4a).From  = 0.5 ps onwards, rapid trailing oscillations are observed which exhibit structured by multiple oscillations (Fig. 4c).Its spectrum is characterized by four main contributions from 0.18 to 2.59 THz and further, less pronounced components at higher frequencies (Fig. 4c, inset).
Moreover, since the large coupling strengths Ω R,=1, to the same cavity mode cause all MP modes to strongly influence each other, the dominant frequencies in their polarization reflect the spectral signatures of all MPPs simultaneously -a unique hallmark of multi-mode non-resonant DSC and the strong mixing of all MPs at the same time (see Supplementary Information).In comparison, the dynamics of the cavity field of a single pair of light-matter coupled modes with Ω R s /ω s = 3.19 is much less structured (Fig. 4d), and its spectrum contains only the lower and upper polariton resonances (Fig. 4d, inset).
An even stronger contrast between the settings of single and multi-mode coupling is evident from the dynamics of the energy redistribution: Whereas energy exchange of a single pair of modes progresses periodically at Ω R s (Fig. 4e-g, dotted curves), here, the large number of participating modes leads to irregularly structured dynamics for the cavity energy, each of the MP modes, and the energy stored in the coupling 27 (Fig. 4e-g, solid curves).Moreover, the extreme nature of the coupling manifests in its vacuum ground state, where strong squeezing of the photonic mode is observed (Fig. 4h).Further characteristics range from a highly non-classical Fock-state probability distribution (see Supplementary Information) to the emission of correlated photon pairs, expected for non-adiabatic modulation of this exotic vacuum ground state 16 .
In conclusion, our work represents a leap forward in c-QED by overcoming the limitations of resonant light-matter coupling with our novel concept of cooperative, multi-mode hybridization, allowing for a significant boost of light-matter coupling strengths.To this end, we designed a maximally compact resonator metasurface that custom-tailors multiple plasmon resonances as well as optical modes to form an ultrabroadband spectrum of Landau cavity polaritons covering 6 optical octaves.This extremely strong coupling results in a highly subcycle vacuum energy exchange and a vacuum ground state hosting a record population of 1.17 virtual photons and 1.06 virtual magnetoplasmon excitations.These figures correspond to a coupling strength of Ω R s /ω s = 3.19 which exceeds previous records of THz c-QED more than twofold.As our structures permit femtosecond optical switching 27 , our tenfold increase of the areal density of virtual excitations will highly benefit the detection of vacuum radiation 29,40 .The extremely strong coupling of multiple electronic modes to the common cavity mode even facilitates hybridization of otherwise orthogonal matter states and can be applied to interactions between a variety of systems such as magnons, phonons or Dirac electrons, including the mixing of entirely different excitations.Combined with the resulting sizeable virtual population of matter and light modes, our concept offers new flexibility and an unprecedented level of control which can be leveraged, for example, in electronic transport 13 , light-induced phase transitions 32 or chemical reactions [19][20][21] , by the interaction with vacuum fields.

Semiconductor heterostructure growth and electron-beam lithography
Our semiconductor heterostructures were grown by molecular-beam epitaxy on undoped (100)-oriented GaAs substrates which were prepared by growing an epitaxial GaAs layer of a thickness of 50 nm followed by an Al0.3Ga0.7As/GaAssuperlattice to obtain a defect-free, atomically flat surface.The GaAs quantum well (QW) stacks were embedded in Al0.3Ga0.7Asbarriers.Si -doping layers were placed symmetrically around the individual QWs, enabling control of the carrier density  QW of the twodimensional electron gases (2DEGs).Throughout the manuscript,  QW is the density per QW, while  2D denotes the combined carrier density of the QW stack, integrated over all QWs along the -direction.
The QW stacks were capped by a GaAs layer of a thickness of 30 nm (20 nm for the 12-QW sample) for protection against oxidation.
The THz resonators were fabricated on the surface of the semiconductor structures by electron-beam lithography and wet-chemical processing of polymethylmethacrylat (PMMA) resist followed by thermal vapour-phase deposition of 5 nm of Ti, improving adhesion of the subsequently deposited Au layer of a thickness of 100 nm.

Wave vector decomposition of resonator near field, and magnetoplasmon modes
The periodicity of our metasurfaces leads to a discretization of the plasmon wave vectors, ,  ∈ ℤ, where   is the unit cell size of the structure in -direction, and  is the plasmon mode index.For each   (), the relevant electric field amplitude of the cavity mode is given by the Fourier component [ℱ(ℰ  )](  ,   = 0), with ℱ(ℰ  ) denoting the 2D Fourier transform of the electric near-field, ℰ  , in -direction.ℰ  is calculated by the finite-element method for the bare resonator on an undoped substrate.The amplitudes are determined separately for each QW separately.Fig. S3 shows these data for the first two resonator modes, as a function of the wave vector and the depth below the metasurface.For large wave vectors, single-particle excitations become energetically possible and the field amplitude moreover drops sharply, limiting the relevant wave vector spectrum up to a maximum index | c | = 10.

Cryogenic THz magneto-spectroscopy
Femtosecond near-infrared pulses (centre wavelength, 807 nm; pulse energy, 5.5 mJ, pulse duration, 33 fs) from a titanium-sapphire amplifier laser (repetition rate, 3 kHz) were used to generate singlecycle THz pulses by optical rectification and to detect the transmitted waveforms by electro-optic sampling.Depending on the required bandwidth, we employed 110-cut ZnTe crystals of a thickness of 1 mm for our structures with up to 12 QWs, and GaP crystals of a thickness of 200 µm for the 24 and 48-QW structures.A mechanical chopper modulated the THz pulses, allowing for differential detection of the transmitted THz electric field, E(t).The sample was kept at a temperature of 10 K in a magnetocryostat with a large numerical aperture, enabling magnetic fields up to 5.5 T applied perpendicularly to the sample surface.The recorded electric field was Fourier transformed and referenced to a measurement without a sample in the cryostat to obtain the transmission spectra.

Sample structure parameters
For each sample, the doping densities, coupling strengths  j for each cavity mode, and the corresponding vacuum photon numbers,  j , are listed in table T1.The design of our gold metasurfaces aims at a maximally compact resonator geometry while providing spectrally well-separated optical modes with low linewidths, high near-field enhancement factors as well as low mode volumes.To this end, we have revisited some of the established design principles for plasmonic resonator structures.As detailed in the following, our new approach enables excess spacing between adjacent resonators to be minimized or completely discarded, increasing the areal resonator density and thus the total oscillator strength of the metasurface by a factor of ~4 as compared to the original design.
The inverted resonator structures are based on a square of 30 µm outer extension and twin current paths feeding a central capacitive gap of a width of 2.5 µm.The unit cell is twice as large as the outer dimensions of the resonator.In the frequency window of interest, the structure supports five optical modes with centre frequencies of 0.52, 1.95, 3.75, 4.6 and 6 THz (see main text, Table T1 and Fig. S2).
The fundamental mode at 0.52 THz is excited by THz radiation linearly polarized in x direction, leading to an oscillating current flow along the paths indicated in Fig. S1a, whereby the current lags behind the driving field by a phase of /2.As the field sweeps charges out of the right inner metal plate (turquoise right arrow), symmetry dictates that a current identical in direction and magnitude drives charges into the corresponding mirror plate on the left side of the structure (left turquoise arrow).Following Kirchhoff's laws, the loop is closed by currents passing through the outer metal plane symmetrically, in  (yellow arrows) and -direction (blue arrows).As a first step, we exploit the anti-symmetry of the oriented currents along the left and right edges of the outer metal plane.As the spacing of adjacent resonators in x-direction,   , is reduced, their near-field currents increasingly overlap.Owing to the opposite orientation of the currents marked by yellow arrows, they locally cancel each other.At the same time, the current flow into and out of the internal resonator area (turquoise arrows) remains unaffected since the role of the yellow currents is increasingly taken over by the current flows out of and into the internal resonator area of the nearest resonator neighbour.As a result, all relevant properties of the optical mode remain unaffected.Since this argument holds for any value of   , we chose   = 0, for which the metal area separating the resonators in x-direction vanishes completely.Here too, the currents exiting the inner metal plates of one resonator seamlessly continue to flow, now feeding the opposite metal plate of the adjacent resonator (Fig. S1b, turquoise arrows) instead of the now absent polarized current paths (Fig. S1a, yellow arrows).Likewise, the -polarized currents flowing along the outer perimeter of the structure are now connected between adjacent resonators (blue arrows) which also renders the -polarized paths unnecessary.
We next consider a reduction of the spacing of the structures in -direction,   .Since the outer currents of adjacent unit cells in -direction carry the same phase, the cancellation exploited above cannot be applied here.However, since the resonance condition for the LC mode favours the shortest overall current paths, the currents are generally concentrated very close to the edges of the metallized layer (Fig. S1), from where they fall off within  ̃ ≈ 3 µm in the direction normal to the edge, where  ̃ is the characteristic decay length, in -direction.As a consequence,   may be significantly reduced as long as the remaining metal stripline is of sufficient width to allow the -polarized current to continue to flow.In the final structure, we chose   = 2.5 µm ≈  ̃, leading to a reduction of the area of the unit cell by a factor close to 4. While the near-field distribution of this structure is virtually identical to that of the original one, the current distribution is markedly different.As Fig. S1 shows, the structure now exhibits two mainly -polarized subcurrents which sweep charges between the adjacent inner metal plates and along the striplines, respectively, and are offset by a phase of .Finally, since the arguments of symmetry and current localization apply to all resonator modes equally, the spectrum of the compacted structure does not differ significantly from the spectrum of the separated structures across the entire spectral range, apart from coupling to surface plasmons [1] and the overall increase in transmission due to the unit cell reduction (Fig. S2).

Magnetoplasmon modes in two-dimensional electron gases
In the two-dimensional electron gas (2DEG) hosted in our quantum wells (QWs), in absence of a magnetic field, plasma excitations obey the dispersion relation where  denotes the in-plane wavevector,  2D is the two-dimensional (2D) electron density,  * is the effective electron mass and  eff () represents the effective dielectric constant for the electron gas [2,3].
Applying a static magnetic bias field oriented perpendicularly to the plane of the 2DEG introduces Landau quantization and the plasmon excitations hybridize with the cyclotron resonance,  c , resulting in the formation of magnetoplasmons (MPs) [4] with a frequency of The effective dielectric constant  eff () for the MPs is set by the dielectric environment of the electron gas in the layers enclosing the 2DEG.In our structures, we consider the dielectric constant of the GaAs substrate below the QWs,  sub = 12.9, the capping layer of a thickness of  and dielectric constant  barrier on top of the QWs, as well as the dielectric function of the top interface to either vacuum, or the gold parts of the resonator structure.The latter two possibilities have been analysed previously [5,3], resulting in corresponding effective dielectric functions,  ungated () for the vacuum interface, and  gated () for a top-metallized 2DEG: Depending on the situation, either of the two functions assumes the role of  eff () in Eq. 1.For more complex, laterally structured samples such as planar metal resonators or gratings, an averaged, effective dielectric function can be constructed, proportionally factoring in the two situations according to a factor  describing the relative metal coverage of the surface: eff,mix () =  gated () + (1 − ) ungated (). (5) Moreover, the multi-QW stack consists of several layers with varying dielectric properties, requiring additional averaging of the effective dielectric function along the growth direction.As was shown previously, the response of the densely packed QWs is well approximated by an effective-medium approach [6].In addition, as the QW stack exhibits a finite extension in growth direction, the charge carriers in the QWs are not an ideal 2D electron gas.More precisely, with increasing QW thickness, their plasma frequency approaches the 3D plasma frequency as an asymptotical upper limit for large wave vectors, whereas the dispersion of an ideal 2D plasma has no such upper bound.To account for this effect, we must include a correction which depends on the QW stack thickness .The effective dielectric function  eff then is given by [7]: This description allows for a MP mode dispersion approaching the 3D plasma frequency asymptotically for || → ∞.However, plasmon and magnetoplasmon excitations are limited in frequency and wave vector by Landau damping which becomes effective when single-particle excitations become relevant [8,4].An upper bound for the frequency of these single-particle excitations is given by  <  F   , with the Fermi velocity  F = ℏ√2 2D  * .
The periodicity of our metasurfaces implies a discretization of the plasmon wave vectors that can be excited by the cavity modes.The corresponding condition for the in-plane wave vectors   is Here,   denotes the unit cell size of the structure in -direction, and  ∈ ℤ is the plasmon mode index.
Linear combinations of plasmon waves with wave vectors −  and   form bright and dark standing waves Ψ b ∝ exp(−  ) + exp(  ) and Ψ d ∝ exp(−  ) − exp(  ), respectively.Since only the bright modes couple to the cavity modes, the dark modes are not further considered.We verify this approach for our complex resonator geometry by calculating the 2D Fourier transform of the -polarized electric near-field component ℰ  of the first two cavity modes (Fig. S3a,b) along the -direction.The field ℰ  was calculated by the finite-element (FEM) method (see chapter 4) [6] and evaluated within -oriented planes that lie within the QW stack, whereby the doping concentration of the QWs was set to zero to obtain the modes of the empty cavity.This choice represents the field component most relevant for light-matter coupling, whereas we neglect the ℰ  -polarized components since they are significantly weaker in most areas of the structure.The resulting amplitudes, integrated within the respective plane, are shown in Fig. S3c,d as a function of the wave vector and the depth below the metasurface.Owing to the spectrum of single-particle excitations (see previous paragraph) and the fact that the field amplitude decays for larger wave vectors, in particular within planes more distant to the metasurface, we limit the magnetoplasmon mode index  to a maximum of || ≤   = 10.The data allows us to model the coupling of the light field to each QW separately and to account for the relative amplitude of each magnetoplasmon resonance.The same analysis is also performed for the higher resonator mode at 1.95 THz.The resulting plasmon dispersions, the wavevectors fulfilling the discretisation and the cutoff wavevector is shown in Fig. S4, whereas the relative coupling strengths for the plasmon modes in shown in the main manuscript Fig. 3a.

Parameter-free FEM model
In addition to our time-domain quantum model, we calculate the optical response of the bare resonator structure by numerical FEFD simulations following the concept of Ref. [6].The calculation yields the response of the coupled structures including the spatially resolved near-field distribution as well as the far-field transmission, without any free fit parameters.Our formalism requires the three-dimensional geometry of the nanostructure including the generally anisotropic dielectric response, implemented as a tensor function ( ⃗,   , ) which depends on the position  ⃗, the cyclotron frequency,  c , and the frequency of the light field, .Within the two-dimensional electron gas,  describes the gyrotropic nature of the cyclotron resonance in  and  direction owing to the -polarized static magnetic bias field.
In the  direction, we employ the background dielectric constant, since the small thickness of the quantum wells leads to a plasma frequency far above the frequencies of interest for our structure.
Additionally, we reduce the numerical complexity by modelling the quantum well stack including its barriers by an effective-medium approach with an effective dielectric tensor.The magnetically invariant

Scaling of the coupling strength
Following Hagenmüller et al. [12], the coupling strength of an ensemble of electronic oscillators coupled to a cavity mode scales with the square root of the total carrier density in the QW stack, Ω R ∝  0.5 .This dependency is perfectly reproduced for our first four structures containing 1, 3, 6 or 12 QWs.However, owing to limited penetration depth of the near-field into the QW stacks, we observe a deviation from this ideal scaling law for the structures with 24 and 48 QWs (Fig. S13).Nevertheless, the coupling strength still increases significantly even including this effect.
2   depends on the unit cell size in -direction,   and the mode index,  ∈ ℤ (see Supplementary Information for the full theory).The small value of   = 30 µm favourably increases the energy spacing as compared to previously investigated structures with larger unit cells7  , leading to more distinct resonances.Plasmon excitation is accounted for up to a maximum index | c | = 10, beyond which the near-field amplitude is negligible.Linear combinations of plasmon pairs (−  ,   ) form bright and dark standing-wave modes which hybridize with the CR at  c = 2 c , forming 2| c | + 1 = 21 magnetoplasmon (MP) resonances (see Supplemental Information).Each MP with a frequency of  MP, ∝ √ P, 2 +  c 2 (Ref. 38) (Fig. 2b) couples to the resonator modes, , with individual vacuum Rabi frequencies Ω R,j,α depending on the MP amplitude (Fig. 2a, blue bars) and the detuning,   −  MP, /2.For any given value of  c , the detuning vanishes only for a single MP mode at a time (Fig. 2b).
photonic mode  is the appropriate figure of merit of DSC.Moreover, the relaxed resonance criterion and the coupling to common cavity modes allows for the total vacuum photon population ∑ 〈  〉 beating patterns (indicated by arrows) resulting from the rapid energy exchange between multiple light and matter modes.The spectrum has a global maximum located near 2 THz and several small, adjacent local maxima resulting from individual MPP modes (Fig.4a, inset).Calculating the internal dynamics by our time-domain theory we find that the electric field of the first cavity mode (Fig.4b, black curve) oscillates on time scales much faster than the cycle duration of the bare mode, (1) -1 , (Fig.4b, greyshaded area).Its spectrum features 8 distinct local maxima of comparable amplitude located between 0.2 THz, where the LP1 is situated, and up to 4.5 THz (Fig.4b, inset).Exceeding the frequency 1 of the uncoupled cavity by almost an order of magnitude, these features directly result from subcycle energy transfer driven by extreme light-matter coupling strengths, Ω R,=1, .Further non-vanishing contributions extend to as low as 0.05 THz and up to 6 THz for a total spectral bandwidth of more than 6 optical octaves.Similarly, the polarization of the first MP mode ( MP,=1 = 1.2 THz) is strongly

Figure 1 |
Figure 1 | Multi-mode light-matter coupling and ultracompact metasurface.a, Illustration of

Figure 2 |Figure 3 |
Figure 2 | Deep-strong light-matter coupling.a, Plasmon frequency  p, and coupling strength Ω R,1, / =1 for each MP mode n in the case of the 48 QW sample.b,c, Illustration of multi-mode

Figure 4 |
Figure 4 | Dynamics and squeezing of extremely strong light-matter coupling.a, Transmitted THz field of the 48-QW sample ( 1 = 2.83) at  c = 0.52 THz (black curve).Inset: Spectral amplitude of the THz field.b, Calculated expectation value for the population of the first cavity mode ⟨ ̂=1 ⟩ after

Fig. S1 |
Fig. S1 | Ultracompact resonator design and current distributions.a, Resonator layout featuring a conventional, large unit cell.Arrows: Current flow of the fundamental optical mode.b, Resonator design compacted in -direction, and current flow.The unit cell is indicated by the dashed rectangle.c, Resonator design, compacted in  and -direction, and current flow.

Figure S2 |
Figure S2 | Transmission for the bare resonators.Transmission spectrum for excitation in x-direction of the resonator array featuring conventional, large unit cell (60 µm x 60 µm, red curve) and of the compacted structure (black curve).

Figure S3 |
Figure S3 | Analysis of magneto-plasmon formation.Electric field component ℰ  in the quantum well plane 50 nm below the resonator structure, for a, the LC mode (j = 1) at 0.52 THz and b, the dipolar mode (j = 2) at 1.95 THz.c, Amplitude components of the Fourier-transformed data of panel a as a function of the wave vector   , for   = 0, and the depth of the plane below the metasurface, .d, Equivalent amplitude components for the data of the dipolar mode in panel b.

Figure S4 |
Figure S4 | Plasmon dispersion for all samples.The dots mark the plasmon frequencies which fulfil the diffraction condition.The black dashed line indicates the cut-off wavevector up to which plasmon modes were considered in the quantum model.

Figure S8 |
Figure S8 | Far-field transmission of the bare metasurface calculated by the FEFD method (solid curve) and transmission obtained from the harmonic oscillator representation of our time domain simulations (dotted curve).

Figure S9 |
Figure S9 | Characteristics of THz excitation used in the time-domain theory.a, Waveform and b, amplitude spectrum.

Figure S10 |
Figure S10 | Switch-off analysis for the optical modes.The time-domain simulation of the 48-QW structure includes a, all modes, b, only the first and c, only the second optical mode of the resonator structure.

Figure S11 |
Figure S11 | Calculated spectra of the expectation values of the polarization of the MP modes 0 ≤  ≤ 10.The spectra are vertically offset by a value of 0.2, for clarity.
background and substrate layers are implemented by the dielectric function  GaAs (2 + TO  , including the optical phonons of GaAs to improve the accuracy of the calculation especially at higher THz frequencies.Here we use  ∞ = 10.87[10], LO /2 = 8.839 THz,  TO /2 = 8.124 THz,  LO = 0.0225 THz and  TO = 0.0255 THz[11].A single unit cell of the metasurface is implemented and periodically extended in  and  direction by periodic boundary conditions.The resulting far-field calculations predict experimental results across the entire spectral range with high accuracy (Fig.S12).

Figure S12 |
Figure S12 | Comparison of experimental spectra and FEFD calculations.a, Experimental THz magneto-transmission as a function of c for the single-QW, b, 3-QW, c, 6-QW, and d, 12-QW structure.e-h Corresponding FEFD THz transmission data.

S13|
Scaling of the coupling strength with the number of electronic oscillators.a, Equivalent coupling strength Ω R,1 / 1 for the first cavity mode as a function of the square root of the total charge carrier density, and linear fit of the data adjusted to the first four data points.b, Equivalent data for the second cavity mode.
1.65 THz, whereas the UP 2, reach up to 2.0 THz.Yet, owing to the coupling of all MP to the same common cavity mode, only one LP forms for each cavity mode.Implementing  QW = 12 (Fig.2g), we observe distinct UP 1, and UP 2, branches at frequencies reaching 2.5 THz and 3.2 THz, several of The authors thank Dieter Schuh and Imke Gronwald for valuable discussions and technical support.We J.M., M.H, D.B., R.H. and C.L. conceived the study.M.P., M.H. and D.B. designed, realised, and characterised the semiconductor heterostructures.J.M., M.H. and L.D. modelled the metasurfaces and fabricated the samples with support from T.I., J.R., D.B. and C.L. J.M., M.H., L.D. and J.R. carried out the experiments with support from R.H. and C.L. The theoretical modelling was carried out by J.M., R.H and C.L. D.B., R.H. and C.L. supervised the study.All authors analysed the data and discussed the results.J.M., R.H. and C.L. wrote the manuscript with contributions from all authors.

Table T1 | Structural and coupling parameters of the 6 QW structures
. 〈 1 〉, 〈 2 〉, 〈〉: vacuum photon numbers for the first and second resonator modes, and for both combined, respectively. 1 and  2 denote the hypothetical coupling strengths of, in each case, a single pair of light and matter modes with a vacuum photon population of 〈 1 〉 or 〈 2 〉, respectively, for comparison.

Table T1 |
Parameters for the cavity modes used to model the frequency response of the resonator structure.