Extreme renormalisations of dimer eigenmodes by strong light-matter coupling

We explore by theoretical means an extreme renormalisation of the eigenmodes of a dimer of dipolar meta-atoms due to strong light-matter interactions. Firstly, by tuning the height of an enclosing photonic cavity, we can lower the energy level of the symmetric `bright' mode underneath that of the anti-symmetric `dark' mode. This is possible due to the polaritonic nature of the symmetric mode, that shares simultaneously its excitation with the cavity and the dimer. For a heterogeneous dimer, we show that the polariton modes can be smoothly tuned from symmetric to anti-symmetric, resulting in a variable mode localisation from extended throughout the cavity to concentrated around the vicinity of the dimer. In addition, we reveal a critical point where one of the meta-atoms becomes `shrouded', with no response to a driving electric field, and thus the field re-radiated by the dimer is only that of the other meta-atom. We provide an exact analytical description of the system from first principles, as well as full-wave electromagnetic simulations that show a strong quantitative agreement with the analytical model. Our description is relevant for any physical dimer where dipolar interactions are the dominant mechanism.


Introduction
The electromagnetic environment can alter the physical and chemical properties of atoms, molecules and other emitters. The ability to tailor the states of emitters through strong light-matter interactions could lead to applications in fields as diverse as chemistry, sensing and photonics [1]. Indeed, there have already been concrete demonstrations such as reduced energy losses in photovoltaics [2], increased conductivity of organic semiconductors [3] and enhanced non-radiative energy transfer between spatially separated molecules [4]. It has even led to the development of new research fields such as 'polaritonic chemistry' [5][6][7][8][9][10]. This paradigm can also be exploited in metamaterials. Traditionally, the optical properties of metamaterials are realised by virtue of designing the underlying geometry of optical components [11]. However, instead of ingenious geometrical designs, it has been demonstrated that one can also induce non-trivial, even topological, changes within an optical medium simply by tuning the light-matter interaction strength [12,13].
In this article, we theoretically demonstrate extreme renormalisations of the eigenmodes of a dimer of dipolar meta-atoms by strong light-matter interactions when placed inside a cavity. We consider dipolar meta-atoms as a prototypical system with the purpose to keep our study as general as possible. Such dipolar meta-atoms can accurately model many artificial systems, such as plasmonic nanoparticles [14,15], microwave helical resonators [13], magnonic microspheres [16], as well as Rydberg [17,18] and cold atoms [19]. They are the simplest building block of more complex systems, yet their interaction with light already presents some fascinating properties. We gain an understanding of such a system with a three-mode analytical model of two point dipoles interacting with the fundamental photonic mode of a cubic metallic cavity.
Our general analytical model is enriched by a deeper numerical analysis of a specific physical system: electromagnetic simulations of a dimer of metallic nanospheroids. The numerics also confirm a strong quantitative agreement between the analytical model and full-wave electromagnetic simulations. We choose to focus on nanoparticle dimers as they have already enjoyed a wealth of investigation. The nanoparticle dimer response hybridises into a 'bright' bonding and 'dark' anti-bonding mode [20], whose frequencies depend on the interparticle seperation [21,22], with a polarisation dependence due to the anisotropic configuration [23]. Interesting effects have been predicted and observed, such as interparticle seperation dependent oscillations in the spectral width [24] and between super-and sub-radiant damping [25], Fano profiles in the near-field coupling [26], and non-linear four-wave mixing [27]. In addition, quantum treatments have studied effects such as tunnelling and screening [28] for closely spaced dimers, as well as the Landau damping that becomes important for very small nanoparticles [29]. We also note that the effect of the photonic environment on the energy splitting between the dark and bright mode has demonstrated to be measurable with current technologies, despite the effect on single particle resonances being very small [30].
We show that the energy of the 'bright' symmetric mode of the dimer can be shifted below that of the 'dark' anti-symmetric mode, which could have profound consequences for the optical properties of the system, as it removes a non-radiative decay channel of the bright mode. Such a mechanism could be used, for example, to improve the photoluminescence of photonic devices [31]. There is another interesting effect when we consider a heterogeneous dimer, comprised of two meta-atoms with inequivalent resonant frequencies. In this case, there is an anti-crossing not only between the cavity mode and the symmetric mode, but also a second anti-crossing between the symmetric mode and the anti-symmetric mode. This results in two polaritonic modes which can be smoothly transformed from symmetric to anti-symmetric. This also manifests in a tunable mode width, from extended throughout the cavity to concentrated in the immediate vicinity of the dimer. Such a property could be used to transduce bulk changes into localised resonance shifts, and exploited in sensor applications. Moreover, there is a critical 'shrouded point' where the polariton mode is neither symmetric nor anti-symmetric. Rather, the field re-radiated by the dimer is near-identical to that of a single meta-atom, as if the other meta-atom was not there . This evokes ideas of metamaterial cloaking such as those investigated by Pendry et al. [32], although the mechanism here is quite distinct. Whilst the renormalisation of dimer energy levels due to a photonic environment has previously been studied [30], the novelty in our work is the demonstration that the bright and dark energy levels can be completely inverted in the strong coupling regime, as well as the existence of a novel 'shrouded point'.
The rest of this article is organised as follows. Section 2 outlines the key concepts, expected qualitative results, and introduces the analytical model. Section 3 introduces the system investigated numerically. A dimer of metallic nanospheroids is used as a case study for deeper analysis, whereas the analytical results can be generalised to many other systems. Section 4 presents the results pertaining to a homogeneous dimer such as the inversion of the energies of the symmetric and anti-symmetric eigenmodes; whereas section 5 discusses the heterogeneous dimer properties, such as polariton modes with a tunable weighting of symmetric and anti-symmetric components. We provide an overview and outlook in section 6. The Supplementary Material [33] contains a detailed derivation of our Hamiltonian The meta-atoms have resonant frequencies ω n , and a Coulomb interaction strength Ω. The cavity mode (ω c ) couples to the meta-atoms with a strength ξ n . (c) Schematic of the results. Top panel: A homogeneous dimer has a symmetric (ω + ) and an anti-symmetric (ω − ) mode. Hybridisation with the photonic mode leads to a symmetric upper (ω U ) and lower (ω L ) polariton mode, whilst the 'dark' anti-symmetric mode remains unchanged. In this way, the ordering of the symmetric and anti-symmetric modes can be inverted. Bottom panel: The photonic mode hybridises with both modes of a heterogeneous dimer (ω ± ). An additional anti-crossing results in three polaritonic modes (ω pol m ) that can be smoothly tuned from symmetric to anti-symmetric.
from first principles, further discussions on the validity of our assumptions, and examines the effects of gain and loss.

