Strong coupling regime and hybrid quasinormal modes of a single plasmonic resonator coupled to a TMDC monolayer

One common goal of nano-optics is to create photonic cavity modes that can significantly enhance light-matter interactions, leading to new capabilities in sensing and quantum optics. While much success has been achieved with semiconductor cavity structures [1], recently substantial attention has been devoted to metallic nanoparticles (MNPs), inspired by their ability to trap light in sub-wavelength spatial scales and allow a broadband enhanced coupling. This can help enable new regimes in quantum plasmonics [2, 3], including the experimental demonstrating of strongly coupled single molecules and MNPs at room temperature [4].

Recently, there has been much interest to increase the interaction strengths even further by coupling to 2D semiconductors, such as transition metal dichalcogenides (TMDCs). Monolayers of TMDCs are direct gap semiconductors with strong light matter interactions which makes them promising for optoelectronic applications [5, 6]. The reduced dimensionality of these atomically thin materials leads to a boost of the Coulomb interaction which is responsible for the formation of tightly bound electron hole pairs (excitons) [7]. The large excitonic binding energies make TMDCs a promising platform to study exciton physics [8, 9]. The TMDC excitons form separate, spectrally well-isolated resonances below the band gap and can couple radiatively to MNPs if their dielectric surrounding is properly designed. This allows one to also achieve signatures of strong coupling from near resonant TMDC-MNP hybrid systems which gives rise to large Rabi splittings.

To date, a variety of MNPs have been used to observe signatures of strong coupling between MNPs and TMDCs, including single MNP nanorods [10, 11, 12], bi-pyramids [13], disks [14], spheroids [15], dimers [16], and arrays of MNPs [17, 18, 19]. Measurements are often performed using reflection/transmission from a low-powered laser or dark-field scattering [10], and are typically performed at room temperature. Impressive spectral Rabi splittings of 500 meV have been achieved using a microcavity coupled to a TMDC sheet with arrays of gold disks [19] (although the splitting arises primarily from the classical mode splitting between the cavity and MNP array), while the average Rabi splittings in recent works range from approximately 50-150 meV. Most if not all of these works are in the regime of so-called “normal mode splitting,” and thus the physics of the splitting can be well described classically [20]. Despite numerous experimental observations, there are no rigorous 3D models in the literature to properly assess the strong coupling regime in these material systems, which is significantly more challenging than typical planar structures, e.g., with quantum wells [21], since the MNPs break translational invariance. The terminology of “normal mode splitting” is also not appropriate, as the MNPs on their own do not support such modes, requiring a more rigorous mode theory from the beginning.

Besides normal mode splitting with planar-like cavity systems, strong coupling phenomena of single quantum dots (which naturally break translational invariance) have been reported in a number of semiconductor microcavity systems [22, 23, 24], where the emitter-cavity coupling rate is defined through $g=(\pi e^{2}f)/(4\pi\epsilon_{r}\epsilon_{0}m_{0}V_{m})^{1/2}$, with $f$ the oscillator strength, $m_{0}$ the free electron mass, and $V_{c}$ the cavity mode volume. In these systems, experimental Rabi splittings up to 140 $\mu$eV have been reported, and the description of $g$ as a cavity-emitter coupling rate is more appropriate since it is describing a single quantum emitter and the light field within a dipole approximation; the single mode approximation is also excellent here since the cavity system is usually very low loss, yielding substantial quality factors ($Q>1000$). The interpretation of these emitter-cavity systems in the strong coupling regime is clear since the emitter-cavity coupling rate is larger than any dissipation rates. However, since the upper limit to $g$ in these systems is bound by the diffraction limit, most, if not all, of the single quantum dot experiments in this field have to work at very low temperatures.

Despite the similarities between point dipole emitters in cavities and TMDC-plasmon systems, to our knowledge, no one has reported a first-principles electromagnetic calculation of the underlying hybrid modes for strongly dissipating metal structures using material models with fully microscopically determined parameters for the TMDC. The understanding of strong coupling is also essential to develop mode quantization theories, which can then allow true signatures of quantum optical effects in these material systems.

