Perspective on real-space nanophotonic field manipulation using non-perturbative light-matter coupling

Confining light below the Abbe diffraction limit [1] by storing a part of the electromagnetic energy in the kinetic energy of electric charges [2] opened the door to a number of groundbreaking real-world applications which has contributed to the great success of the field of nanophotonics. In a nanophotonic device, the high energy density of the electromagnetic field makes it relatively easy to couple with different kinds of localised material excitations and reach the strong light-matter coupling regime, originally achieved in cavity quantum electrodynamics (CQED) atomic systems [3]. In such a regime, light and matter degrees of freedom hybridise, leading to novel, polaritonic excitations of mixed light-matter character [4, 5].

Standard theoretical models used to describe strong coupling consider a single optically active matter transition coupled to a single photonic mode. Although some care has to be used when performing calculations on such a reduced Hilbert space [6, 7, 8, 9], this single-mode approximation has enabled modelling of a wide range of CQED systems with remarkable easiness and generality. However, the requirement is that the energy spacing between the considered resonances and the neglected ones is much larger than the strength of light-matter coupling, thus permitting to integrate out excited modes with negligible populations.

However, the ongoing race for record coupling strengths [10, 11] has led to situations in which higher-energy electronic states cannot be neglected, requiring a model which considers the coupling of multiple matter excitations to the same photonic mode. We refer to this regime as the very-strong coupling (VSC) regime, first predicted by Khurgin in 2001 [12]. The hybridization of multiple excited matter states has an important consequence: the matter component of the polariton, represented itself by a linear superposition of different bare matter wavefunctions, has a wavefunction different from each of the bare states [13]. Following a 2013 proposal [14], such an effect was observed for the first time in 2017 [15], as a modification of approximately $30\%$ of the Wannier exciton Bohr radius in GaAs microcavities. Larger numbers of matter states which can be hybridised by the coupling with the photonic field could correspond to a broader design space for the resulting electronic wavefunction. This idea led to the study of systems with a continuum of ionised excitations [16, 17] and eventually to the discovery of novel bound excitons stabilised by the photonic interaction [18], and to novel polaritonic loss channels [19].

In this Article we theoretically investigate the possibility of both multi-mode electronic as well as multi-mode photonic hybridization, leading to a modification of the spatial electromagnetic profiles of the resulting polariton modes. Given the possibility of fast [20, 21] in-situ tuning of the light-matter interaction by optical and electrical means, subcycle multiwave mixing nonlinearities between different polariton states [22] or even all-optical subcycle switching [23, 24], such an approach could open the door to dynamical manipulation of subwavelength fields, with potential disruptive applications for, e.g., on-chip optical tweezers [25].

Although to the best of our knowledge it was never explicitly discussed, the regime of photonic VSC has been already reached in various systems, as superconducting qubits coupled to microwave photons in a long transmission-line resonator [26, 27]. Moreover, it has been theoretically [28] and experimentally [29] demonstrated in microcavities, where the coupling strength becomes larger than the bare excitation frequencies. In such a regime, the diamagnetic term of the Hamiltonian creates a dominant real-space repulsive interaction localised at the dipole position, which expels the electromagnetic field and may even lead to light-matter decoupling [28, 29]. It has also been experimentally observed that, in plasmonic nanocavities, the greatly enhanced coupling between molecular excitons and gap plasmons causes a significant modification of the plasmonic modes profile [30].

Here we focus on Landau polaritons, where the giant electronic dipoles of cyclotron resonances (CRs) of two-dimensional electron gases (2DEGs) are coupled to strongly enhanced light fields of subwavelength THz resonators. After inital predictions in Ref. [31], multiple experimental realizations followed, some of which established world-records for the largest light-matter coupling ever achieved in any CQED system [32, 33, 29, 34].

In the first part of the paper we will develop a general theory describing multi-mode light-matter coupling of CRs of a planar 2DEG which is Landau-quantised by a perpendicular static magnetic field, and multiple photonic resonances. Our approach highlights the main electronic and optical features observable for this multi-mode coupling. In the second part, we apply our formalism to structures based on planar plasmonic metasurfaces. To this end, we perform numerical simulations which verify the predictions of our theory and demonstrate how multi-mode photonic hybridization can lead to a modification of the electromagnetic spatial profile of the polariton modes.