Analytical model
For simplicity, we model the photonic cavity as a cubic metal box with height L, see Figure 1(a). We consider the meta-atoms as two point dipoles, separated by a distance d, centred in the cavity at r n = (L/2)(1, 1, 1+(−1) n d/L), where n ∈ {1, 2} labels the meta-atoms. The dipoles are associated with bosonic excitations with a resonance frequency ω n , and corresponding creation operators b † n . We consider small cavities such that there is no spectral overlap between the resonances of different electromagnetic eigenmodes, with the spacing between these resonances much larger than the splitting between the dimer eigenfrequencies [33]. In this way, we can consider the cavity as single mode and tune the cavity height such that the fundamental cavity mode, with frequency ω c = πc √ 2/L (with c the speed of light in vacuum) and creation operator c † , is in near-resonance with the dipole resonant frequencies, see Figure  1(b). In the rotating-wave approximation [34], and in the basis Ψ = (b 1 , b 2 , c), the matrix Hamiltonian reads Here, Ω = √ Ω 1 Ω 2 parameterises the quasistatic Coulomb interactions, where Ω n = (ω n /2)(a n /d) 3 , with a n a length scale that characterises the strength of dipolar excitations. A full derivation of the Hamiltonian from first principles is laid out in the Supplementary Material [33]. The coupling constants ξ n between the cavity mode and each meta-atom (labelled by n) can be analytically expressed as [33] ξ n = 8π Ω n ω c We note that if the cavity height becomes too small than the dipole approximation breaks down. For closely space particles we thus require L d. Under this assumption, ξ n is a monotonically decreasing function of L. For the parameters we will consider in this work, ξ n ∝ 1/L to a good degree of accuracy. We must also be careful not to extend this model to arbitrarily large cavity heights where the single cavity mode approximation fails [33]. Gain and loss can be easily incorporated with e.g. a Lindblad master equation approach [33]. In this case, the (Hermitian) eigenvalue spectrum still dictates the absorption peaks in the spectra, but the intensity and widths are determined by the details of driving and dissipation; and of course one must ensure that dissipation does not induce a transition from strong to weak coupling.