In this work, we introduce a rigorous electromagnetic theory and model to describe the regime of strong coupling in TMDC-MNP systems, including a full analysis of the underlying hybrid modes. In particular, we develop a quasinormal mode (QNM) theory of such hybrid systems, enabling a full 3D description of the modes, light-matter interactions, and the strong coupling behavior of these modes. Quasinormal modes are known to be highly accurate for describing such hybrid material systems, and have recently shown much success in modelling a wide range of resonators in nano-optics and plasmonics [25, 26, 27, 28, 29]. Here we explore the QNMs of a MoSe₂ monolayer TMDC encapsulated in hBN coupled to a gold nanorod MNP. By tuning the dimensions of the MNP to get the desired resonance frequency, we show strong coupling between the second lowest mode of the MNP, which is the next mode in frequency space from the fundamental mode, and the TMDC for various gap sizes between the two structures.

We then compare the QNM analysis to the full-dipole numerical calculations (with no approximations), and show excellent agreement to justify the validity of the QNM analysis over a wide range of frequencies and spatial points near the resonators. We also compare the optical response for selected temperatures of 4, 77, and 300 K with microscopically calculated material parameters [30, 31] and show the anti-crossing behavior when detuning the resonances by varying the size of the MNP, which is a clear signature of strong-coupling. We then compare the splitting observed from the spectral Purcell response and the splitting between the complex poles of the hybrid QNMs, with no fitting parameters. Lastly, we also study how the transmission spectra from an embedded dipole differs from the near field Purcell factor, which show complimentary but different signatures of the strong coupling regime. The dipole transmission requires the full two-space-point Green function and accounts for spatial quenching between the dipole and the detection point, while the Purcell factor is a measure of the projected local density of states at the dipole location, which is also contained within the Green function with equal space points. Conveniently, we construct an analytical solution to the full two-space-point Green function through an analytical expansion of the hybrid QNMs, which critically contains the QNM phase, and can be applied to a wide range of problems in nanophotonics.

We first summarize some common models and results reported in the literature, and discuss the various formalisms currently being used to explain the TMDC-MNP strong coupling behavior. For this specific work we will refer to “strong coupling” in the semiclassical sense of QNM spectral splitting when the two hybrid modes are suitably coupled together. As is well known, strong coupling in the linear regime can also be explained from the classical perspective of linear dispersion theory [32], and this is essentially the regime we also explore below, though we adopt a rigorous QNM approach.

We begin by introducing the usual textbook model of coupled harmonic oscillators, also commonly used to describe strong coupling and (sic) “normal mode splitting” in TMDC-MNP systems [20, 33]:

where $\mathbf{P}$ is the polarizability of the TMDC exciton due to the dominant exciton of interest, $\mathbf{E}$ is the electric field of the interacting cavity mode, $g$ is the exciton-cavity coupling rate, $\gamma_{x,c}$ is the phenomenological decay rate of the TMDC exciton or cavity mode, and $\omega_{ex,c}$ is the resonant frequency of the TMDC or cavity. The fields $\mathbf{P}$ and $\mathbf{E}$ are both 3D vectors in general, but the excitons live primarily in-plane to the TMDC, so these vectors are effectively 2D, although out-of-plane excitons with comparably weak radiative interaction can exist using stacked monolayer designs [15].

The strong coupling regime here characterizes the coupled mode splitting of some spectral signature such as the emission spectra, which can be seen as representing classical oscillations of particles and fields in the time domain (if adopting a two level atom model for one of the oscillators). This also allows one to construct a quantized field theory in order to invoke the quantum language of strong coupling (i.e., vacuum Rabi oscillations) to connect to quantum optics formalisms like the Jaynes-Cummings (JC) model. However, in contrast to spatially extended systems such as MNP-TMDC hybrids, the JC model (or quantum Rabi model, if a rotating wave approximation is not made [34]) either assumes point coupling (such as in atom cavity-QED), or constant coupling as a funbction of space, such as with quantum wells in DBR (distributed Bragg reflector) mirrors.

Phenomenologically, as a first approximation, the eigenfrequencies of this simplified hybrid system can be found by assuming a harmonic time-dependence for $\mathbf{P}$, $\mathbf{E}$ $\propto\exp{(\pm i\omega t})$. When the system is on resonance, such that $\omega_{0}\equiv\omega_{c}=\omega_{x}$, and assuming that $\gamma_{x},\leavevmode\nobreak\ \gamma_{c}\ll g,\leavevmode\nobreak\ \omega_{c},\leavevmode\nobreak\ \omega_{x}$, the complex eigenfrequencies are readily obtained as:

