Introduction

Strong coupling of light and matter is governed by the exchange of energy via virtual excitations, at the Rabi frequency, \({\varOmega }_{{{{{{\rm{R}}}}}}}\). In the ultrastrong1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16 (\(0.1\le {\varOmega }_{{{{{{\rm{R}}}}}}}/{\omega }_{{{{{{\rm{c}}}}}}} < 1\)) and deep-strong7,9,14 coupling (DSC) regime (\({\varOmega }_{{{{{{\rm{R}}}}}}}/{\omega }_{{{{{{\rm{c}}}}}}} \ge 1\)), the vacuum ground state is characterised by an unusually large population of virtual excitations17,18 with non-classical properties, including significant single or two-mode squeezing17,18. This results in spectacular effects of cavity quantum electrodynamics (c-QED)19 including the vacuum Bloch-Siegert shift10, polaritonic chemistry20,21,22, the creation of photon-bound excitons23, vacuum-field induced charge transport13,15, strong nonlinearities24, tunnelling25 and coherent polariton scattering26. Landau polaritons4,5,7,27 of plasmonic THz resonators were the first optical systems to enter the DSC regime7,9,14 with \({\varOmega }_{{{{{{\rm{R}}}}}}}/{\omega }_{{{{{{\rm{c}}}}}}}=1.43\) and a large ground state population of 0.37 virtual photons7. Moreover, they enable non-adiabatic switching28 of the extremely squeezed quantum vacuum ground state, which may enable the observation of Unruh-Hawking radiation29,30 with a further boost of \({\varOmega }_{{{{{{\rm{R}}}}}}}/{\omega }_{{{{{{\rm{c}}}}}}}\). While previous investigations have demonstrated remarkable progress in this direction, they also discovered fundamental barriers given by light-matter decoupling31 or dissipation32. One of the most significant restrictions for increasing \({\varOmega }_{{{{{{\rm{R}}}}}}}/{\omega }_{{{{{{\rm{c}}}}}}}\), however, is the paradigm of resonant light-matter interaction of a single photonic and a single electronic mode, limiting the maximum contributing oscillator strength (Fig. 1a).

Fig. 1: Multi-mode light-matter coupling and ultracompact metasurface.
figure 1

a Illustration of resonant ultrastrong coupling of a single cavity mode (upper parabola) to a single matter excitation (bottom parabola) with a vacuum Rabi frequency \({\varOmega }_{{{{{{\rm{R}}}}}}}\). The weak population by virtual excitations in the vacuum ground state is indicated by semi-transparent spheres. b Illustration of deep-strong coupling of one light mode (upper parabola) with a frequency of \({\omega }_{j=1}\) to multiple matter excitations (bottom parabolas) with frequencies \({\omega }_{\alpha }\), by vacuum Rabi frequencies \({\varOmega }_{{{{{{\rm{R}}}}}},j,\alpha }\) under off-resonant conditions. Owing to the extremely large light-matter coupling, a significant number of virtual excitations are present. c Three-dimensional cut-away illustration of a conventional metasurface structure (gold shape) and its electric near-field distribution \(\widetilde{E}\), for the fundamental LC mode at ν1 = 0.52 THz. The unit cell (size: \(60\,{{{\rm{\mu m}}}}\times 60\,{{{\rm{\mu m}}}}\)) is indicated by the dashed line. QW: quantum well stack. d Highly compacted metasurface with a unit cell size of \(30\,{{{\rm{\mu m}}}} \times 32.5\,{{{\rm{\mu m}}}}\). e Measured THz transmission of the bare resonator array (blue circles), and calculation (red curve).

Full size image