Electromagnetic simulations
The analytical model is augmented by focusing on a specific physical system of metallic nanoparticles, analysed using full-wave finite element simulations (implemented with JCMsuite software). We have checked that strong light-matter coupling, and the results discussed herein, can be achieved by embedding the nanoparticles in a metallic cavity made of gold (Au) Fabry-Pérot mirrors with a SiO 2 spacer layer [35]. However, for the sake of simplicity and easy comparison with the analytical model, we model the cavity with perfect electric conductor boundary conditions. For dipolar meta-atoms, we consider metallic ellipsoidal particles (with semi-axes R x = 15nm and R y = R z = 10nm), where the permittivity is given by a simple Drude model with parameters ω p and γ chosen such that permittivity at a wavelength of 500nm is −3.6 + 0.01i. Center-to-center distance between the particles is 35nm. The system is excited with two dipole sources with equal dipole moments, placed at (±77.5nm, 0, 0). These dipoles can be either parallel or antiparallel to each other, thus allowing to excite either the bright or the dark mode of the system. We then perform frequency sweeps for each cavity height, and by using the simulated fields, calculate the total power absorbed by the particles. This allows us to identify the system resonances and plot field distributions of various modes.

Homogeneous dimer
Let us first consider an isolated homogeneous dimer, made from meta-atoms with the same resonance frequencies ω n ≡ ω 0 and length scales a n ≡ a, and thus the same light-matter couplings ξ n ≡ ξ. In this case, the dipolar Hamiltonian H dp in equation 2 is diagonalised by the operators with eigenfrequencies ω ± = ω 0 ± Ω, corresponding to a higher energy symmetric mode with aligned dipoles (↑↑) and a lower energy anti-symmetric mode with anti-aligned dipoles (↑↓). Often, these modes are referred to as the 'bright' and 'dark' modes due to their finite and vanishing total dipole moment respectively in the far-field.
Indeed, the anti-symmetric 'dark' (ω − ) mode has no coupling to the cavity mode and remains an eigenmode of the full light-matter Hamiltonian H Ψ . This can be intuitively understood with symmetry arguments, by considering that an excitation of the symmetric cavity mode is incapable of subsequently exciting an anti-symmetric dimer mode. On the other hand, the hybridisation between the cavity and symmetric dipolar (ω + ) mode results in a Rabi-splitting of magnitude Ω R = (ω + − ω c ) 2 /4 + 2ξ 2 and a (symmetric) upper and lower polariton eigenmode with frequencies ω U,L = (1/2)(ω + + ω c ) ± Ω R . Note that the conventional strong-coupling criteria Ω R > (γ 1 + γ 2 )/2 is satisfied for the damping rates chosen in the numerical model . The Rabi splitting causes the symmetric lower polariton mode to decrease lower in energy as we increase the cavity height, see Figure 2(a). In fact there is a critical cavity height L c where the symmetric (lower polariton) and anti-symmetric (dark) mode become degenerate. This occurs when the Rabi frequency equals a certain detuning ∆ = (1/2)(ω + +ω c )−ω − , namely the detuning between the dark mode and the average frequency of the bright and photon mode. For larger cavity heights, L > L c , the energies of the symmetric and anti-symmetric mode become inverted.  Note that there are no fitting parameters here. Rather we extract ω 0 and ω + from the resonance peaks of simulations for a single particle and dimer in free space respectively, which via the expression for Ω also tells us the value of a. The colour scale corresponds to the absorbed power in the metallic nanospheroids as a function of the driving frequency ω dr of symmetrically driven point dipoles which are used to excite the system. We see that the peaks of these data follow the same trend as the (lower polariton) symmetric mode ω L . There is no absorption peak at the energy of the dark mode as this anti-symmetric mode is not excited by the symmetric driving. Cuts at fixed cavity height are displayed in Figure 2(d) showing in more detail how the intensity spectrum shifts with cavity height. Such data is likely the most convenient to compare to possible experimental implementations. In addition, Figure 2(c) shows the electric field distribution Re(E y ) for a system driven resonantly in a symmetric (bottom panel) and anti-symmetric (top panel) configuration. We clearly see the hybrid light-matter nature of the symmetric mode, with the field intensity pattern extending across the whole cavity; whereas the anti-symmetric mode has a field pattern that is near-identical to its free-space counterpart, and concentrated in the vicinity of the dimer.