In this section we develop a theory for the light-matter coupling between $M$ photonic resonator modes and the CRs of a 2DEG with a charge carrier density $N_{\textrm{2DEG}}$ and an effective mass $m^{*}$. Following Kohn’s theorem [35], we neglect Coulomb interactions between the electrons, which manifest in the nonlinear susceptibility of strongly driven Landau electron systems [36], but have no role in the determination of the optical resonance. Moreover, while our theory focuses on a single quantum well (QW) hosting the entire electron distribution, it is equally applicable to densely packed multi-QW structures as usually employed in experiments, where the intensity of the electromagnetic field doesn’t vary significantly within the thickness of the multi-QW stack. Following the elegant theory from Ref. [34] we can write the Hamiltonian of our system as

	$\displaystyle\hat{H}$	$\displaystyle=$	$\displaystyle\hat{H}_{\text{cav}}+\sum_{j=1}^{N}\hbar\omega_{c}\hat{c}_{j}^{\dagger}\hat{c}_{j}+\frac{e^{2}}{m^{*}}\sum_{j=1}^{N}\hat{A}_{-}(\mathbf{r}_{j})\hat{A}_{+}(\mathbf{r}_{j})$		(1)
		$\displaystyle+$	$\displaystyle i\sqrt{\frac{\hbar\omega_{c}e^{2}}{m^{*}}}\sum_{j=1}^{N}\left[\hat{c}_{j}^{\dagger}\hat{A}_{+}(\mathbf{r}_{j})-\hat{c}_{j}\hat{A}_{-}(\mathbf{r}_{j})\right],$

where $\hat{c}_{j}$ is the bosonic lowering operator for the electrons, leading to a transition from the $j$th to the $(j-1)$th Landau level with a transition energy $\hbar\omega_{c}$, and $\hat{H}_{\text{cav}}$ is the Hamiltonian describing the bare electromagnetic field in the resonator. In case of high electron density and strong in-plane confinement of both the 2DEG and the electromagnetic field, plasmonic modes hosted by the system can be non-negligible and lead to the formation of magneto-plasmon modes with a renormalized frequency of $\tilde{\omega}_{c}=\sqrt{\omega_{c}^{2}+\omega_{P}^{2}}$, where $\omega_{P}$ is the 2D plasmon frequency for the 2DEG [19]. However, a correct estimation of $\omega_{P}$ not only takes into account the in-plane confinement of the 2DEG, but also includes the screening of the metallic resonator in proximity of the electrons, leading to a reduction of the plasmon energy [37, 38, 39]. For our structures, this effect strongly limits the extent of renormalization such that we disregard plasmon effects.

In Eq. 1 we introduced the non-Hermitian vector potentials written in terms of the in-plane component of the vector potential $\hat{\mathbf{A}}(\mathbf{r})$ as

$\displaystyle\hat{A}_{\pm}(\mathbf{r})$

$\displaystyle=$

$\displaystyle\frac{\hat{A}_{x}(\mathbf{r})\mp i\hat{A}_{y}(\mathbf{r})}{\sqrt{2}}.$

(2)

The full vector potential can be expressed as a sum of photonic modes obtained by solving Maxwell’s equations for the bare cavity, with dimensionless spatial field profiles $\bf{f}_{\nu}(\textbf{r})$, frequencies $\omega_{\nu}$, and second-quantized bosonic annihilation operators $\hat{a}_{\nu}$ as

$\displaystyle\hat{\mathbf{A}}(\textbf{r})$

$\displaystyle=$

$\displaystyle\sum_{\nu}\sqrt{\frac{\hbar}{2\epsilon_{0}\epsilon_{r}(\mathbf{r})\omega_{\nu}\mathbb{V}_{\nu}}}\mathbf{f}_{\nu}(\textbf{r})\left(\hat{a}_{\nu}^{\dagger}+\hat{a}_{\nu}\right).$

(3)

Here, the classical modes are orthogonal over the full integration domain $\mathbb{V}$ [40]

$\displaystyle\int_{\mathbb{V}}\bf{f}^{*}_{\nu}(\mathbf{r})\bf{f}_{\mu}(\mathbf{r})d\mathbf{r}=\mathbb{V}_{\nu}\delta_{\nu,\mu},$

