Multimode ultrastrong coupling in three-dimensional photonic-crystal cavities

Polarization-dependent multimode coupling in a 3D-PCC

3D photonic crystals are dielectric structures that exhibit periodic refractive index modulations in all three dimensions. Here, we consider the woodpile structure where the rod arrays are stacked in an alternating orthogonal pattern, which is known to exhibit a complete photonic band gap^20,21. We used photolithography and deep reactive ion etching to fabricate a rod array in silicon wafers (Fig. 1a, Supplementary Fig. 9, see Methods for more details). A woodpile lattice^15,22 is formed by stacking patterned silicon wafers sequentially (Fig. 1b). To induce photonic modes in the photonic band gap, a planar defect (a 60-μm-thick GaAs layer) was inserted at the center of the woodpile lattice, which breaks discrete translational invariance in the stacking direction z (Fig. 1b). We performed numerical simulations with COMSOL Multiphysics to calculate the transmittance spectrum of the woodpile structure under normal incidence, as shown in Fig. 1d. To maintain high peak amplitude, only four silicon layers were placed on each side of the GaAs defect layer. Four cavity modes can be observed within the photonic band gap. Those cavity modes exhibit finite quality factors that depend on their position in the band gap: The closer they are to the center of the band gap, the larger their quality factor is. Specifically, the quality factors of the modes with frequencies ω_p = 1/2π = 338 GHz, ω_p = 2/2π = 382 GHz, ω_p=3/2π = 417 GHz, and ω_p = 4/2π = 468 GHz, as computed by finite difference time-domain (FDTD) simulations, are 72, 70, 6200, and 1540, respectively.

a A photo of a fabricated architecture containing the designed woodpile structure. b Side view of the unit cell of the 3D-PCC. A GaAs wafer is sandwiched by a pair of 3D photonic crystals, which consists of four layers with a rod array in the z-direction. Each silicon rod had a width of w = 0.26a and a height of h = 0.3a where a = 333 μm is the lattice constant. The rods in each successive layer of the 3D photonic crystal are rotated by 90°. The rods in the third and fourth layers are displaced by 0.5a relative to the rods in the first and second layers, respectively. The THz radiation propagates along the  −z-direction. c Schematic of a 2DEG embedded in a 3D woodpile cavity. d Transmittance spectrum of the bare 3D-PCC from numerical simulations. The shaded region corresponds to the photonic band gap computed from the band structure of an infinite woodpile lattice (in all directions) without defect. The blue arrows mark the doubly degenerate cavity modes (with respect to the two polarizations) that are labeled by p = 1, 2, 3, 4 throughout the paper. The inset is a magnified view of the spectrum.

The woodpile 3D-PCC exhibits four mirror symmetry planes per unit cell in the xy plane. In the central unit cell, those mirror symmetry planes are x = 0; ±a/2 and y = 0; ±a/2 (see Fig. 1b). Depending on the incident polarization, σ, two families of cavity modes can be distinguished due to these symmetries. The modes excited by THz light polarized along x are even (odd) with respect to y = 0; ±a/2 (x = 0; ±a/2). These modes are mostly linearly polarized in the xz-plane and are referred to as σ = x. On the other hand, the modes excited by THz light polarized along y are odd (even) with respect to y = 0; ±a/2 (x = 0; ±a/2). These modes are mostly linearly polarized in the yz-plane and are referred to as σ = y. Furthermore, the 3D-PCC features another symmetry: a combination of a 4-fold rotational symmetry with respect to the z-axis and an inversion symmetry with respect to the plane z =  0. This symmetry results in the degenerate frequency for the σ = x and σ = y modes in the spectrum²³ (see Supplementary Figs. 10–12 and Supplementary Section 3). However, the lack of inversion symmetry leads to asymmetric cavity profiles along the z-direction. For instance, Fig. 2a, b display the expectation value of the electric energy density in the vacuum state \(\left\vert 0\right\rangle\), \(\left\langle 0\right\vert {I}_{p,\sigma }(z)\left\vert 0\right\rangle\), for the cavity modes p = 1 with σ = x and σ = y. Therefore, despite having the same resonance frequency, the cavity modes excited by orthogonal polarizations exhibit different spatial profiles. Importantly, this structure is crucial in investigating the role of cavity mode profiles in multimode coupling because the cavity modes excited by orthogonal polarizations exhibit the same resonance frequency and comparable light–matter coupling strengths, with their primary distinction being the spatial profiles. It should be noted that all cavity modes are localized in the vicinity of the defect GaAs layer with a standard deviation below one Si rod thick for p = 3, 4 and about two logs thick for p = 1, 2 (see Supplementary Section 1).