Heterogeneous dimer
There is no completely 'dark' mode of a heterogeneous dimer, where ω 1 = ω 2 in equation 2. In this case, the cavity mode hybridises with both dimer eigenmodes, leading to three polariton bands with an additional anti-crossing and some new physics to consider. To start our investigation, it is enlightening to rotate the full Hamiltonian H Ψ into the basis of symmetric and anti-symmetric operators. These operators are the very same which diagonalised the homogeneous dipolar Hamiltonian. Thus in the basis α = (α + , α − , c) the Hamiltonian reads where we now define ω 0 = (1/2)(ω 1 + ω 2 ) as the average of the resonance frequencies of the two metaatoms, and δ = (1/2)(ω 1 − ω 2 ) is the deviation from the average. We see that the light-matter coupling is large for the symmetric mode and small for the anti-symmetric mode. In this sense the symmetric mode is 'brighter' whilst the other is 'darker'.
As with the homogeneous dimer, the Rabi-splitting between the photonic mode and symmetric dimer mode causes the (originally) symmetric polariton mode to decrease in energy as the cavity height increases, whilst the (originally) anti-symmetric polariton mode has a small unimportant light-matter coupling. The key qualitative difference is that the frequency difference between the single particle resonances leads to an additional anti-crossing between the symmetric and anti-symmetric dimer mode. This can be clearly seen in the inset of Figure 3(b), which shows the dark anti-symmetric mode and the symmetric lower polariton mode of the homogeneous dimer (dashed lines), as well as the polariton bands that arise from the hybridisation between these two modes for a heterogeneous dimer (solid lines). We note that the 3rd polariton mode (not shown) is near-identical to the upper polariton mode of the homogeneous dimer (blue curve in Figure 2(a)), as a result of the weak coupling to the anti-symmetric dimer mode.
The result of the hybridisation between the anti/symmetric modes is that the polariton eigenmodes are neither fully symmetric nor anti-symmetric, and in-fact can be smoothly tuned from symmetric to anti-symmetric. To quantify how (anti-)symmetric a mode is we introduce the measure where v m = (v m+ , v m− , v mc ) is the eigenvector of the Hamiltonian H α corresponding to the m-th eigenvalue (ordered in increasing energy). A value of S m = ±1 corresponds to a fully symmetric and anti-symmetric mode respectively. Figure 3(a) explicitly shows how S m varies from positive to negative (and vice versa) for the 1st and 2nd polariton modes. Figure 3(b) shows electromagnetic simulations of the response of the system to both a symmetric and an anti-symmetric driving. Here, the nanoparticles have the same dimensions, and we alter the plasma frequency to obtain inequivalent single particle resonances; see equation 4. It would correspond to a dimer made from different materials [36]. Specifically, Figure 3(b) shows the absorption under a symmetric drive (red), overlaid on the absorption under an anti-symmetric drive (black). By plotting the absorbed power in both cases we observe how the anti/symmetric makeup of the polariton eigenmode influences the response of the system depending on the symmetry of the driving. Figure 3(c) further highlights the tunable polariton modes. As the cavity height is increased, the 2nd polariton mode changes from symmetric to anti-symmetric, resulting in an electric field distribution that contracts from extended throughout the cavity to concentrated in the vicinity of the dimer.
The 3rd panel of Figure 3(c) shows the critical cavity height L = 339nm when the left meta-atom becomes 'shrouded'. No electric field is excited on the left nanospheroid in response to the symmetric driving field, and the field re-radiated from the dimer is that only of a single meta-atom. Which metaatom becomes shrouded is determined by that one with a higher (respectively lower) single meta-atom resonance frequency, when driving the system at the resonance of the lower m = 1 (respectively middle m = 2) polariton branch. In the analytical model this corresponds to S = 0, when the polariton eigenmode has an equal weight of both the symmetric and anti-symmetric mode (v 2+ = −v 2− and thus is a linear superposition of b 2 and c). Interestingly this means that the cavity has renormalised the field of a dimer to resemble that of a single meta-atom.