(4)

with $\mathbb{V}_{\nu}$ the mode volume of the $\nu$th photon mode and $\epsilon_{r}(\mathbf{r})$ the background, non-resonant dielectric constant. The amplitudes of the non-Hermitian vector potentials then take the form

	$\displaystyle\hat{A}_{-}(\textbf{r})$	$\displaystyle=$	$\displaystyle\sum_{\nu}\sqrt{\frac{\hbar}{2\epsilon_{0}\epsilon_{r}(\mathbf{r})\omega_{\nu}\mathbb{V}_{\nu}}}f_{\nu}(\textbf{r})\left(\hat{a}_{\nu}^{\dagger}+\hat{a}_{\nu}\right),$		(5)
	$\displaystyle\hat{A}_{+}(\textbf{r})$	$\displaystyle=$	$\displaystyle\sum_{\nu}\sqrt{\frac{\hbar}{2\epsilon_{0}\epsilon_{r}(\mathbf{r})\omega_{\nu}\mathbb{V}_{\nu}}}f^{*}_{\nu}(\textbf{r})\left(\hat{a}_{\nu}^{\dagger}+\hat{a}_{\nu}\right),$

with

$\displaystyle f_{\nu}(\mathbf{r})$

$\displaystyle=$

$\displaystyle\frac{f_{\nu,x}(\mathbf{r})+if_{\nu,y}(\mathbf{r})}{\sqrt{2}}.$

(6)

Crucially, the orthogonality condition in Eq. 4 holds only if the integral is performed over the entire three-dimensional space, while the integral of two orthogonal modes over any sub-domain does not vanish in general. This concept is illustrated in Fig. 1 for the model case of a planar microcavity (a,b) and for a split-ring resonator (c,d), integrated over either the full three-dimensional volume (a,c) or a thin, quasi-two-dimensional surface (b,d). In both cases, two orthogonal modes (red and blue arrows) become non-orthogonal when the integral is performed over a quasi-two-dimensional slice of the overall volume. This finding is relevant for our systems because the third term of Eq. 1, the so-called diamagnetic term of the light-matter interaction Hamiltonian, contains generally non-vanishing expressions of the form

$\displaystyle\sum_{j=1}^{N}f^{*}_{\nu}(\mathbf{r}_{j})f_{\mu}(\mathbf{r}_{j})$

$\displaystyle=$

$\displaystyle N_{\textrm{2DEG}}\int_{\mathbb{S}}f^{*}_{\nu}(z,\mathbf{r}_{\|})f_{\mu}(z,\mathbf{r}_{\|})d\mathbf{r}_{\|},$

where $\mathbb{S}$ is the sample surface, $z$ is the out-of-plane position of the 2DEG and $\mathbf{r}_{\|}$ is the in-plane position. Placing a 2DEG at the center of the planar microcavity, or below the split-ring resonator, will thus result in an interaction of different photon modes which is mediated and modulated by the coupling to the electrons.

We elucidate this insight by developing the first $M$ photonic functions $f_{\nu}(z,\textbf{r}_{\|})$ restricted over the quantum well plane as linear superpositions of at most $M$ orthonormal basis functions over the sample surface $\mathbb{S}$,

$\displaystyle\int_{\mathbb{S}}\phi^{*}_{\nu}(\textbf{r}_{\|})\phi_{\mu}(\textbf{r}_{\|})d\mathbf{r}_{\|}=\delta_{\nu,\mu}$

(8)

and obtain

$\displaystyle f_{\nu}(z,\textbf{r}_{\|})=\sum_{\mu\leq\nu}\alpha_{\nu,\mu}\phi_{\mu}(\textbf{r}_{\|}).$

(9)

It is always possible to choose the basis such that $\alpha_{1,1}$ is real and $\alpha_{\nu,\mu}=0$ if $\nu<\mu$. Using Eq. 9 the degree of non-orthonormality between the resonator modes with respect to the quantum well plane can be captured by defining the overlap matrix