Expectation value \(\left\langle 0\right\vert {I}_{p,\sigma }(z)\left\vert 0\right\rangle\) of the electric energy density along z in the vacuum state \(\left\vert 0\right\rangle\) for a the σ = x and b σ = y polarizations and the mode p = 1 (see Supplementary Section 1). The blue dashed line denotes the 2DEG location. The standard deviation \(\sqrt{\left\langle 0\right\vert {E}_{p,\sigma }^{2}({{{{\boldsymbol{\rho }}}}})\left\vert 0\right\rangle }\) of the in-plane electric field in the vacuum state (see Supplementary Section 1) at z = z_2DEG for the σ = x (c, d) and σ = y (e, f) polarizations. The upper and lower panels correspond to the mode p = 1 (ω₁/2π = 338 GHz) and p = 2 (ω₂/2π = 382 GHz), respectively. The dashed rectangle at the center illustrates the silicon rod that is right above the 2DEG layer.

To study the coupling between the cavity modes of the 3D-PCC and the CR of a 2DEG in GaAs, with dipole moment in the xy -plane, the GaAs defect layer is replaced by a GaAs substrate containing a multiple quantum well heterostructure near the top surface. The total thickness of the sample remains at around 60 μm. The standard deviation of the in-plane electric field in the electromagnetic vacuum state, \(\sqrt{\langle 0\vert {E}_{p,\sigma }^{2}({{{{\boldsymbol{\rho }}}}})\vert 0\rangle }\) (see Supplementary Section 1 for definition), for the first two bare cavity modes p = 1, 2 and σ = x, y is displayed in Fig. 2c–f at the vertical location of the 2DEG, z_2DEG. For σ = x, the cavity electric field is tightly confined at the edges of the rods (Fig. 2c, d). Conversely, the mode profiles for σ = y distribute more evenly within the unit cell (Fig. 2e, f). It should be noted that while the two families of modes σ = x, y provide a good way to classify the cavity modes without 2DEG, those two polarizations are coupled by the off-diagonal elements of the dielectric tensor describing the CR (see Supplementary Section 2). However, it turns out that such a coupling is small and the polariton eigenmodes can thus still be well identified by their polarization index σ, as shown below.

We conducted numerical simulations with COMSOL Multiphysics for the 3D-PCC with a GaAs 2DEG under an external magnetic field, B, along the z-direction (Supplementary Section 2). Figure 3a, b shows transmittance spectra for the x and y polarizations, respectively. For simplicity, we restrict ourselves to the first two modes p = 1, 2 in the following discussion, since the modes p = 3, 4 with extremely narrow line widths were not observed in the experimental data due to the limited resolution of THz spectroscopy. The cyclotron frequency, ω_c = eB/m_eff, is directly proportional to B (Supplementary Fig. 8). Four polariton branches, two upper polaritons (UPs) and two lower polaritons (LPs), are observed for the x polarization (Fig. 3a), accompanied by a splitting spanning the space around the line of the bare CR frequency (white dashed), ω_c, that separates the UPs and LPs. By contrast, this splitting is not visible for the y polarization (Fig. 3b). Instead, we observed an S-shaped polariton branch that crosses the ω_c line and behaves as an LP for cavity mode p = 2 at high B. The broad transmittance background at low frequencies is the photonic band edge, as shown in Fig. 1d.