Here, we present a conceptually different approach to DSC which overcomes the limitations of resonant coupling and establishes record coupling strengths: for sufficiently large spatial overlap of light and matter polarisation 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 strong 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 \({\varOmega }_{{{{{{\rm{R}}}}}}}^{{{{{{\rm{s}}}}}}}/{\omega }_{{{{{{\rm{s}}}}}}}=3.19\) for a hypothetical single pair of light and matter modes, for reference. Surpassing previous values almost twofold, this opens up unique possibilities for a range of non-adiabatic c-QED phenomena such as vacuum-field modified transport13,15, control of cavity chemistry20,21,22, and Unruh-Hawking radiation29,30. Moreover, our cavity mediates extremely strong coupling of plasmons that are orthogonal in their bare state - a distinct hallmark of the regime of deep-strong multi-mode coupling, which may in the future be exploited to drive phase transitions33 in quantum materials, solely by the vacuum field.

Results

Our resonator design capitalises on the principles of subwavelength near-field confinement34,35 of THz vacuum fields and current distributions in metasurfaces36, exploiting their great flexibility for designing the spectral shape, resonance frequencies, and near-field structure of optical modes. The structures are specifically optimised 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 layout37 (Fig. 1c). As a key advance, 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 nor strictly be classified as an inverse or a direct design. More precisely, its unit cell is approximately only a quarter as large as compared to conventional designs4,7, while a similar near-field enhancement is achieved. Its \(x\)-polarised modes include the fundamental LC resonance with a centre frequency of \({\nu }_{j=1}=0.52\,{{{{{\rm{THz}}}}}}\), while a higher-order mode lies at \({\nu }_{j=2}=1.95\,{{{{{\rm{THz}}}}}}\) (Fig. 1e), with j denoting the cavity mode index.

We experimentally evaluate our multi-mode concept by fabricating the gold resonators by electron-beam lithography on top of n-doped GaAs multi-quantum well (QW) heterostructures which host two-dimensional electron gases with a nominal charge carrier concentration of \({\rho }_{{{{{{\rm{QW}}}}}}}=1.8\times {10}^{12}{{{{{\rm{c}}}}}}{{{{{{\rm{m}}}}}}}^{-2}\), and are separated by AlGaAs barrier layers7 (see Methods). A magnetic bias field B of up to 5.5 T is oriented perpendicularly to the QW plane and Landau-quantises the electrons, achieving a cyclotron resonance (CR) with a tuneable frequency \({\nu }_{{{{{{\rm{c}}}}}}}={e|}{{{{{\bf{B}}}}}}|/2\pi {m}^{*}\), where e is the elementary charge, m* = 0.07 me is the electron effective mass and me is the free electron mass. Owing to the comparably large effective mass of electrons in gold of \({m}_{{{{{{\rm{Au}}}}}}}^{*}/{m}_{{{{{{\rm{e}}}}}}}\) ⪆ 1 (ref. 38), the CR of electrons in the resonator material itself is neglected.