$$\mathcal{F}_{\nu,\mu}=\int_{\mathbb{S}}f^{*}_{\nu}(z,\textbf{r}_{\|})f_{\mu}(z,\textbf{r}_{\|})d\mathbf{r}_{\|}=\sum_{\gamma\leq\min(\nu,\mu)}\alpha^{*}_{\nu,\gamma}\alpha_{\mu,\gamma},$$

(10)

and its normalized version

$$\eta_{\nu,\mu}=\frac{\mathcal{F}_{\nu,\mu}}{\sqrt{\mathcal{F}_{\mu,\mu}\mathcal{F}_{\nu,\nu}}},$$

(11)

both of which may assume values from $0$ to $1$. These matrices quantify the spatial overlap of any pair ($\mu,\nu$) of photon modes over the quantum well plane. A diagonal matrix $\eta_{\nu,\mu}\propto\delta_{\nu,\mu}$ implies vanishing overlap between the photon modes, while a fully populated matrix corresponds to a strong overlap.

By introducing a set of collective bosonic matter operators

$\displaystyle\hat{b}_{\mu}=\frac{1}{\sqrt{N_{\text{2DEG}}}}\sum_{j=1}^{N}\phi_{\mu}(\textbf{r}_{j,\|})\hat{c}_{j},$

(12)

with the in-plane position $\mathbf{r}_{j,\|}$ of the $j$th electron, we can finally write the Hamiltonian in Coulomb Gauge as

	$\displaystyle\hat{H}$	$\displaystyle=$	$\displaystyle\sum_{\nu}\hbar\omega_{\nu}\hat{a}_{\nu}^{\dagger}\hat{a}_{\nu}+\sum_{\nu}\hbar\omega_{c}\hat{b}_{\nu}^{\dagger}\hat{b}_{\nu}$
		$\displaystyle+$	$\displaystyle\sum_{\nu,\mu}h_{\nu,\mu}\left(\hat{a}_{\nu}^{\dagger}+\hat{a}_{\nu}\right)\left(\hat{a}_{\mu}^{\dagger}+\hat{a}_{\mu}\right)$
		$\displaystyle+$	$\displaystyle\sum_{\nu}\sum_{\mu\leq\nu}\left[\left(g_{\nu,\mu}b_{\mu}+g^{*}_{\nu,\mu}b_{\mu}^{\dagger}\right)\left(\hat{a}_{\nu}^{\dagger}+\hat{a}_{\nu}\right)\right].$

Here,

	$\displaystyle g_{\nu,\mu}$	$\displaystyle=$	$\displaystyle\alpha_{\nu,\mu}\sqrt{\frac{\hbar^{2}\omega_{c}N_{\text{2DEG}}e^{2}}{2m^{*}\epsilon_{0}\bar{\epsilon}_{r}\omega_{\nu}\mathbb{V}_{\nu}}},$		(14)
	$\displaystyle h_{\nu,\mu}$	$\displaystyle=$	$\displaystyle\sum_{\gamma\leq\nu,\mu}\frac{g_{\nu,\gamma}g_{\mu,\gamma}}{\hbar\omega_{c}},$

represent coupling parameters consisting of the background dielectric constant $\bar{\epsilon}_{r}$ of the quantum well material. This Hamiltonian presents some important features. First, it is bosonic and quadratic, which allows us to determine its eigenmodes by Hopfield diagonalization [41]. Second, the interaction term displays cross-interactions between different photon and matter modes.

Following the Hopfield approach, we diagonalise the Hamiltonian by introducing the hybrid multi-mode polariton operators,

$\displaystyle\hat{p}_{\mu}$

$\displaystyle=$

$\displaystyle\sum_{\nu}\left(x_{\nu,\mu}\hat{a}_{\nu}+w_{\nu,\mu}\hat{b}_{\nu}+y_{\nu,\mu}\hat{a}_{\nu}^{\dagger}+z_{\nu,\mu}\hat{b}_{\nu}^{\dagger}\right),$

(15)

whereby $\left(x_{\nu,\mu},w_{\nu,\mu},y_{\nu,\mu},z_{\nu,\mu}\right)$ are real-valued Hopfield coefficients. The dressed polariton frequencies $\omega^{p}_{\mu}$ are the eigenvalues of the polariton eigenequation

$\displaystyle\hbar\omega_{\mu}^{p}\hat{p}_{\mu}$