Transmittance (T) spectra as a function of cyclotron frequency, ω_c/(2π), obtained from a, b numerical simulations and c, d experiments. The top (bottom) panels are for σ = x (σ = y) polarizations. The white dots denote the peak frequencies extracted from experimental data using longer time-domain traces (Supplementary Section 3 and Supplementary Fig. 15). The error bars of the white dots are discussed in Supplementary Section 3. The white dashed line shows the bare cyclotron frequency, ω_c/(2π) (Supplementary Fig. 8). The LP–UP splitting is marked.

We performed THz time-domain spectroscopy measurements to obtain the transmittance spectra under the same conditions. All measurements were taken at 4 K in a cryostat with a superconducting magnet. As the measurements directly collected the transmitted THz waveform in the time domain, the frequency resolution of the data is limited by the length of the time-domain traces that were used in the Fourier transformation. Note that the optical components in the system, for example, cryostat windows and crystals for THz generation and detection, induced echoes of the THz pulses at longer time delays. Including the echoes in the Fourier transformation causes artificial dips in the spectrum (the Fabry–Pérot effect)¹². Therefore, the time-domain traces used to plot the color plots in Fig. 3c, d were truncated before the arrival of THz echoes (33 ps) and zero padded to avoid the artificial dips. The frequency resolution of the color plot was limited. Furthermore, the white dots denote the B-dependent polariton frequencies that were extracted from longer time-domain traces through fittings. More details about the differences between short and long time-domain traces (Supplementary Figs. 5 and 13) and the extraction procedures (Supplementary Fig. 15) are discussed in Supplementary Section 3. The experimental results show good agreement with the simulations, namely, we observed a polariton splitting that separates the UPs and LPs for the x polarization and no splitting for the y polarization. The UPs at high frequencies were not observed in experiments likely because the transmission is low for resonances near the center of the photonic band gap and also because there exist additional losses induced by imperfections in the fabricated woodpile lattice.

Theoretical analysis

To gain a better understanding of the observed polarization-dependent transmittance spectra, we developed a microscopic quantum model described by the Hamiltonian \(\hat{H}={\hat{H}}_{{{\rm{cav}}}}+{\hat{H}}_{{{\rm{CR}}}}+{\hat{H}}_{{{{\rm{int}}}}}+{\hat{H}}_{{A}^{2}}\) (Supplementary Section 1). The free-photon Hamiltonian is \({\hat{H}}_{{{\rm{cav}}}}={\sum }_{p,\sigma } \hslash {\omega }_{p}{\hat{a}}_{p,\sigma }^{{{\dagger}} }{\hat{a}}_{p,\sigma }\), with \({\hat{a}}_{p,\sigma }\) the annihilation operator of a photon in mode p with polarization σ and frequency ω_p. The effective CR Hamiltonian \({\hat{H}}_{{{{{\rm{CR}}}}}}=\hslash {\omega }_{{{{{\rm{c}}}}}}\int\frac{d{{{{\boldsymbol{\rho }}}}}}{a}{\hat{b}}^{{{{\dagger}}} }({{{{\boldsymbol{\rho }}}}})\hat{b}({{{{\boldsymbol{\rho }}}}})\) is written in terms of the collective CR operators at the in-plane position ρ in the woodpile unit cell. The associated creation operator,