Discussion
We have demonstrated an extreme renormalisation of the energies of an interacting dimer of meta-atoms due to strong light-matter interactions, simply by tuning the height of an enclosing photonic cavity. For a homogeneous dimer this resulted in an inversion of the ordering of the symmetric and anti-symmetric modes; a striking fundamental demonstration that one can impart non-trivial changes to the optical properties of metamaterials, without altering the underlying geometry or symmetry of the constituents. Indeed a similar mechanism has recently been proposed to significantly enhance the photoluminescence efficiency of carbon nanotubes by removing non-radiative decay channels to a dark exciton state [31].
In the case of a heterogeneous dimer we have seen that the polariton modes can be smoothly tuned from symmetric to anti-symmetric, resulting in a tunable mode width that varies from extended throughout the cavity to concentrated around the vicinity of the dimer. There is also a critical cavity height where one of the meta-atoms becomes 'shrouded'. The mode is neither symmetric nor antisymmetric, in fact the electric field re-radiated by the dimer resembles that of a single meta-atom, almost as if one of the meta-atoms was not there.
These phenomena should be realisable in any physical dimer where the interactions between the emitters are predominantly dipolar. In fact, qualitatively similar results should be possible in any system where one can independently tune the cavity-emitter interactions and emitter-emitter interactions. As one example, we provided full electromagnetic simulations of metallic nanoparticles. These results were in strong quantitative agreement with the analytical model derived from first principles, with no fitting parameters, showing that the simple point dipole model captured the key physics.
Our results highlight the incredible tunability offered by strong light-interactions of meta-atoms; which could find use in e.g. plasmonic transducers and sensors [37,38]. Indeed, our scheme offers opportunities for sensing devices that transduce a change in a bulk property affecting the extended cavity mode, to that of a localised resonance shift of the dimer. Moving forward, it could be interesting to investigate quantum effects such as a possible enhancement or modulation of entanglement effects, similarly to that reported for magnons [39]. We outline the derivation of the Hamiltonian in the main text from first principles. During this process we also derive the general Hamiltonian for all cavity modes (including non-resonant terms) which could be useful for future investigations. We show how the inclusion of only the fundamental mode = (0, 1, 1) is valid for the parameters considered in the main text. Finally we extend the model to include image dipoles, and show that the omission of these and non-resonant terms was valid in the main text.