$\displaystyle=$

$\displaystyle\left[\hat{p}_{\mu},\hat{H}\right].$

(16)

The Hopfield transformation can subsequently be inverted as

$\displaystyle\left(\hat{a}_{\nu}+\hat{a}_{\nu}^{\dagger}\right)$

$\displaystyle=$

$\displaystyle\sum_{\mu}\left(x_{\nu,\mu}-y_{\nu,\mu}\right)\left(\hat{p}_{\mu}+\hat{p}_{\mu}^{\dagger}\right),$

(17)

allowing us to find the coupled electric field components corresponding to the non-Hermitian vector potential

	$\displaystyle\hat{E}_{-}(\textbf{r})$	$\displaystyle=$	$\displaystyle\sum_{\nu,\mu}\sqrt{\frac{\hbar\omega_{\nu}}{2\epsilon_{0}\bar{\epsilon}_{r}\mathbb{V}_{\nu}}}f_{\nu}(\textbf{r})\left(x_{\nu,\mu}-y_{\nu,\mu}\right)\left(\hat{p}_{\mu}^{\dagger}+\hat{p}_{\mu}\right),$
	$\displaystyle\hat{E}_{+}(\textbf{r})$	$\displaystyle=$	$\displaystyle\sum_{\nu,\mu}\sqrt{\frac{\hbar\omega_{\nu}}{2\epsilon_{0}\bar{\epsilon}_{r}\mathbb{V}_{\nu}}}f^{*}_{\nu}(\textbf{r})\left(x_{\nu,\mu}-y_{\nu,\mu}\right)\left(\hat{p}_{\mu}^{\dagger}+\hat{p}_{\mu}\right).$

From Eq. II we can clearly see that, as expected from our initial discussion, the electric field corresponding to the polaritonic mode $\hat{p}_{\mu}$ is a linear combination of all bare electromagnetic mode profiles $f_{\nu}(\textbf{r})$, each weighted by the Hopfield coefficients.

In order to highlight the role of the normalised overlap factors for the coupling strength, we now assume a single pair of photonic modes ($M=2$) with frequencies $\omega_{1}$ and $\omega_{2}$ and mode volumes $\mathbb{V}_{1}$ and $\mathbb{V}_{2}$. Their non-orthogonality is quantified by a single overlap parameter $\eta_{2,1}$. By expliciting Eq. 10, we arrive at

	$\displaystyle\mathcal{F}_{1,1}$	$\displaystyle=$	$\displaystyle\alpha_{1,1}^{2},$
	$\displaystyle\mathcal{F}_{2,2}$	$\displaystyle=$	$\displaystyle|\alpha_{2,1}|^{2}+|\alpha_{2,2}|^{2},$		(19)
	$\displaystyle\mathcal{F}_{2,1}$	$\displaystyle=$	$\displaystyle\alpha^{*}_{2,1}\alpha_{1,1},$

which leads to

	$\displaystyle\alpha_{1,1}$	$\displaystyle=$	$\displaystyle\sqrt{\mathcal{F}_{1,1}},$
	$\displaystyle\alpha_{2,1}$	$\displaystyle=$	$\displaystyle\frac{\mathcal{F}_{2,1}^{*}}{\sqrt{\mathcal{F}_{1,1}}}=\sqrt{\mathcal{F}_{2,2}}\eta^{*}_{2,1},$		(20)
	$\displaystyle\alpha_{2,2}$	$\displaystyle=$	$\displaystyle\sqrt{\mathcal{F}_{2,2}-\frac{|\mathcal{F}_{2,1}|^{2}}{\mathcal{F}_{1,1}}}=\sqrt{\mathcal{F}_{2,2}}\sqrt{1-|\eta_{2,1}|^{2}}.$