In the assembled structure, the periodicity of the field localisation of the metasurface quantises 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 resonances32. The plasmons feature a frequency of \({\omega }_{{{{{{\rm{P}}}}}}}\left({{{{{{\bf{q}}}}}}}_{{{{{{\boldsymbol{x}}}}}}}\right)=\sqrt{\frac{{\rho }_{2{{{{{\rm{D}}}}}}}{e}^{2}}{2{m}^{*}{\epsilon }_{0}{\epsilon }_{{{{{{\rm{eff}}}}}}}\left({{{{{\bf{q}}}}}}\right)}\left|{{{{{{\bf{q}}}}}}}_{{{{{{\boldsymbol{x}}}}}}}\right|}\) (refs. 39,40) (Fig. 2a, red graph), where \({\rho }_{2{{{{{\rm{D}}}}}}}\) is the total carrier density integrated over all QWs, and \({\epsilon }_{0}\) and \({\epsilon }_{{{{{{\rm{eff}}}}}}}\left({{{{{\bf{q}}}}}}\right)\) are the vacuum permittivity and effective dielectric function, respectively. The plasmon wave vector \({{{{{{\bf{q}}}}}}}_{x}\left(\alpha \right)=\frac{2\pi }{{L}_{x}}\alpha\) depends on the unit cell size in \(x\)-direction, \({L}_{x}\) and the mode index, \(\alpha {\mathbb{\in }}{\mathbb{Z}}\) (see Supplementary Information for the full theory). The small value of \({L}_{x}=30\,{{{\rm{\mu 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 (\({{{{{{\boldsymbol{-}}}}}}{{{{{\bf{q}}}}}}}_{x},{{{{{{\bf{q}}}}}}}_{x}\)) form bright and dark standing-wave modes which hybridise with the CR at \({\omega }_{{{{{{\rm{c}}}}}}}=2\pi {\nu }_{{{{{{\rm{c}}}}}}}\), forming 2|αc| \(+1=21\) magnetoplasmon (MP) resonances (see Supplementary Information). Each MP with a frequency of \({\omega }_{{{{{{\rm{MP}}}}}},\alpha }\propto \sqrt{{\omega }_{{{{{{\rm{P}}}}}},\alpha }^{2}+{\omega }_{{{{{{\rm{c}}}}}}}^{2}}\) (ref. 41; Fig. 2b) couples to the resonator modes, \(j\), with individual vacuum Rabi frequencies \({\varOmega }_{{{{{{\rm{R}}}}}},j,\alpha }\) depending on the MP amplitude (Fig. 2a, blue bars) and the detuning, \({\nu }_{j}-{\omega }_{{{{{{\rm{MP}}}}}},\alpha }/2\pi\). For any given value of \({\nu }_{{{{{{\rm{c}}}}}}}\), the detuning vanishes only for a single MP mode at a time (Fig. 2b).

Fig. 2: Deep-strong light-matter coupling.
figure 2

a Plasmon frequency \({\nu }_{{{{{{\rm{p}}}}}},\alpha }\) and coupling strength \({\varOmega }_{{{{{{\rm{R}}}}}},1,\alpha }/{\omega }_{j=1}\) for each MP (magnetoplasmon) mode \(\alpha\) in the case of the sample with 48 QWs (quantum wells). b, c Illustration of multi-mode coupling of one cavity mode to several matter modes, MP0 = CR, MP1 and MP2, as a function of the cyclotron frequency, \({\nu }_{{{{{{\rm{c}}}}}}}\). b uncoupled modes. c Coupled modes comprising of one lower polariton (LP1) and three upper polaritons (UP1,1, UP1,2 and UP1,3). d THz magneto-transmission \({{{{{\mathcal{T}}}}}}\) as a function of \({\nu }_{{{{{{\rm{c}}}}}}}\) of the single-QW structure. The continuous white curves trace the polariton modes obtained from the multi-mode Hopfield model for the first resonator mode (coupling strength: \({\eta }_{1}\) = 0.55). The dotted white curves represent the polaritons linked to the higher mode ν2 (\({\eta }_{2}\) = 0.13). h Spectrum obtained from time-domain quantum model and identical polariton frequencies, for comparison. e Transmission of the 3-QW structure (\({\eta }_{1}\) = 0.76) and i simulation. f Transmission of the 6-QW structure (\({\eta }_{1}\) = 1.34), and j simulation. g Transmission of the 12-QW structure (\({\eta }_{1}\) = 2.32), and k simulation.

Full size image

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 orthogonal16, only a single cavity mode. We thus reuse the cavity mode index \(j\) and introduce the index \(\beta\) for the resulting 22 magnetoplasmon cavity polaritons (MPPs). These are categorised into a lower polariton, \({{{{{\rm{L}}}}}}{{{{{{\rm{P}}}}}}}_{j,(\beta=0)}\), and, as a signature of multi-mode coupling, several upper polaritons, \({{{{{\rm{U}}}}}}{{{{{{\rm{P}}}}}}}_{j,1\le \beta \le {\alpha }_{{{{{{\rm{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 (\(\beta < 0\)) are included in our model but not further discussed here (see Supplementary Information). As shown in the simplified example in Fig. 2c, all MPP frequencies rise monotonically when \({\nu }_{{{{{{\rm{c}}}}}}}\) is increased from zero. The \({{{{{\rm{L}}}}}}{{{{{{\rm{P}}}}}}}_{j=1}\) asymptotically emerges from below the CR (dotted line) near \({\nu }_{{{{{{\rm{c}}}}}}}\approx 0\) and converges towards the cavity frequency for large \({\nu }_{{{{{{\rm{c}}}}}}}\). The \({{{{{\rm{U}}}}}}{{{{{{\rm{P}}}}}}}_{j=1,\beta }\) modes, for \({\nu }_{{{{{{\rm{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 field17. All UP modes bend upwards with increasing \({\nu }_{{{{{{\rm{c}}}}}}}\) and asymptotically converge towards the CR.

In a first experimental campaign, we investigate this distinct 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 \({\nu }_{{{{{{\rm{c}}}}}}}\) (see Methods). For the single-QW structure (Fig. 2d), the single LP1 mode is located below the CR for vanishing \({\nu }_{{{{{{\rm{c}}}}}}}\) and approaches the frequency of the bare first cavity mode at 0.5 THz, as \({\nu }_{{{{{{\rm{c}}}}}}}\) is increased. For \({\nu }_{{{{{{\rm{c}}}}}}}=0\), the \({{{{{\rm{U}}}}}}{{{{{{\rm{P}}}}}}}_{1,\beta }\) are observed as a dense fan of partially overlapping modes distributed from \(\sim 0\,{{{{{\rm{THz}}}}}}\) upwards. The most strongly visible UP mode, highlighted by the uppermost semi-transparent white curve, is highest in frequency with \(\nu=0.78\,{{{{{\rm{THz}}}}}}\). Increasing \({\nu }_{{{{{{\rm{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 \(\nu=1.65\,{{{{{\rm{THz}}}}}}\) at \({\nu }_{{{{{{\rm{c}}}}}}}=1.9\,{{{{{\rm{THz}}}}}}\). The \({{{{{\rm{U}}}}}}{{{{{{\rm{P}}}}}}}_{2,\beta }\) ensemble (upper dotted curves) is dominated by a single spectral feature centred near \(\nu=2.0\,{{{{{\rm{THz}}}}}}\) for \({\nu }_{{{{{{\rm{c}}}}}}}=0\,{{{{{\rm{THz}}}}}}\) while slightly shifting upwards with increasing \({\nu }_{{{{{{\rm{c}}}}}}}\).

For \({n}_{{{{{{\rm{QW}}}}}}}=3\), the overall electronic dipole moment is boosted and the frequency spacing \({\omega }_{{{{{{\rm{P}}}}}}}\left({{{\alpha }}}\right)\) of the uncoupled MPs and, as a result, that of the coupled MPPs is increased (Fig. 2e). For \({\nu }_{{{{{{\rm{c}}}}}}}=0\) THz, the \({{{{{\rm{U}}}}}}{{{{{{\rm{P}}}}}}}_{1,\beta }\) extend from \(\sim 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 \({{{{{\rm{U}}}}}}{{{{{{\rm{P}}}}}}}_{2,\beta }\) modes further up while lowering the frequencies of the \({{{{{\rm{L}}}}}}{{{{{{\rm{P}}}}}}}_{j}\) modes. These trends continue for \({n}_{{{{{{\rm{QW}}}}}}}=6\) (Fig. 2f), where for \({\nu }_{{{{{{\rm{c}}}}}}}=0\) THz, the MPP fan manifests in well-separated \({{{{{\rm{U}}}}}}{{{{{{\rm{P}}}}}}}_{1,\beta }\) modes which extend up to 1.65 THz, whereas the \({{{{{\rm{U}}}}}}{{{{{{\rm{P}}}}}}}_{2,\beta }\) reach up to 2.5 THz. Yet, owing to the coupling of all MP to the same common cavity mode, only one LP forms for each cavity mode. Implementing \({n}_{{{{{{\rm{QW}}}}}}}=12\) (Fig. 2g), we observe distinct \({{{{{\rm{U}}}}}}{{{{{{\rm{P}}}}}}}_{1,\beta }\) and \({{{{{\rm{U}}}}}}{{{{{{\rm{P}}}}}}}_{2,\beta }\) branches at frequencies reaching 2.5 THz and 3.2 THz, several of which are simultaneously ultrastrongly coupled with centre frequencies reaching up to \(\approx\)5 times the bare cavity frequency, \({\nu }_{j=1}=0.52\,{{{{{\rm{THz}}}}}}\).

Such extremely strong, off-resonant coupling of multiple light and matter modes with highly disparate frequencies represents an advanced setting of c-QED in which strongly subcycle exchange of vacuum energy between the coupled modes is expected. For a quantitative understanding, we developed a multi-mode theory of DSC based on the following Hopfield Hamiltonian27,42:

$$\widehat{{{{{{\mathcal{H}}}}}}}= \mathop{\sum }\limits_{j}{{\hslash }}{\omega }_{j}{\widehat{a}}_{j}^{{{\dagger}} }{\widehat{a}}_{j}+\mathop{\sum }\limits_{\alpha }{{\hslash }}{\omega }_{{{{{{\rm{MP}}}}}}}\left(\alpha \right){\widehat{b}}_{\alpha }^{{{\dagger}} }{\widehat{b}}_{\alpha }+\mathop{\sum }\limits_{\alpha,j}{{{\hslash }}\varOmega }_{{{{{{\rm{R}}}}}},j,\alpha }\left({\widehat{a}}_{j}^{{{\dagger}} }+{\widehat{a}}_{j}\right)\left({\widehat{b}}_{\alpha }^{{{\dagger}} }+{\widehat{b}}_{\alpha }\right) \\ +\mathop{\sum }\limits_{\alpha,j}\frac{{{\hslash }}{\varOmega }_{{{{{{\rm{R}}}}}},j,\alpha }^{2}}{{\omega }_{{{{{{\rm{MP}}}}}}}(\alpha )}{\left({\widehat{a}}_{j}^{{{\dagger}} }+{\widehat{a}}_{j}\right)}^{2}+{\widehat{{{{{{\mathcal{H}}}}}}}}_{{{{{{\rm{ext}}}}}}}.$$
(1)

The first two terms describe the bare cavity and MP resonances with frequencies \({\omega }_{j}=2\pi {\nu }_{j}\) and \({\omega }_{{{{{{\rm{MP}}}}}}}\left(\alpha \right)\) and bosonic annihilation operators \({\widehat{a}}_{j}\) and \({\widehat{b}}_{\alpha }\), respectively. The third term describes their mutual coupling, quantified by a set of vacuum Rabi frequencies, \({\varOmega }_{{{{{{\rm{R}}}}}},j,\alpha }\), and the fourth term represents the diamagnetic cavity blue shift while \({\widehat{{{{{{\mathcal{H}}}}}}}}_{{{{{{\rm{ext}}}}}}}\) implements coupling to external fields.

The individual vacuum Rabi frequencies, \({\varOmega }_{{{{{{\rm{R}}}}}},j,\alpha }\), are chosen proportional to the near-field amplitude component at the corresponding MP wave vector obtained by Fourier transform of the cavity field, and are multiplied by a common factor \(\lambda={E}_{j}/\sqrt{{V}_{j}}\) which scales all \({\varOmega }_{{{{{{\rm{R}}}}}},j,\alpha }\) with the amplitude \({E}_{j}\) and the volume \({V}_{j}\) of each cavity mode5 and is adjusted for best agreement with the experiment. We include the decay of the near-field amplitude in growth direction, thereby accounting for the coupling to each QW individually (Supplementary Information). Solving Heisenberg’s equations of motion for our Hamiltonian (see Methods) unveils the temporal evolution and, by Fourier transform, the spectral shape of the coupled modes. The theory accurately reproduces the measured spectra (Fig. 2h–k) regarding the line widths and oscillator strengths of all modes over the complete experimentally accessible spectral range from <0.1 THz to 6 THz, with quantitative precision. With multiple matter modes off-resonantly coupled to multiple light modes with extreme strengths, the notion of the anti-crossing point that maximises the coupling strength of a single pair of modes becomes irrelevant. More specifically, when \({\varOmega }_{{{{{{\rm{R}}}}}},j,\alpha }/{\omega }_{{{{{{\rm{MP}}}}}}}\left(\alpha \right)\approx 1\) is simultaneously reached for several MP modes \(\alpha\), the field fluctuations of off-resonant modes influence the vacuum ground state |G〉 only slightly less than resonant ones (see Supplementary Information). In this setting, the number of virtual photons \(\langle {N}_{j}\rangle=\langle G|{\hat{a}}_{j}^{{{\dagger}} }{\hat{a}}_{j}|G\rangle\) of each photonic mode \(j\) 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 \({\sum }_{j}\left\langle {N}_{j}\right\rangle\) to be increased almost arbitrarily by adding electronic oscillator strength within a spectral window up to several octaves wide.

For our 1-QW structure, we calculate \(\left\langle {N}_{1}\right\rangle=0.07\) and \(\left\langle {N}_{2}\right\rangle=4\times {10}^{-3}\). For comparison with previous demonstrations, the hypothetical coupling strength of a single pair of resonant light and matter modes that yields an equivalent vacuum photon number, \({\eta }_{j}={\varOmega }_{{{{{{\rm{R}}}}}},{{{{{\rm{j}}}}}}}^{{{{{{\rm{s}}}}}}}/{{{{{{\rm{\omega }}}}}}}_{j}\), amounts to \({\eta }_{1}=0.55\) and \({\eta }_{2}=0.13\) for the first two cavity modes. The calculation for the 3-QW structure results in an increased number of virtual photons and effective coupling strength for the cavity modes (Tab. 1). The model also reproduces the transition from the merged ensemble of \({{{{{\rm{U}}}}}}{{{{{{\rm{P}}}}}}}_{1,\beta }\) modes (Fig. 2h) to clearly distinguishable individual MPP resonances (Fig. 2i–k). In addition, the calculation replicates the data for \({n}_{{{{{{\rm{QW}}}}}}}=6\) (Tab. 1). For \({n}_{{{{{{\rm{QW}}}}}}}=12\) we reach \(\left\langle {N}_{1}\right\rangle=0.76\), \(\left\langle {N}_{2}\right\rangle=0.08\) and \({\eta }_{1}=2.32\),\({\eta }_{2}=0.60\). While these values already exceed previous records significantly7,14, we push the limits of our approach even further with two additional structures employing 24 and 48 QWs, respectively. With the comparably large extension of these QW stacks in the \(z\) direction of up to ~2 μm, the lowermost QWs are located at a distance from the surface where the near-field amplitude of the bare cavity is significantly reduced relative to the plane of the resonators. However, owing to superradiant coupling prevalent at high densities of Landau electrons43, the polarisation components of all Landau electrons throughout the QW stack are synchronised in phase to a high degree. The optical mode thus couples to this collective polarisation field of all QWs in the entire stack, increasing the effective mode volume. While the coupling strength is hereby reduced, its increase by the boosted net oscillator strength by far dominates. Both effects are included in our FEFD model. The deep-strongly coupled polariton modes of these optimised structures extend over an even wider frequency range and create a highly structured spectrum (Fig. 3a,c). The increased interaction strength moreover causes a small signature of coupling of \(x\) and \(y\)-polarised components by the gyrotropic nature of the CR. As a result, we observe an additional, faint polariton branch located ~200 GHz above the LP2 mode for \({\nu }_{{{{{{\rm{c}}}}}}}=\) 2 THz, which can be attributed to a \(y\)-polarised resonator mode and is invisible for our \(x\)-polarised probe pulses at lower coupling strengths.

Fig. 3: Extremely strong, multi-octave light-matter coupling.
figure 3

a THz magneto-transmission \({{{{{\mathcal{T}}}}}}\) of the 24-QW sample as a function of the cyclotron frequency, \({\nu }_{{{{{{\rm{c}}}}}}}\) (see Fig. 2). The extended Hopfield model yields coupling strengths of \({\eta }_{1}=2.80\) and \({\eta }_{2}=0.85\) for the first and second resonator mode, respectively. Calculated polariton frequencies (solid & dotted curves) with distinct resonances labelled. b Calculated transmission and identical polariton frequencies, for comparison. c Transmission of the 48-QW structure. Coupling strengths: \({\eta }_{1}=2.83\), \({\eta }_{2}=0.88\). d Calculated transmission and identical polariton frequencies.

Full size image

Our calculations (Fig. 3b,d) confirm yet higher vacuum photon populations and coupling strengths for \({n}_{{{{{{\rm{QW}}}}}}}=24\) (Tab. 1), finally achieving \(\left\langle {N}_{1}\right\rangle=1.00\), \(\left\langle {N}_{2}\right\rangle=0.17\) as well as \({\eta }_{1}=2.83\),\({\eta }_{2}=0.88\) for \({n}_{{{{{{\rm{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, \(\left\langle N\right\rangle=\left\langle {N}_{1}\right\rangle+\left\langle {N}_{2}\right\rangle=1.17\), comfortably surpasses it. The latter value would be found for an effective coupling strength of \({\varOmega }_{{{{{{\rm{R}}}}}}}^{{{{{{\rm{s}}}}}}}/\omega_{{{\rm{s}}}}=3.19\) – exceeding existing structures with \(\eta=1.43\) (ref. 7) and \(\eta=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 transport13,15, chemical reactions20,21,22 or phase transitions33, 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 polarisation dynamics of our structure by the transmitted THz field of the 48-QW sample. At \({\nu }_{{{{{{\rm{c}}}}}}}=0.52\) THz, the experimental waveform consists of a pronounced initial cycle at a delay time of \(t=0\) ps (Fig. 4a). From \(t=0.5\) ps onwards, rapid trailing oscillations are observed which exhibit 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, grey-shaded 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, \({\varOmega }_{{{{{{\rm{R}}}}}},j=1,\alpha }\). 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 polarisation of the first MP mode (\({\nu }_{{{{{{\rm{MP}}}}}},\alpha=1}=1.2\) THz) is strongly structured by multiple oscillations (Fig. 4c). Its spectrum is characterised by four main contributions located between 0.18 and 2.59 THz, as well as further, less pronounced components at higher frequencies (Fig. 4c, inset). Moreover, since the large coupling strengths \({\varOmega }_{{{{{{\rm{R}}}}}},j=1,\alpha }\) to the same cavity mode cause all MP modes to strongly influence each other, the dominant frequencies in their polarisation 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 \({\varOmega }_{{{{{{\rm{R}}}}}}}^{{{{{{\rm{s}}}}}}}/{\omega }_{{{{{{\rm{s}}}}}}}=3.19\) is much less structured (Fig. 4d), and its spectrum contains only the lower and upper polariton resonances (Fig. 4d, inset).

Fig. 4: Dynamics and squeezing of extremely strong light-matter coupling.
figure 4

a Transmitted THz field of the 48-QW sample (\({\eta }_{1}\) = 2.83) at \({\nu }_{{{{{{\rm{c}}}}}}}=0.52\) THz (black curve). Inset: Spectral amplitude \({{{{{\mathcal{A}}}}}}\) of the THz field. b Calculated expectation value for the population of the first cavity mode \({{{{\mathrm{Re}}}}}\left\langle {\hat{a}}_{j=1}\right\rangle\) after excitation (black curve) by a broadband pulse (grey curve). The shading marks one oscillation period of the bare cavity mode. Inset: Corresponding spectra. c Calculated expectation value of the polarisation of the first MP mode, \({{{{\mathrm{Re}}}}}\left\langle {\hat{b}}_{\alpha=1}\right\rangle\), and spectrum (inset). d \({{{{\mathrm{Re}}}}}\left\langle {\hat{a}}_{j=1}\right\rangle\) for the same coupling strength as in b, yet only for a single pair of light and matter modes. Inset: Corresponding spectrum. e Energy of the first cavity mode for the full calculation (solid curve) and the single-mode reference (dotted curve) as in panel d. f Energies of the first MP mode (solid curve) and CR (dotted curve) for the two cases. g Corresponding coupling energies between cavity mode and MP mode (solid curve) or CR (dotted curve). h Wigner-function representation for the photonic state of a deep-strongly coupled system with \({\varOmega }_{{{{{{\rm{R}}}}}}}/{\omega }_{{{{{{\rm{c}}}}}}}\) = 3.0 showing the quasi-probability \({{{{{\mathcal{A}}}}}}\) as a function of the field quadratures \({{{{\mathrm{Re}}}}}\left(\alpha \right)\) and \({{{{{\rm{Im}}}}}}\left(\alpha \right)\).

Full size image

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 \({\varOmega }_{{{{{{\rm{R}}}}}}}^{{{{{{\rm{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 coupling28 (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 state17.

In conclusion, our work represents a leap forward in c-QED by overcoming the limitations of resonant light-matter coupling with our innovative concept of cooperative, multi-mode hybridisation, 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 \({\Omega }_{{{{{{\rm{R}}}}}}}^{{{{{{\rm{s}}}}}}}/{\omega }_{{{{{{\rm{s}}}}}}}=3.19\) which exceeds previous records of THz c-QED more than twofold. As our structures permit femtosecond optical switching28, our tenfold increase of the areal density of virtual excitations will highly benefit the detection of vacuum radiation30,44. The extremely strong coupling of multiple electronic modes to the common cavity mode even facilitates hybridisation 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 high flexibility and an extensive level of control which can be leveraged, for example, in electronic transport13,15, light-induced phase transitions33 or chemical reactions20,21,22, by the interaction with vacuum fields.

Methods

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/GaAs superlattice to obtain a defect-free, atomically flat surface. The GaAs quantum well (QW) stacks were embedded in Al0.3Ga0.7As barriers, creating rectangular potential wells. Si δ-doping layers were placed symmetrically around the individual QWs, enabling control of the carrier density \({\rho }_{{{{{{\rm{QW}}}}}}}\) of the two-dimensional electron gases (2DEGs) formed in each QW. The spacing of the QWs and the doping density have been optimised following the strategy discussed in ref. 7. In particular, we took care to exclude coupling between adjacent QWs. Throughout the manuscript, \({\rho }_{{{{{{\rm{QW}}}}}}}\) is the density per QW, while \({\rho }_{2{{{{{\rm{D}}}}}}}\) denotes the combined carrier density of the QW stack, integrated over all QWs along the \(z\)-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. Table 1 contains information on the QW and barrier thicknesses.

Table 1 Structural and coupling parameters of the 6 QW structures
Full size table

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, \({{{{{{\bf{|q}}}}}}}_{x}\left(\alpha \right)|=\frac{2\pi }{{L}_{x}}\alpha,\alpha {\mathbb{\in }}{\mathbb{Z}}\), where \({L}_{x}\) is the unit cell size of the structure in \(x\)-direction, and \(\alpha\) is the plasmon mode index. For each \({{{{{{\bf{q}}}}}}}_{x}\left(\alpha \right)\), the relevant electric field amplitude of the cavity mode is given by the Fourier component \([{{{{{\mathcal{F}}}}}}\left({{{{{{\mathcal{E}}}}}}}_{x}\right)]({{{{{{\bf{q}}}}}}}_{x},{{{{{{\bf{q}}}}}}}_{y}=0)\), with \({{{{{\mathcal{F}}}}}}\left({{{{{{\mathcal{E}}}}}}}_{x}\right)\) denoting the 2D Fourier transform of the electric near-field, \({{{{{{\mathcal{E}}}}}}}_{x}\), in \(x\)-direction. \({{{{{{\mathcal{E}}}}}}}_{x}\) is calculated by the finite-element method for the bare resonator on an undoped substrate. The amplitudes are determined separately for each QW. 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 single-cycle 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 magneto-cryostat 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 \({\eta }_{{{{{{\rm{j}}}}}}}\) for each cavity mode, and the corresponding vacuum photon numbers, \({N}_{{{{{{\rm{j}}}}}}}\), are listed in Table 1.