Dipolar Hamiltonian
In this section we detail the derivation of the dipolar Hamiltonian. We consider the quasistatic Coulomb interactions between two dipolar meta-atoms. We consider meta-atoms composed of physical objects whose natural frequency of the charge density oscillations along the x-direction are tuned in resonance with the fundamental cavity mode. Whereas the natural frequencies along the other directions are detuned and not relevant to the description. In this way we need only consider the one transverse polarization along the x-axis. We consider a dimer of meta-atoms labelled by n ∈ {1, 2}, whose positions are r n = (1/2) (L x , L y , L z + (−1) n d). They have masses M n , natural frequencies ω n , and charges Q n = N n e 2 , where e is elementary charge constant. We denote the displacement of the electronic centre of mass by h n , and the conjugate momentum by Π n . Including the quasistatic Coulomb interactions the Hamiltonian of the meta-atoms reads [29] H dp = n Π 2 n 2M n + M n ω 2 where is the dielectric constant of the surrounding medium. Throughout this work we work in units such that = 1 and = 1 but retain the symbolic representation of for completeness. We switch to the second quantisation representation by the substitutions which gives where Ω = √ Ω 1 Ω 2 parameterises the strength of Coulomb interactions. Here, Ω n = (ω n /2)(a n /d) 3 where a n = (Q 2 n /4π M n ω 2 n ) 1/3 is a length scale characterising the strength of the dipolar excitations.

Magnetic vector potential of a box with perfect metallic walls
Before deriving the light-matter Hamiltonian in the subsequent section, first we must deduce the magnetic vector potential in a box with perfect metallic walls. We follow the paper by K. Kakazu and Y.

Deriving A(r, t) in full generality
We consider a box with perfect metallic walls. The cavity has sides of length L x , L y and L z . The boundary conditions are that the tangential component of E and the normal component of B vanish at the cavity walls. Maxwell's equations in free space leads to the wave equation Thus for both the magnetic and electric field, for each component we need to solve an equation of the form Using separation of variables and assuming a harmonic time dependence, the solution is of the form ψ(x, y, z, t) = X(x)Y (y)Z(z)e −iωt . Thus we have X + k 2 X = 0 with general solution and similarly for Y and Z. Using Maxwell's equations, particularly ∇ × E − ∂B ∂t = 0, and the boundary conditions allows us to deduce the specific form of the electric field eigenmodes as v x = N x cos(k x x) sin(k y y) sin(k z z)x, (13) v y = N y sin(k x x) cos(k y y) sin(k z z)ŷ, where N i = 8/(V (1 + δ i )) is the normalisation coefficient, δ i ≡ δ( i ) is the Kronecker delta of i , and k i = π i /L i is the 'wavenumber'. Here, i ∈ {x, y, z} labels the Cartesian coordinates, and we label the photonic modes with the list of integers = { x , y , z }. The normalisation coefficient arises from imposing orthonormality cavity dV v i v i = δ ii δ . Note also that not all the combinations of (l x , l y , l z ) contribute non-zero eigenvectors (e.g. (0, 0, 1) and its permutations). These can be easily rejected.
Thus we can construct a general solution of the electric field from a linear combination of these field eigenmodes where D i are the expansion coefficients, and the phase factor is introduced by convention. The transversality condition (derived from ∇ · E = 0) is This shows us that we are over specifying the system, and in fact only need two orthogonal unit vectors to resolve the components in the expansion of E. Therefore let us perform a rotation to a coordinate frame whereê 3 =k and the other two unit vectors,ê 1 andê 2 , are orthogonal to k . First write E as where D and V are vectors formed from the components D i and v i in the obvious way (note that technically V is a vector of vectors). Then we introduce the orthogonal 'change of basis' matrix Q which (because orthogonal matrices preserve the dot product when acting on both of the vectors) allows us to write E as where U = Q V and C = Q D ; their components are u λ and C λ respectively, where λ ∈ {1, 2, 3} (we will shorten the set of λ in a moment). Going back to the summation notation we have Finally, we choose the orthogonal matrix where θ = arctan k z , k 2 x + k 2 y and φ = arctan(k x , k y ) are the polar and azimuthal angles of k . The transversality requirement is now easily satisfied by setting C 3 = 0. Thus from now on summations over the different modes can be restricted to the set λ ∈ {1, 2}. Finally we quantise the system by promoting the expansion coefficients to annihilation operators where ω = c|k |, and c λ is the annihilation operator obeying dc λ /dt = −iω c λ . The normalisation coefficient is chosen so that the Hamiltonian has the simple form where the global energy shift of λ (1/2) has been dropped. Getting back to the task in hand, we now write the electric field with operators as From this we can automatically use Maxwell's equation, ∂B/∂t = −∇ × E, to write the magnetic field, As this expression already contains a curl, we readily obtain the expression for the magnetic vector potential as In the expression for A we have implicitly chosen the Coulomb gauge by omitting the gradient of an arbitrary function. The mode functions are given explicitly by where the prefactors due to the two photon polarisations are