where $\hbar\Omega_{R}\equiv\hbar\sqrt{4g^{2}-(\gamma_{ex}-\gamma_{c})^{2}}$ is the spectral mode splitting. Thus, the coupling constant can be determined by the observed splitting, using

which can be thought of as a JC type cavity-emitter coupling rate, though this is usually defined for a single two-level atom or point dipole emitter [35].

The usual criterion for observing the strong coupling regime, for vacuum dynamics or linear excitation, is often given by the condition: $\Omega_{R}=2g>{|\gamma_{c}-\gamma_{ex}|}$. However, for large damping rates of the cavity mode, a more accurate measure is given by [20, 11, 12]

but even this criterion is only a crude approximation for spatially extended TMDC-MNP systems. Clearly one needs to understand the underlying properties of the modes before defining what is meant here by strong coupling, especially for the development of a quantum model. In particular, at best, $g$ is an effective coupling constant, but it does not tell us anything about the spatially dependent coupling.

For the strong coupling regime, in addition to satisfying the approximate relation, $\Omega_{R}>\gamma_{x}+\gamma_{c}$, the other usual assumption is that $g/\omega_{0}<0.1$, otherwise the rotating-wave approximations used in the derivation of the JC models does not work. Above this threshold, we enter the so-called “ultrastrong” coupling (USC) regime [34, 36]. However, at the semi-classical level, such effects are in fact easy to account for, namely we do not have to invoke a rotating wave approximation in the theory.

To date, almost all works to describe strong coupling between TMDCs and MNPs in the literature use this simple coupled oscillator model for explaining experimental data, with reasonably good success in fitting the locations of the eigenfrequencies [37, 38, 11, 12, 13, 39, 40, 41]. For example, the coupled oscillator model was used to describe the strong coupling behavior between WSe₂ and a gold bipyramid MNP [13]; the work explored the coupling as a function of number of TMDC layers, which resulted in a saturation of coupling as the number of layers increased. Simulation fits with the simple model and experimental data are good, but they fail to account for the actual effective mode volume, $V_{\rm eff}$, of the MNP being a function of space, and the meaning of $V_{\rm eff}$ is also ambiguous in general, since it is defined from an effective theoretical coupling constant,

where $d$ is the exciton dipole moment, $N$ is the number of excitons, $c$ is the speed of light in vacuum, and $\lambda$ is the wavelength; here $V_{c}$ is clearly a spatially-averaged effective mode volume, but in reality it changes as a function of position as does $g$. This simple theory likely relates back to older works that exploit a simple coupled mode theory with translational invariance, and also without any dissipation in general, since the theory exploits the properties of “normal modes”, which are not the correct cavity modes of open resonators [42, 43]. In addition, it is not clear how to choose the ratio of number of excitons $N$ per mode volume $V_{C}$ Reference [44] also uses a bipyramid nanoparticle coupled to MoSe₂, which employs the same assumption of a fixed and heuristic effective mode volume from a two oscillator model. Moreover, most reports in this field [12] do not properly calculate the mode volume of a lossy mode system; e.g., in this case, it was assumed that $V_{c}\propto(\epsilon^{R}+2\omega_{c}\epsilon^{I}/\gamma_{c})|E_{c}|^{2}$, and $N$ is also not well defined. This is because the localized plasmons excited by the 3D nanoparticle break translational invariance, even if the TMDC-MNP hybrid system is excited with a plane wave.

A recent review paper [45] gives a detailed critique on the use and warnings of using $N$ and $V_{c}$ in many nanophotonic systems, highlighting the limitations of Eq. (5) and the problems that can arise. For example, with collective exciton systems, usually a total coupling constant is determined as $g=\sqrt{\sum_{i}g(\mathbf{r}_{i})^{2}}$, where $i$ is the $i^{\rm th}$ exciton, and the exciton Bohr radius is assumed to be constant over the entire TMDC sheet. The authors highlighted that for a TMDC coupled to plasmonic system, this is obviously not applicable, and a more appropriate calculation would include a more detailed theory of excitons as spatially extended composite particles.