Defining the renormalised mode volume as $\tilde{\mathbb{V}}_{\nu}=\frac{\mathbb{V}_{\nu}}{\mathcal{F}_{\nu,\nu}}$, Eq. 14 leads to expressions for the coupling strengths

	$\displaystyle g_{1,1}$	$\displaystyle=$	$\displaystyle\sqrt{\frac{\hbar^{2}\omega_{c}N_{\text{2DEG}}e^{2}}{2m^{*}\epsilon_{0}\bar{\epsilon}_{r}\omega_{1}\tilde{\mathbb{V}}_{1}}},$		(21)
	$\displaystyle g_{2,1}$	$\displaystyle=$	$\displaystyle\sqrt{\frac{\hbar^{2}\omega_{c}N_{\text{2DEG}}e^{2}}{2m^{*}\epsilon_{0}\bar{\epsilon}_{r}\omega_{2}\tilde{\mathbb{V}}_{2}}}\eta_{2,1},$		(22)
	$\displaystyle g_{2,2}$	$\displaystyle=$	$\displaystyle\sqrt{\frac{\hbar^{2}\omega_{c}N_{\text{2DEG}}e^{2}}{2m^{*}\epsilon_{0}\bar{\epsilon}_{r}\omega_{2}\tilde{\mathbb{V}}_{2}}}\sqrt{1-|\eta_{2,1}|^{2}}.$		(23)

For the given basis, the interpretation of these coefficients is that the photonic mode $\nu=1$ is coupled to only a single matter mode, $\mu=1$. In contrast, the coupling strength for the photonic mode $\nu=2$ originates from simultaneous coupling to both matter modes owing to the non-vanishing overlap parameter $\eta_{2,1}$.

In order to show the peculiar spectroscopic features expected in systems with non-negligible overlap between the photonic modes, we plot in Fig. 2 the spectra obtained by diagonalising the Hamiltonian in Eq. II for two resonator modes. The three cases concern settings of vanishing overlap ($\eta_{2,1}=0$, panels a,b), medium overlap ($\eta_{2,1}=0.5$, panels c,d), and maximum overlap ($\eta_{2,1}=1$, panels e,f), whereby in each case the left and right panel show spectra as a function of the cyclotron frequency, and electron density $N_{\textrm{2DEG}}$, respectively.

We can point out two characteristic signatures for the overlap. First, we consider varying the cyclotron frequency (panels a,c,e). For vanishing mode overlap $\eta_{2,1}=0$ (panel a), we observe the opening of separate polariton gaps for each pair of photonic mode and matter excitation. On the contrary, maximum overlap of $\eta_{2,1}=1$ (panel e) leads to the emergence of an S-shaped resonance (blue curve). In this case, the mode structure originates from the coupling of a single matter excitation $\mu=1$ to both photonic modes $\nu=1,2$, simultaneously, leading to three polariton branches in total. The S-shaped center mode is confined between the cavity frequencies $\omega_{1}$ and $\omega_{2}$, thus manifesting a double-mode nature. Second, we analyze the mode structure as a function of electron density (panels b,d,f). Here, we see that at larger densities and thus larger couplings, two modes blue-shift in the case of vanishing overlap, while a single mode blue-shifts in the presence of substantial overlap. We attribute this behavior to the contribution of the diamagnetic term which, being of higher order in $N_{\textrm{2DEG}}$, becomes dominant at very large densities and tends to blue-shift the upper polariton of each set of polaritonic solutions, taking into account that polaritonic modes never cross their bare components [28, 42]. Nevertheless, in the case of maximum overlap, the diamagnetic term between the two photonic modes leads to a repulsion of the upper polaritons, leading to an anti-crossing behaviour above a certain critical value of the electronic density.

In order to explore the relevance of our theory for experiments with Landau polaritons, we performed numerical FEM simulations of a metamaterial coupled to a 2DEG, following Bayer et al. [29]. To compute the complex field distribution and transmission spectra without any fitting parameters, our formalism models the dielectric environment of the nanostructure with a dielectric constant of $\epsilon_{\rm Au}=10^{5}+10^{5}i$ for the gold metamaterial [29]. The cyclotron resonance of the two-dimensional electron gas is implemented as a gyrotropic medium, where the dielectric tensor of a plasma of charge carriers magnetically biased along the z-direction describes the two-dimensional polarization response of the cyclotron resonance in the plane perpendicular to the magnetic field. In the z direction we only employ the background dielectric constant, as the confinement inhibits a plasma response. Additionally, to reduce numerical complexity of modelling several quantum wells and corresponding barriers between them with finite thickness, we employ an effective medium for the complete quantum well stack with an effective dielectric tensor [43]. In the x-y-direction, we employ periodic boundary conditions to reflect the array character of our structure and solve Maxwell’s equations numerically. The resulting far-field calculations predict experimental results across the entire spectral range with high accuracy.