The component of A aligned with the dipole polarisation at the metaatom sites
We now derive the expression for the component of the magnetic vector potential aligned with the dipole polarisation (x-axis), exactly at the location of the meta-atoms r n = (1/2)(L x , L y , L z + τ n d), where τ 1,2 = ±1 labels the meta-atom. By simple substitution, we can first write the x-component of the magnetic vector potential at the nanoparticle positions as Several of the terms in this sequence are zero, and thus it would be useful to relabel the system to remove this redundancy. This has the practical advantage that we do not have to keep track of which modes actually have any coupling to the meta-atoms and which do not. Thus we exploit the following identities n∈N cos πn 2 f (n) = n∈N (−1) n f (2n), n∈N sin πn 2 f (n) = n∈N (−1) n f (2n + 1), sin πn 2 (1 ± x) = (±1) n+1 sin πn 2 (1 + x) .
These identities allows us to write the magnetic vector potential as where the sum over is restricted to the set P where x is even, y is odd, and z is any natural number. Here, the parity of the mode with polarisation λ at nanoparticle n is encoded in the quantity P n = (−1) x+ y (τ n ) z +1 , where obviously P n ∈ {1, −1}.

Light matter interaction Hamiltonian
Now that we have derived the magnetic vector potential, we can move on to obtaining the light-matter interaction Hamiltonian. Neglecting diamagnetic terms which only introduce a negligible self-energy renormalisation of the photon frequencies, the light-matter coupling in the minimal coupling formalism is which can be expressed in terms of the bosonic operators as The light-matter coupling strength is parameterised by 3.1 Including only the fundamental cavity mode We tune our system such that only the fundamental cavity mode = {0, 1, 1} is in near-resonance with the meta-atom resonance frequency and has a significant light-matter coupling. We will justify this in later sections. In this case the interaction Hamiltonian reads where Here, c † creates a photon in the mode = {0, 1, 1} with polarisation λ = 2 and frequency ω c ≡ ω {0,1,1} (the other polarisation has no component in the x-direction). The photonic Hamiltonian is trivially