To help address this problem in terms of the underlying open cavity modes, plasmonic cavity systems can be rigorously described in terms of QNMs, which are the complex eigenmode solutions for resonators with open boundary conditions. The use of QNMs allows one to obtain the correct position-dependent mode volume (which simply characterizes the normalized mode strength squared) and coupling constant as functions of space which directly involves the excitonic wave function. The same theory also gives the photonic Lamb shifts and collective effects if required, e.g., coupling to multiple quantum dots at different spatial locations, and important propagation effects. These QNM can also be fully quantized [46] allowing for the exploration of new physics beyond the simple JC models [47], and accounting for important QNM phase effects. Thus, it is highly desirable to connect to the underlying QNMs of the MNP, before and after coupling to the TMDC system. In the presence of TMDC coupling, the “hybrid” QNMs fully characterize the classical electromagnetic coupling. Notably, these hybrid QNMs are valid modes regardless of whether we work in the weak or strong coupling regime, so they offer a unified description. Just like two coupled atoms, which can form super-radiant and sub-radiant states, the coupling between photonic resonators can also form analogues of these hybrid states, and are the natural dressed-states of the system.

For the MNP, we use a simple cylindrical gold nanorod with rounded ends (see Fig. 2) where the radius (R) and length (L) can be tuned to achieve a spectral resonance that overlaps with the TMDC resonance. The second lowest mode (with respect to frequency) of the MNP is chosen to couple to the TMDC as it is typically much larger in quality factor (thus, smaller decay rate), compared to the first order mode, which can aid in obtaining strong coupling behaviour with the TMDC; note that Ref. [11] uses the third lowest mode of a silver nanorod to achieve high-$Q$ low-mode-volume conditions, requiring a much larger MNP. All calculations discussed below are fully 3D in nature without any approximations or fitting parameters, other than the microscopic parameters already stated in the models.

For our models, we will first compute the MNP QNM as well as the hybrid QNMs formed from the coupled TMDC-MNP system, at different gap separations and for different temperatures. The QNMs are spatially dependent and contain important phase effects and frequency information (e.g., from the complex poles). Using the QNMs in a Green function expansion, we can easily compute a range of important light-matter interactions at any spatial position over a wide range of frequencies, and analytically. We also stress again, that these QNMs can also be fully quantized for use in a rigorous open-system quantum optics model [46] if needed. In typical experiments, the absorption, transmission, or reflection from an optical source is usually what can be used to detect this mode hybridization, e.g., in combination with near field optical techniques. One can also add additional emitters such as single quantum dots, and study the modified spontaneous emission (Purcell effect), as well as light propagation and photon exchange effects.

As an application of the theory, we will report on two frequency-dependent functions that are relevant to probing strong coupling. The first is the Purcell factor (PF) as a function of frequency, which is related to the modified spontaneous emission of a point dipole at some position ${\bf r}_{d}$. The Purcell factor depends on the projected local density of states (LDOS), which is completely contained within the imaginary part of the (photonic) Green function, ${\rm Im}\left[{\bf G}({\bf r}_{d},{\bf r}_{d},\omega)\right]$. Since we are interested in the strong coupling between the TMDC and the MNP, we consider an immutable dipole (weakly coupled and in a linear response, though this is not a model restriction). Second, we consider the emitted spectrum at some surface away from a dipole source, which can be probed, e.g., by near field microscopy techniques and localized detectors. For this latter case, we require the two-space-point Green function, ${\bf G}({\bf r},{\bf r}_{d},\omega)$, projected on a surface point ${\bf r}$ at a different position away from the dipole. This scenario differs in the sense that it also accounts for quenching and propagation between two spatial points, unlike the LDOS. These spectral function also relate to the QNM phase (at different positions), and are complimentary but subtly different. In a traditional normal mode approach, these spectral function would typically be identical, unless additional filtering is included in the input-output theory.

To numerically implement and fully corroborate our QNM theory and model assumptions, we carry out full numerical electromagnetic simulations in COMSOL Multiphysics [58]. The numerical generalized PF (denoted as $F_{\rm P}$), which includes the homogeneous background contribution of $\epsilon_{\rm B}=1$ as well the scattered contributions of all of the modes present in the nanostructure, is calculated by performing full dipole simulations and examining surface-integrated Poynting vector around the dipole at $\mathbf{r}_{\rm d}$ (cf. Fig. 2), from

where the surface $s$ is of a small sphere centred around the finite-size dipole oriented in the $i^{\rm th}$ direction; the small sphere is approximately 1 nm in radius, and $\hat{\bf{n}}$ is the normal vector directed outward relative to the dipole. This allows us to capture the power flow which we have also checked is quantitatively accurate against known solutions in free space. In experiments, this would be measured as the ratio of the spontaneous emission rate of a dipole emitter in the full scattering structure to the rate in a homogeneous background medium. COMSOL can also be used to obtain the QNMs of the system which allows us to obtain the PF everywhere in the system, not just at a single point as in Eq. (13).