promotes an electron from the highest-occupied Landau level (LL) ν−1 with momentum k to the lowest-unoccupied LL ν and momentum k−G_y, with \(\hat{c}\) and \({\hat{c}}^{{{{\dagger}}} }\) the electron annihilation and creation operators, G_j = m_j × 2π/a (m_j = 0, 1, 2, . . . ) the j = x, y component of the in-plane reciprocal lattice vector of the 3D-PCC, \({l}_{{{{{\rm{c}}}}}}=\sqrt{\hslash /eB}\) the magnetic length, e the electron charge, ν the LL filling factor, and \({{{{\mathcal{N}}}}}\) the LL degeneracy (see Supplementary Section 1). The collective CR operators approximately satisfy bosonic commutation relations \([\hat{b}({{{{\boldsymbol{\rho }}}}}),{\hat{b}}^{{{{\dagger}}} }({{{{{\boldsymbol{\rho }}}}}}^{{\prime} })]=\delta ({{{{\boldsymbol{\rho }}}}}-{{{{{\boldsymbol{\rho }}}}}}^{{\prime} })\) in the dilute regime. We emphasize that those degenerate CR modes at each position ρ are extra degrees of freedom that are introduced in the effective Hamiltonian \({\hat{H}}_{{{{{\rm{CR}}}}}}\) because they happen to be the modes coupled to the electromagnetic field (see below).

The next term in the Hamiltonian,

is the linear coupling between the CR and the electromagnetic field, whose strength is \({g}_{p,\sigma,j}({{{{\boldsymbol{\rho }}}}})={E}_{p,\sigma,j}({{{{\boldsymbol{\rho }}}}},{z}_{{{{{\rm{2DEG}}}}}})\sqrt{{e}^{2}{\omega }_{{{{{\rm{c}}}}}}{n}_{{{{{\rm{e}}}}}}/(4{\varepsilon }_{0}{m}_{{{{{\rm{eff}}}}}}{\omega }_{p}a)}\). Here j = x, y denotes the different components of the (real, dimensionless) electric field mode functions E_p,σ,j(ρ, z_2DEG) at the in-plane position ρ and vertical location of the quantum well, m_eff is the electron mass in GaAs, ε₀ is the vacuum permittivity, and n_e is the total electron density. The coupling of the CR excitations to each cavity mode (p, σ) is weighted by the electric field mode functions since the electric field in the 3D-PCC is not uniform in the xy plane, in contrast to 1D-PCCs^11,12. Nevertheless, since the electric field mode functions are quasi-uniform over a cyclotron orbit, i.e., a ≪ l_c, the dipole approximation can still be used to calculate the coupling matrix elements.

The last term in the Hamiltonian, the so-called A² term,

is responsible for direct coupling between the cavity modes. This term scales with \({D}_{p,{p}^{{\prime} };\sigma,{\sigma }^{{\prime} }}={\sum}_{j}\int(d{{{{\boldsymbol{\rho }}}}}/{a}^{2}){g}_{p,\sigma,j}({{{{\boldsymbol{\rho }}}}}){g}_{{p}^{{\prime} },{\sigma }^{{\prime} },j}({{{{\boldsymbol{\rho }}}}})/{\omega }_{{{{{\rm{c}}}}}}\propto {\sum}_{j}\int\,d{{{{\boldsymbol{\rho }}}}}\,{E}_{p,\sigma,j}({{{{\boldsymbol{\rho }}}}},{z}_{{{{{\rm{2DEG}}}}}}){E}_{{p}^{{\prime} },{\sigma }^{{\prime} },j}({{{{\boldsymbol{\rho }}}}},{z}_{{{{{\rm{2DEG}}}}}})\), and thus, involves an overlap integral between the in-plane spatial profiles of the modes. It should be emphasized that the cavity modes are normalized over the entire volume of the 3D-PCC, i.e., \(\int\,d{{{{\boldsymbol{\rho }}}}}dz\,\varepsilon ({{{{\boldsymbol{\rho }}}}},z){{{{{\bf{E}}}}}}_{p,\sigma }({{{{\boldsymbol{\rho }}}}},z)\cdot {{{{{\bf{E}}}}}}_{{p}^{{\prime} },{\sigma }^{{\prime} }}({{{{\boldsymbol{\rho }}}}},z)= {a}^{3}{\delta }_{p,{p}^{{\prime} }}{\delta }_{\sigma,{\sigma }^{{\prime} }}\), with ε(ρ, z) the inhomogeneous dielectric profile. Thus, the cavity modes are orthogonal if the entire volume of the 3D-PCC is considered. Nevertheless, the in-plane overlap of the cavity mode profiles, which governs the mixing of the cavity modes mediated by the CR, is a priori finite, as was recently reported in THz metamaterial resonators^24,25,26. When the in-plane overlap is small, each cavity mode couples to the CR independently, and no matter-mediated intermode coupling exists. Each cavity mode leads to two polariton branches, resulting in a total of 2N branches^{24,26,27,28,29} (see Fig. 4a for N = 2). A splitting between the LPs and UPs along the bare CR frequency appears. When the spatial overlap of the cavity modes increases, the cavity modes can couple to each other through the matter, and the gap between the LPs and UPs decreases (see Fig. 4b for N = 2). When the different cavity modes within the active region highly overlap, these modes efficiently hybridize through the matter. The LP–UP splitting vanishes, producing N + 1 polariton branches with an S-shaped middle polariton (MP) and an additional dark state along the CR line (see Fig. 4c for N = 2). Furthermore, when the light–matter coupling strength is comparable or larger than the frequency difference between the cavity modes, the system enters the superstrong coupling (SSC) regime^{30,31,32,33,34,35,36}, and the mixing between different cavity modes increases, especially close to the inflection point of the S-shaped MP (Fig. 4b).