Our structure is a negative resonator (cut from a gold surface, Fig. 3) [44] of hexagonal shape, fabricated on top of a gallium arsenide (GaAs) substrate (white area in panels a,d), with the cyclotron resonances hosted in doped GaAs quantum well structures (light red region in panels a,d). In order to explore the direct impact of the overlap over the optical spectrum, we consider two types of quantum well designs. In the first layout the 2DEG covers the whole unit cell area (panel a). We refer to this design as unstructured. A second layout is instead realised by in-plane confinement of the 2DEG within a small rectangular patch at the center of the resonator (panel d). We refer to this layout as structured.

The numerical transmissions for the four samples are shown in panels (b,e) as a function of the cyclotron frequency, and in panels (c,f) as a function of the electron density $N_{\textrm{2DEG}}$. The simulation is performed considering an exciting electromagnetic wave which is linearly polarized along the gap (x) direction, and incident perpendicularly to the metamaterial plane. From the transmission spectrum at low electronic density, shown in Fig. 4 (a), we recognize $M=2$ active photon resonances within the given frequency range, whose in-plane field profiles along the gap direction are plotted in Fig. 4 (b,c).

In the structured case the patch acts as a Fabry-Pérot resonator for the quasi-2D plasmonic excitations of the electron gas in the QWs [45]. This leads to a non-vanishing frequency for the fundamental plasmonic mode to which the lower polariton in panel (e) of Fig. 3 would converge for a vanishing cyclotron frequency. We estimated the fundamental plasmon mode frequency using the formula [37]

$\displaystyle\omega_{P}^{0}=\sqrt{\frac{N_{\textrm{2DEG}}e^{2}\pi}{2m^{*}\epsilon_{0}\epsilon_{\mathrm{eff}}W}},$

(24)

with $W$ the patch width and $\epsilon_{\mathrm{eff}}$ the effective permittivity taking into account the screening of the gold resonator by averaging the screened and unscreened portions of the QW area. The resulting value is $\omega_{P}^{0}\approx 0.2\,$THz. Although for the sake of completeness we did use such a value in our simulations for the structured QWs, we notice that for such low frequencies the polaritons have vanishing photonic components and the transmission spectra are not noticeably affected by the exact value of $\omega_{P}^{0}$.

Once we calculated the overlap parameter as in Eq. 11 for the two configurations, we employed our multimode theory to fit simultaneously the resonances for both the spectra of the $\omega_{c}$-sweep and the $N_{\textrm{2DEG}}$-sweep, considering the matter resonance as the magnetoplasmon mode $\tilde{\omega}_{c}$ and treating the normalised mode volumes $\tilde{\mathbb{V}}_{\nu}$ as fitting parameters.

From the discussion in the previous section we expect that passing from the unstructured to the structured sample, as the integration surface is reduced, not only the normalised mode volumes $\tilde{\mathbb{V}}_{\nu,\mu}$ will vary, but also the modes will become less orthogonal, thus increasing the overlap parameter $\eta_{2,1}$. This is indeed the case as can be seen by the calculated values of $\eta_{2,1}=0.32$ for the unstructured sample in Fig. 3 (a) and of $\eta_{2,1}=0.95$ for the structured sample in Fig. 3 (d) derived by Eq. 11.

Comparing the transmission spectra for the different configurations allows us to recognize, albeit in attenuated form, the main differences in the spectral features predicted by the theory (marked on the plots by blue solid lines). At a first glance, we notice that the single polariton anticrossings are well resolved in the unstructured case, as they mainly arise from one-to-one coupling of photonic modes to orthogonal matter excitations. In the structured platform instead, we observe a reduction of the polariton splitting, and the appearance of an S-shaped resonance. The reduction of the polariton splitting is mainly due to the fact that reducing the integration area for the single mode leads to a larger normalised mode volume $\tilde{\mathbb{V}}_{\nu}$, and as such to a smaller coupling strength. On the other hand, the confinement of the 2DEG around the central gap of the resonator increases the overlap between the modes, which becomes close to 1, leading to the appearance of the characteristic S-shaped polariton.