The gold nanorod (on the hBN substrate) was designed such that the peak PF of the second lowest mode was located near 1.6 eV to spectrally overlap with the TMDC response. The final dimensions of the nanorod were determined to be $L=148$ nm and $R=10$ nm. The small frequency shift due to the presence of the TMDC, which alters the effective background permittivity, was not considered in this design, but will be discussed later when examining the anti-crossing behavior with respect to detuning.

We next describe how the hybrid QNMs can be obtained to assess and understand the strong coupling regime of the coupled TMDC MNP system and the spectral interaction of QNMs, i. e., through “QNM splitting”. Numerically, we again use COMSOL, since it can work in complex frequency space. We adopt an inverse Green function approach, where the normalized QNM are obtained from the solution to a scattering problem in complex frequency space [59, 60].

The optical QNMs are solutions to the vector Helmholtz equation in the complex frequency domain, which are the proper mode solutions of cavity structures with open-boundary conditions. Assuming non-magnetic ($\mu_{r}=1$) permeability, the electric field QNM eigenvalue equation is

where $\tilde{\mathbf{f}}_{{\mu}}$ are the QNMs, $\tilde{\omega}_{\mu}=\omega_{\mu}-i\gamma_{\mu}$ are the complex eigenfrequencies, and $\mu$ is the mode number. The permittivity, $\mathbf{\epsilon}(\mathbf{r},\omega)$, is completely general and is composed of the Drude gold (Eq. (6)), the anisotropic TMDC model (Eq. (12)), the background permittivity ($\epsilon_{\rm B}$), and hBN substrate permittivity ($\epsilon_{\rm hBN}$). The QNM cavity factors are given by $Q_{\mu}=\omega_{\mu}/2\gamma_{\mu}$.

The electric-field Green function can then be defined in terms of an expansion over the QNMs and a complex frequency prefactor [29],

where the QNMs only depend on spatial position. Note that the form of the Green function uses an unconjugated product, and this is a general feature of QNMs, including the MNP on its own as well as the hybrid modes. The unconjugated product is essential to capture complex interference effects arising from the QNM phase of the hybridized modes, which can also gives rise to Fano-like resonances [61, 62]. The above expansion is known to be highly accurate for spatial locations within and near the resonator [25, 28, 43]. For positions far from the resonator, then the above solution can be used with the Dyson equation to obtain regularized (non-divergent) QNMs [63], even for spatial positions far from the resonator; or one can also use near-field to far-field transformations [64] (see also Ref. 65).

As mentioned earlier, the concept of effective mode volume here represents a measure of $\tilde{\bf f}_{\mu}^{2}$, and is thus a spatially dependent quantity, and it is also complex (for QNMs). The complex mode volume, for mode $\mu$, for an emitter at position ${\bf r}_{\rm d}$ in $\epsilon_{\rm B}$, is given by

The effective mode volume for use in Purcell’s formula can be calculated as the real part of Eq. (16),

Thus, for one QNM, $\mu=c$ (positive frequency pole), we recover the usual one mode solution from the Purcell factor,

where we assume the dipole is aligned with the dominant polarization of the QNM, and $n_{\rm B}=\sqrt{\epsilon_{\rm B}}$. With several modes, there are complex interference effects stemming from the concept of a complex effective mode volume, which is a consequence of the QNM phase. These effects can also be obtained directly from the single QNM Green function, which is in fact much more convenient with coupled QNMs.

For several QNMs, the generalized PF can be obtained from the two QNM Green function ($\mu=1,2$) relative to the homogeneous solution,

where ${\rm Im}[G_{\rm hom}]=\frac{n_{\rm B}}{6\pi\epsilon_{0}}\left(\frac{\omega}{c}\right)^{3}$. The plus one appears due to the total Green function being the sum of the scattered and the homogeneous parts [66].

Figure 4 shows the numerical full-dipole PF (markers) as well as the individual hybrid QNMs (dashed/dotted lines) and the sum of the two QNMs (solid line) for the hybrid TMDC-MNP system at various temperatures of interest (4 K, 77 K, 300 K) as well as the single MNP on hBN (independent of temperature). Considering no other modes are used in the sum (i.e., higher and lower order plasmonic modes that are not spectrally isolated from the frequency range of interest), the agreement between the full numerical simulations and the two QNM simulations (with no fitting parameters) is overall in very good agreement.