Significance of the other cavity modes
To get a feel for the action and significance of the different modes, let us assume that we have (and that it is possible to have) tuned the cavity such that we can neglect all other modes apart from (an arbitrary) and λ, and consider a homogenous dimer (ω 0 ≡ ω 1 = ω 2 ). In this case the Hamiltonian in the rotating wave approximation, in the basis where we have defined P 1,2 ≡ P 1,2 , as well as ξ ≡ ξ 1 λ = ξ 2 λ and ω c ≡ ω . Let us write this Hamiltonian in the basis γ † = (γ † B , γ † D , c † ) where the 'bright' and 'dark' operators are γ B,D = (1/ where we have defined ω B,D = ω 0 ± P 1 P 2 Ω. Interestingly, we see that the mode which couples to an individual cavity mode is not necessarily the higher energy symmetric dipolar mode; the one we typically consider as the 'bright' mode in free-space. For greater clarity we can write the new true bright and dark eigenmodes (in this context) as  We characterise the extent to which a cavity mode renormalises the bright mode of a heterogeneous dimer, when considering just that cavity mode in the model. We plot the percentage change of the polariton mode ω pol away from the uncoupled bright-mode ω B for x ∈ {0, 2, 4, 6}, y ∈ {1, 3, 5, 7}, z ∈ {0, 1, 2, 3}, and both polarisations. We show the mode with = {0, 1, 1} and λ = 2 in red, which is exactly the mode considered in the main text. All other parameters are the same as in the main text. The cyan region corresponds to the range of L considered in the main text.
We can also separate the Hamiltonian into two block diagonal pieces; the dark mode ω D , and the Hamiltonian of the photon and bright modes. In the basis (γ † B , c † ) this reads The polariton eigenmodes are where we have defined the average frequency ω Bc = (1/2)(ω B +ω c ) and detuning δ Bc = (1/2)(ω B −ω c ). In figure 4 we plot the percentage difference between this polariton mode and the uncoupled bright-mode, for several cavity modes. This gives us a measure of how significant each cavity mode is in the full general light-matter Hamiltonian of equation 34. We see that in the range of cavity dimensions considered in the main text the mode = {0, 1, 1} is the only mode which causes a significant renormalisation of the bright mode.

Non-resonant terms and image dipoles
In the main paper we have neglected non-resonant terms (such as b † 1 b † 1 ) as their contributions are negligible for Ω, ξ ω 0 . In both the main paper and this Supplementary Material we have also neglected image dipoles. We will reintroduce both these contributions in this sections and show that they alter the spectra in insignificant ways.
The boundary conditions of the metallic box renormalise the purely dipolar electric field of the point dipole. This can be represented by image dipoles, leading to a image dipole Hamiltonian where I (s) = =0 (−1) y + z C x L x , y L y , z L z + dmod 2 ( z + s) and C(x, y, z) = y 2 + z 2 − 2x 2 [x 2 + y 2 + z 2 ] 5/2 .
Here, the summation is over all α except when they are all zero. Thus the full non-resonant Hamiltonian including image dipoles is H = H dp + H im + H ph + H int , where the component Hamiltonians are given in equations 9, 42, 23 and 33 respectively. In figure 5 we plot the dispersion of this Hamiltonian (for the mode = (0, 1, 1)) including non-resonant and image dipoles (red lines) as well as that of the Hamiltonian from the main text which neglects these terms (dashed black lines). We see that the dispersions are near-identical.

Effects of gain and loss
To investigate dynamics, as well as the effects of gain and loss, we can take the usual approach of forming the Lindblad master equation for the density matriẋ where X m ∈ (b 1 , b 2 , b c ) and γ m are the damping rates. We consider the full non-resonant Hamiltonian (neglecting image dipoles) where the first three terms are given by equations 9, 37 and 35 respectively; and where we consider the symmetric (and anti-symmetric) driving f 1 (t) = ±f 2 (t) = exp(iω D t). We introduce the (dimensionless) dipole moment p n = b n + b † n and conjugate momentum π n = i b n − b † n ; and correspondingly define p c = c + c † and π c = i c − c † . The equations of motion for these expectation values can be obtained using Ȯ = Tr(ρO), and exploiting the cyclic properties of the trace. They arė p n = −ω 0 π n − 2ξp c − γ n p n , (48) p c = −ω c π c − γ c p c , π n = ω 0 p n + 2Ωp n+1 − γ n π n + 2f n , π c = ω c p c − 2ξ(π 1 + π 2 ) − γ c π c + 2f c .
Using a complex field amplitude representation, we can solve for the steady-state solutions h n = A n exp(iω D t), and hence obtain the total dipole strength, |A 1 | + |A 2 |, as a function of cavity height and driving frequency ω D . As a first approximation we consider this dipole strength proportional to the absorbed field at the meta-atoms, which we plot in Figure 6.