We report in Figs. 5 the in-plane electric field distributions along the gap direction for the coupled modes of the hexagonal resonator platform in the structured configuration. The reported data sets are extracted from the FEM simulation and calculated by our multimode theory, respectively plotted at the right and left sides of panels (a-e).

The simulation field maps are extracted at the resonance and at the value of $\omega_{c}$ marked by the green arrows on the central $\omega_{c}$-sweep transmission plot. The corresponding theoretical electric field profiles are instead obtained by Eq. II as linear combinations of the numerically extracted fields of the uncoupled resonances, shown in Fig. 4, weighed by the photonic coefficients displayed by panels in 5 (f-i). Note that the linear superposition is calculated at a cyclotron frequency within half of the resonance linewidth from the nominal one. This procedure has been necessary because the value of the Hopfield coefficients, and thus the predicted shape of the electromagnetic mode, varies strongly within the linewidth of the broadened modes forming the upper anticrossing.

By observing the field maps we can notice that these refer to three different cases: panels (a,c,e) display map distributions similar to the uncoupled ones, as the weight of one of the two modes is greatly dominant over the other. Panel (d) displays a case in which the two photonic weights are comparable, and the electric field map is noticeably different from either of the bare ones. Finally, in panel (b) our theory predicts the field of the bare photonic mode mainly localised in the central gap, while the simulation shows the electric field diffracting in the far-field of the plasma waves, although remaining confined on the area of the quantum well patch. This effect is related to the one recently investigated in Ref. [19]. Here, the authors point out how the electromagnetic field, confined in the resonator gap, can excite a continuum of propagative high-wavevector plasmonic waves leaking away energy from the polaritonic resonances. In our case the main difference is that the patch acts as a Fabry-Pérot resonator. Even if higher-order discrete modes are quasi-resonant with one more polaritonic branches, the energy of the excited modes remains confined in the patch and has thus only a limited effect on the polaritonic resonances [46]. Our two-mode Hopfield model misses this effect, which could nevertheless be correctly described expanding the basis to include many discrete plasmonic modes of the patch [16] or alternatively using a theory able to deal with continuum spectra [17].

Finally, Fig. 6 highlights the modification of the in-plane electric field driven by the multi-mode hybridization for the specific case of the second polariton mode, from the anticrossing point to its saturation at the higher photonic frequency $\omega_{2}=3.55$ THz. The calculated electric field maps in (b) refer to the cyclotron frequency values marked by the vertical black dashed lines in panel (a) (same as Fig. 5 (g)). We can clearly see how changing the cyclotron frequency varies the electric field map, displacing the minimum of the field across the central gap, a feature suggestive of potential applications in sub-wavelength sensing and optical tweezers.

Our results thus demonstrate that, by optimizing the resonator-2DEG structure, we are able to dynamically modify the sub-wavelength electromagnetic field profile, moving its maxima by varying the applied magnetic field.

In conclusion, we theoretically investigated the multi-mode coupling between the cyclotron resonances of a 2DEG and highly-confined THz-resonator modes. We developed a general theory describing multi-mode coupling taking into account the non-orthogonality of the electromagnetic modes. We highlighted specific spectral features due to the presence of multiple photonic modes and demonstrated the possibility to tune the level of inter-mode coupling by lateral confinement of the 2DEG. Finally using these effects opens up the possibility to dynamically tailor the spatial profile of sub-wavelength electromagnetic modes by varying the applied static magnetic field. This approach can potentially be used to realize sub-wavelength optical tweezers to trap and move nanoparticles over sub-micron distances.

The theoretical results encourage us to explore novel experimental methods and setups allowing to observe the predicted modification of the electric field profiles, driven by the coupling. Moreover, we aim to investigate further different resonator configurations, in order to maximise the effects of the multi-mode hybridization, heading towards novel quantum technological applications, based on a controllable and potentially dynamical tuning of high confined electromagnetic fields.

S.D.L. is a Royal Society Research Fellow. S.D.L. and E.C. acknowledge funding by the Leverhulme Trust (Philip Leverhulme Prize) and the Royal Society (grant RGF\EA\181001). C. L. and J. M. acknowledge funding by the DFG through Grant No. LA 3307/1-2.