Having now confirmed that the dipole PF is well described by the two QNM Green function expansion, we stress that Eq. (15) can be used to obtain the PF at all locations near the resonator, highlighting the remarkable power of the QNM technique. This can be recognized from Figure 4 showing the absolute value of the $x$-component of the individual QNMs, $|\hat{\mathbf{x}}\cdot\tilde{\mathbf{f}}_{\mu}(\mathbf{r})|$, in the $xz$-plane ($y=0$) and $xy$-plane ($z=-0.5$ nm, in the middle of the gap), for a gap of 1 nm. The MNP, substrate, and TMDC are outlined as white lines. The two hybridized modes look qualitatively similar to the bare MNP QNM, but are indeed new modes, not just perturbations of the original. Figure 5 shows the anti-crossing behavior in the total PF, calculated using Eq. (19), for the MNP-TMDC system with a gap of 2 nm which is expected in the presence of strong coupling, showing a minimum Rabi splitting of 84.6 meV and 86.1 meV for temperatures of 4 K and 300 K, respectively. Panel (c) shows the value of the Rabi splitting as a function of detuning, which is achieved by modifying the length of the MNP. In addition to the anti-crossing behavior, $\rm\hbar\Omega_{R}>\hbar(\gamma_{p}+\gamma_{ex})$ is sometimes used as a benchmark for the onset of strong coupling; in this case, these values are 45.1 meV and 71.7 meV, respectively, which are much less than the observed splitting.

The minimum Rabi splitting is observed at a detuning of approximately 10-15 meV, rather than 0 meV. This is due to finite QNM losses as well as the procedure used to design the MNP; the radius and length of the MNP (on a substrate of hBN) were varied until the peak of the PF matched the resonance frequency of the TMDC. However, due to the presence of the TMDC in the system, it slightly shifts the background effective permittivity which blue shifts the MNP resonance. The peak permittivity is larger for lower temperature, explaining why the shift is slightly larger at 4 K relative to 300 K. A 10 meV shift can be achieved by changing the effective substrate permittivity from 4.5 to 4.6, which is reasonable magnitude to expect from the presence of a TMDC layer.

The effect of the gap size between the TMDC and MNP, $z_{\rm gap}$, is shown in Fig. 6 at 4 K and 300 K. The individual QNM PFs are shown as dashed/dotted lines and the total PFs are shown as solid lines. As expected, the observed Rabi spitting increases as the TMDC gets closer to the MNP since the strength of the electric field due to the plasmonic resonance increases dramatically near the metal surface. Interestingly, the Fano-like behavior of the individual QNMs is most pronounced at either large gap sizes or low temperatures. This can be explained by the QNM phase interference between the bare modes, which is more drastic with the inequality of the linewidths of the two modes, and is more pronounced at elevated temperatures.

Next, Fig. 7(a) displays a summary of the complex poles of the hybrid QNMs as a function of temperature and gap size. Also shown is the resonance frequency of the TMDC alone (dashed) and the QNM pole of the MNP alone (marker, ‘$+$’). Since the TMDC alone does not support a QNM, the TMDC resonance is represented by a vertical line. As the gap between the TMDC and MNP increases, the modes become less coupled, resulting in the modes moving towards either the MNP (upper mode) or TMDC (lower mode), until they are completely uncoupled. The upper mode is more MNP-like, while the lower mode is a more pronounced hybridization of the two original modes. The difference in the imaginary part of the modes contributes to greater changes in phase which in-turn result in greater Fano-like lineshapes seen in Fig. 6. In panel Fig. 7(b), the Rabi splitting as a function of gap size for 4 K and 300 K is shown using two different calculations: (i) the splitting due to the poles (dashed line) is calculated as the difference between the real part of the poles; and (ii) the splitting due to the PF (solid line) is calculated as the difference between the frequency locations of the two maximums of the total PF. Since the data is obtained for discrete values of $\rm z_{\rm gap}$, the data is fit to an exponential curve to make the readability of the information more clear. The standard deviations of the fits are shown by the shaded region, which is negligible for the pole-calculated fit. We can see that, overall, the splitting is larger when observed from the PF than from the poles.