a–c Sketch of the typical polariton branches (black solid lines) as a function of the magnetic field for two cavity modes (N = 2, red dashed lines). The CR is depicted as a blue dashed line, ω_c. a When the spatial overlap between the cavity mode profiles is low, the cavity modes coupled independently to the CR, resulting in 2N polariton branches (“decoupled scenario”). A splitting occurs between the LPs and UPs. b As the spatial overlap increases, the cavity modes start mixing with each other via the CR, and the splitting gradually decreases. c When the cavity mode profiles exhibit high spatial overlap, the splitting vanishes, leading to N + 1 polariton branches and an additional dark state along the ω_c line (“coupled scenario”). In this case, UP₁ and LP₂ merge to form an MP, which is a mixture of the two cavity modes close to the inflection point. UP upper polariton, LP lower polariton, MP middle polariton.

To better understand how the two regimes with either mode coupling or mode decoupling emerge from the microscopic description, we introduce new CR excitation operators following the spatial profiles of the cavity modes, \({\hat{b}}_{p,\sigma }=\int(d{{{{\boldsymbol{\rho }}}}}/a)\,({\widetilde{g}}_{p,\sigma }({{{{\boldsymbol{\rho }}}}})/{\Omega }_{p,\sigma })\hat{b}({{{{\boldsymbol{\rho }}}}})\), with \({\widetilde{g}}_{p,\sigma }({{{{\boldsymbol{\rho }}}}})={g}_{p,\sigma,y}({{{{\boldsymbol{\rho }}}}})-i{g}_{p,\sigma,x}({{{{\boldsymbol{\rho }}}}})\) and the effective coupling strength Ω_p,σ defined by \({\Omega }_{p,\sigma }^{2}=\int(d{{{{\boldsymbol{\rho }}}}}/{a}^{2})\,| {\widetilde{g}}_{p,\sigma }({{{{\boldsymbol{\rho }}}}}){| }^{2}\). By construction, the linear coupling term (2) written in terms of those operators takes the simple form

The commutation relations between the new CR modes read \([{\hat{b}}_{p,\sigma },{\hat{b}}_{{p}^{{\prime} },{\sigma }^{{\prime} }}^{{{{\dagger}}} }]={\xi }_{p,{p}^{{\prime} };\sigma,{\sigma }^{{\prime} }}\), where

is proportional to the in-plane spatial overlap of the different cavity modes, similarly to the off-diagonal contributions of the A² term. If the new CR modes \({\hat{b}}_{p,\sigma }\) were orthogonal, which is a priori not the case, one would have \({\xi }_{p,{p}^{{\prime} };\sigma,{\sigma }^{{\prime} }}={\delta }_{p,{p}^{{\prime} }}{\delta }_{\sigma,{\sigma }^{{\prime} }}\). The deviation of the parameter \({\xi }_{p,{p}^{{\prime} };\sigma,{\sigma }^{{\prime} }}\) with respect to \({\delta }_{p,{p}^{{\prime} }}{\delta }_{\sigma,{\sigma }^{{\prime} }}\), therefore, is a measure of how much coupling between the cavity modes is mediated by the CR.

First, we find that \({\xi }_{p,{p}^{{\prime} };x,y} \, < \, 0.36\) for all \(p,{p}^{{\prime} }\), indicating that the cavity modes with orthogonal polarizations are weakly coupled by the CR. Thus, the polarization index σ remains quite a good quantum number to classify the polariton eigenmodes, as stated before. For a given polarization σ, \({\xi }_{p,{p}^{{\prime} };\sigma,\sigma }=1\) (\({\xi }_{p,{p}^{{\prime} };j,j}=0\)) corresponds to perfect overlap (no overlap) between the cavity modes p and \({p}^{{\prime} }\), in which case the cavity modes are coupled (decoupled) via the CR. For our 3D-PCC system, we find a rather small overlap between the spatial profiles of the σ = x modes p = 1 and p = 2 (ξ_1,2;x,x = 0.29), suggesting that the intermode coupling is weak. Conversely, for σ = y modes, the spatial profiles of cavity modes p = 1 and p = 2 highly overlap (ξ_1,2;y,y = 0.91), suggesting that they are strongly coupled through the CR.

When \({\xi }_{p,{p}^{{\prime} };\sigma,{\sigma }^{{\prime} }}={\delta }_{p,{p}^{{\prime} }}{\delta }_{\sigma,{\sigma }^{{\prime} }}\), and neglecting the coupling between different polarizations induced by the A² term, one has \({D}_{p,{p}^{{\prime} };\sigma,\sigma }\approx ({\Omega }_{p,\sigma }{\Omega }_{{p}^{{\prime} },\sigma }/{\omega }_{c}){\xi }_{p,{p}^{{\prime} };\sigma,\sigma }=({\Omega }_{p,\sigma }^{2}/{\omega }_{c}){\delta }_{p,{p}^{{\prime} }}\), and the full Hamiltonian can then be simplified to a “decoupled” Hamiltonian of the form \(\hat{H}={\hat{H}}_{{{{{\rm{cav}}}}}}+{\hat{H}}_{{{{{\rm{CR}}}}}}+{\hat{H}}_{{{{{\rm{int}}}}}}+{\hat{H}}_{{A}^{2}}\), with \({\hat{H}}_{{{\rm{cav}}}}={\sum}_{p,\sigma } \hslash {\omega }_{p}{\hat{a}}_{p,\sigma }^{{{\dagger}}} {\hat{a}}_{p,\sigma }\), \({\hat{H}}_{{{{{\rm{CR}}}}}}=\hslash {\omega }_{{{{{\rm{c}}}}}}{\sum}_{p,\sigma }{\hat{b}}_{p,\sigma }^{{{{\dagger}}} }{\hat{b}}_{p,\sigma }\), \({\hat{H}}_{{{{{\rm{int}}}}}}\) given by Eq. (4), and

In contrast to the full Hamiltonian, such a decoupled Hamiltonian can be diagonalized in each subspace (p, σ) independently²⁸.