As discussed earlier, another experimental observable could be the emitted spectrum to some detector position that is away from a local excitation (such as from a dipole source). The theory to describe the detected field at these points requires the full two-space-point Green function (or propagator), rather than just the projected LDOS. In normal mode theory, such a spectrum would be identical to the one at the dipole location (apart from a scalign constant), but in QNM theory the phase of the QNM changes as a function of position. In addition, there are propagation and quenching effects between two spatial position, than can render these two function different in plasmonic and TMDC systems. Thus, for completeness, here we also describe how the Poynting vector (and thus, the dipole transmission through some arbitrary surface) can be reconstructed using the QNMs in the relevant spectral region of interest (in a spatial region that is not too far from the resonators).

The Poynting vector is defined as,

where $\mathbf{E}$ is the electric field, $\mathbf{H}$ is the auxillary magnetic field ($\mathbf{H}=\mu_{0}^{-1}\mathbf{B}$, where $\mu_{0}$ is the permeability of free space and $\mathbf{B}=(i\omega)^{-1}\bm{\nabla}\times\mathbf{E}$ is the magnetic field). The electric field due to a point dipole source with dipole moment ${\mathbf{d}}$ is constructed from the Green function via,

where we note that the Green function is a tensor of rank 2 (3 by 3 matrix). The total power flow through a surface $T$ is then given by,

The upper surface ($T^{+}$) and lower surface ($T^{-}$) values are shown in Fig. 2, which are monitors in the $xy$ plane with normal vectors pointing outwards. The observed spectral splitting due to the power transmission of a dipole for a gap of 2 nm and temperature of 4 K is 74.1 meV. In contrast, the splitting observed in Fig. 7 due to the complex poles and the PF were 59.5 meV and 73.8 meV, respectively. There is a slightly greater splitting when observed in the transmission of a dipole emitter compared to the Purcell enhancement of the emission of the same dipole emitter. This is a consequence of the changing QNM from the dipole location to the surfaces for the power flow.

Power transmission can typically be measured experimentally by exciting the system with a plane wave or tightly focussed beam [67] and comparing the optical transmission with and without the optical scatterer. This is different from what we have presented here, since our source is a dipole emitter, but the transmission should be similar to a background free measurement. In contrast, experimentally determining the PF requires measuring the change in the radiative lifetime for a quantum emitter of either the photoluminescence spectra [68] or Raman spectra [69]. These differences are subtle but fully captured in our model. The two-space-point Green functions can also be used to compliment other formalism for systems containing two and more quantum emitters, since they account for important photon exchange effects [70, 71]

In summary, we have presented a QNM approach for the description of the modal physics of a hybrid TMDC-MNP system in the strong coupling regime. In this fully 3D approach, there is no fixed $g$, but rather the splitting behaviour is captured from the underlying hybrid modes, using the full 3D geometry of hybrid cavity structure.

We also examined the effect of temperature, gap size, and detuning to illustrate the characteristics of the coupling between the TMDC and MNP. Spectral splittings as large as 90 meV are observed for a gap size of 0.5 nm. Fano-like contributions from the individual QNMs also manifest in the splitting associated with strong coupling. These hybrid modes are indeed the desired dressed-states states of the hybrid system. Knowing the QNMs of the system allows us to examine the coupling as a function of frequency and space, which can be used to examine more exotic and complex systems such as coupled arrays of dipoles or quantum emitters. Moreover, the use of a QNM model is critically important for the quantization of these lossy modes in quantum optics since they contain the full dispersive nature of the modes that arises from complex eigenfrequencies, unlike traditional normal modes that require phenomenologically added dephasing contributions in Jaynes-Cummings-type models [46, 47]. Our approach thus paves the way to describe such complex cavity structures, including the coupling of more quantum emitters, fully accounting for collective effects, as well as radiative and nonradiative decay and transport as mediated from the hybrid modes.

Note added: during completion of this paper, we became aware of related work on the coupling between TMDCs and MNPs [72]. Similar to our work, they consider a MNP cavity that is described by a single QNM which is then coupled to a TMDC with in-plane Coulomb interaction coupled exciton states. The coupling to the TMDC sheets is then carried out using a quantum reaction-coordinate approach to exciton-resonator interactions, and master equations, as well as semiclassical calculations for the scattered field using a Lippmann-Schwinger equation. Unlike our approach which is based on a microscopic material theory to calculate the underlying hybrid QNMs for linear optics (and also for use for quantized QNM master equations [46, 47]), Ref. [72] develops a coupled oscillator model with microscopically determined parameters that can be used to also describe a non-linear optical response.