We computed transmittance spectra by extending the input–output model of ref. ³⁷ in a simple planar geometry to a multimode PCC cavity (Supplementary Section 1). Dissipation is included through the (Markovian) coupling of cavity modes and CR excitations to phenomenological bosonic reservoirs, with associated quality factors of the bare cavity modes computed by FDTD, and of the intrinsic CR decay rate, respectively (Supplementary Section 1). The input–output model allows to select the polarization (x, y) of the input THz field to probe the polariton eigenmodes for both σ = x and σ = y. The calculated spectra using the full Hamiltonian (Fig. 5a, b) are consistent with the simulations and experimental results (Fig. 3), except for the transmittance background at low frequencies due to the lower edge of the photonic band gap that is not included in the model. In particular, the model replicates the existence of an LP–UP splitting close to the ω_c line and the emergence of an S-shaped MP that adiabatically transforms from the UP of mode 1 to the LP of mode 2 for σ = x and σ = y modes, respectively. Note that since the input field in the model is assumed to be independent of the cavity mode it couples to (Supplementary Section 1), the transmittance at high frequencies is overestimated as compared to the experimental data and numerical simulations. Since ξ_1,2;x,x (σ = x) is small, the decoupled model provides a faithful description of the system, and thus the calculated spectra (Fig. 5c) are almost the same as that using the full Hamiltonian (Fig. 5a). For σ = y, however, the decoupled model is inaccurate as it fails to predict the emergence of the S-shaped MP (Fig. 5d). This is due to the large overlap ξ_1,2;y,y = 0.91, which means that the CR operators \({\hat{b}}_{1,y}\) and \({\hat{b}}_{1,y}^{{{{\dagger}}} }\) are not orthogonal, i.e., \([{\hat{b}}_{1,y},{\hat{b}}_{2,y}^{{{{\dagger}}} }] \, \ne \, 0\).

Transmittance (T) spectra as a function of cyclotron frequency, ω_c/(2π), obtained from calculations using a, b the full Hamiltonian and c, d the decoupled Hamiltonian, respectively (Supplementary Section 1). The top (bottom) panels are for σ = x (σ = y) modes. The white dots denote the peak frequencies extracted from experimental data using longer time-domain traces (Supplementary Section 3 and Supplementary Fig. 15). The error bars of the white dots are discussed in Supplementary Section 3. The white dashed line represents ω_c/(2π). The yellow arrows in (b, d) mark the UP that significantly deviates from the experimental data points in the calculation using the decoupled Hamiltonian. Calculated ground-state (virtual) correlations \(\langle {\hat{a}}_{p,\sigma }^{{{{\dagger}}} }{\hat{a}}_{{p}^{{\prime} },\sigma }\rangle\) using the full Hamiltonian for e σ = x and f σ = y modes. The intermode ground-state correlations between the photonic modes p = 1 and p = 2 are significant for the σ = y modes.

In order to assess in which regime of light–matter coupling our system operates, we extracted the coupling strengths for each cavity mode from the theoretical model and evaluated the relevant figures of merit. The USC regime, which has attracted widespread interest due to the compelling prospect of manipulating basic material properties by engineering the vacuum electromagnetic field surrounding the material inside a cavity^12,38,39, is achieved when the effective coupling strength for each cavity mode Ω_p,σ reaches about 10% of the bare mode frequency at resonance^18,19. Here we find Ω_1,x/ω₁ ≈ 0.21 and Ω_1,y/ω₁ ≈ 0.2 at B = 0.81 T, and Ω_2,x/ω₂ ≈ 0.21 and Ω_2,y/ω₂ = 0.17 at B = 0.92 T, respectively. This shows that our system operates in the USC regime for both cavity modes and polarization σ.

It is interesting to investigate unique quantum features that can emerge in our multimode system. Since the ground state in the USC regime is known to be a dressed ground state that has a finite photon population, \(\langle {n}_{p,\sigma }\rangle=\langle {\hat{a}}_{p,\sigma }^{{{{\dagger}}} }{\hat{a}}_{p,\sigma }\rangle\)^37,40, an intriguing question is whether the virtual photons (cavity vacuum fields) in different modes can mix through simultaneous USC with the matter, resulting in finite intermode correlations in the ground state \(\langle {\hat{a}}_{p,\sigma }^{{{{\dagger}}} }{\hat{a}}_{{p}^{{\prime} },\sigma }\rangle \, > \, 0\). We computed the photon ground-state correlations \(\langle {\hat{a}}_{p,\sigma }^{{{{\dagger}}} }{\hat{a}}_{{p}^{{\prime} },\sigma }\rangle\) for both σ = x, y by inverting the Bogoliubov transformation that diagonalizes the full (quadratic) Hamiltonian. The photon number for each mode p and polarization σ in the ground state (without external driving) is shown in Fig. 5e. While \(\langle {\hat{a}}_{p,\sigma }^{{{{\dagger}}} }{\hat{a}}_{p,\sigma }\rangle\) reaches its maximum at resonance between the CR and the cavity mode (p, σ), the intermode correlation \(\langle {\hat{a}}_{1,\sigma }^{{{{\dagger}}} }{\hat{a}}_{2,\sigma }\rangle\) is peaked in between the two resonances. For σ = x, we find that \(\langle {\hat{a}}_{1,x}^{{{{\dagger}}} }{\hat{a}}_{2,x}\rangle \ll \langle {\hat{a}}_{p,x}^{{{{\dagger}}} }{\hat{a}}_{p,x}\rangle\). However, for σ = y, \(\langle {\hat{a}}_{1,y}^{{{{\dagger}}} }{\hat{a}}_{2,y}\rangle\) is comparable to \(\langle {\hat{a}}_{p,y}^{{{{\dagger}}} }{\hat{a}}_{p,y}\rangle\) (Fig. 5f), which suggests that the mixing of virtual photons in different photonic modes is also directly related to the spatial overlap of the cavity profiles. Those intermode vacuum correlations should be accessible in our system from, e.g., first-order equal-time correlation measurements.

Based on the scaling of the intermode ground-state correlations with the overlap between the cavity modes and the mode frequencies (see Supplementary Fig. 1b), we extend the standard figure of merit for single-mode USC to the multimode case by introducing the parameter

which does not depend on the magnetic field and whose diagonal elements (\(p={p}^{{\prime} }\)) coincide with the standard figure of merit. While intermode (off-diagonal) coupling for σ = x is smaller than the diagonal coupling (η_12,x = 0.11 as compared to η_11,x ≈ η_22,x ≈  0.21), we find that intermode coupling for σ = y is of the same magnitude (η_12,y = 0.17 as compared to η_11,y ≈ 0.19 and η_22,y ≈  0.16), consistently with the intermode ground-state correlations.

On the other hand, the SSC regime is achieved when the coupling strength becomes comparable to the frequency difference between the uncoupled modes, thereby signaling the breakdown of the single-mode approximation^{30,31,32,33,34,35,36}. In this regime, the hybridization of cavity mode profiles³⁰ and complex multimode dynamics^32,41 have been discussed. However, our results highlight that such a criterion is not sufficient to characterize the SSC regime because the mixing between cavity modes is also influenced by the spatial profiles of the cavity modes. For example, two independent sets of LPs and UPs are observed for σ = x as a result of the weak overlap of the cavity modes. Thus, the single-mode approximation still holds in this case even though the standard criterion is satisfied. As explained above, this “mode decoupling” situation would generally occur in all systems where quantum emitters fill the entire cavity mode volume^{24,26,27,28,29}. To circumvent this issue, we extend the standard figure of merit for the SSC regime to take into account the spatial overlap between the cavity modes,

This parameter quantifies the strength of the intermode coupling with respect to the frequency difference ω₂−ω₁, including the effect of the cavity mode overlap. This parameter is B-independent as ω_c is included in the denominator, which balances the ω_c term in the expression for \({\tilde{g}}_{1,\sigma }{\tilde{g}}_{2,\sigma }^{*}\). We find that Λ_x ≈ 0.33 and Λ_y ≈ 0.49, which indicates that the polarization σ = y is further in the SSC regime than the polarization σ = x. Deep into the SSC regime (Λ_σ ~ 1), the S-shaped MP would become a hybrid mode equally composed of the two cavity modes. We show in Supplementary Section 1 that the SSC figure of merit governs the weight of the MP onto the different cavity modes (Supplementary Fig. 1c), which in turn governs the intermode correlations in polaritonic excited states.