Dipolar and Quadrupolar Excitons Coupled to a Nanoparticle-on-a-Mirror Cavity

TL;DR

This paper analyzes how dipolar and quadrupolar excitons interact with a nanoparticle-on-a-mirror cavity, revealing differences in strong coupling behavior and emission dynamics using analytical and numerical methods.

Contribution

It introduces a fully analytical two-dimensional transformation optics model to study plasmon-exciton interactions beyond the quasi-static approximation.

Findings

Finite-size effects influence exciton charge distribution.

Strong coupling phenomena differ between bright and dark excitons.

Far-field scattering spectra match numerical simulations.

Abstract

We investigate plasmon-emitter interactions in a nanoparticle-on-a-mirror cavity. We consider two different sorts of emitters, those that sustain dipolar transitions, and those hosting only quadrupolar, dipole-inactive, excitons. By means of a fully analytical two-dimensional transformation optics approach, we calculate the light-matter coupling strengths for the full plasmonic spectrum supported by the nanocavity. We reveal the impact of finite-size effects in the exciton charge distribution and describe the population dynamics in a spontaneous emission configuration. Pushing our model beyond the quasi-static approximation, we extract the plasmonic dipole moments, which enables us to calculate the far-field scattering spectrum of the hybrid plasmon-emitter system. Our findings, tested against fully numerical simulations, reveal the similarities and differences between the strong coupling phenomenology for bright and dark excitons in nanocavities.

Equations

$$\textbf{Q}=\left(\begin{array}[]{cc}Q_{xx}&Q_{xz}\\ Q_{zx}&Q_{zz}\end{array}\right)=\frac{Q}{\sqrt{2}}\left(\begin{array}[]{cc}\sin 2\alpha&\cos 2\alpha\\ \cos 2\alpha&-\sin 2\alpha\end{array}\right).$$

$$P_{\mu }(\omega )$$

$$P_Q(\omega )$$

$$W_{\mu }^{(0)}(\omega )$$

$$W_Q^{(0)}(\omega )$$

$$\varrho^{\prime}=\ln\bigg{(}\frac{2iD\sqrt{\rho(1+\rho)}}{\varrho-is}+1\bigg{)},$$

$$J_i(\omega )=\frac{{\gamma }_i}{2\pi }{P}_i(\omega ),$$

$$(\sqrt{\rho}+\sqrt{1+\rho})^{2n}=\bigg{|}\text{Re}\bigg{(}\frac{\epsilon_{m}(\omega)-\epsilon_{d}}{\epsilon_{m}(\omega)+\epsilon_{d}}\bigg{)}\bigg{|},$$

$$J_i(\omega )=\sum _{n=1}^{\infty }\sum _{\sigma =\pm 1}\frac{(g_i^{n,\sigma }{)}^2}{\pi }\frac{{\gamma }_m/2}{(\omega -{\omega }_{n,\sigma }{)}^2+({\gamma }_m/2{)}^2},$$

$${\omega }_{n,\sigma }=\frac{{\omega }_p}{\sqrt{{ϵ}_{\infty }+{ϵ}_d{\xi }_{n,\sigma }}}$$

$$g_i^{\mathrm{PS}}=\sqrt{\sum _{\sigma =\pm 1}\sum _{n=n_{\text{min}}}^{\infty }(g_i^{n,\sigma }{)}^2},$$

$$P_{\mu }^{\mathrm{ext}}(\omega )$$

$$P_Q^{\mathrm{ext}}(\omega )$$

$${\hat{H}}_{\mathrm{sys}}$$

$$\frac{\partial \hat{\rho }}{\partial t}=i[\hat{\rho },{\hat{H}}_{\mathrm{sys}}]+\sum _{n,\sigma }\frac{{\gamma }_m}{2}{\mathcal{L}}_{{\hat{a}}_{n,\sigma }}[\hat{\rho }],$$

$${\hat{H}}_{\exp }={\hat{H}}_{\mathrm{sys}}+E_{\mathrm{L}}{e}^{-i{\omega }_{\mathrm{L}}t}{\hat{M}}^{†}+E_{\mathrm{L}}{e}^{i{\omega }_{\mathrm{L}}t}\hat{M}$$

$$\frac{\partial {\hat{\rho }}^{\prime }}{\partial t}=i[{\hat{\rho }}^{\prime },{\hat{H}}_{\exp }^{\prime }]+\sum _{n,\sigma }\frac{{\gamma }_{n,\sigma }}{2}{\mathcal{L}}_{{\hat{a}}_{n,\sigma }}[{\hat{\rho }}^{\prime }]+\frac{{\gamma }_i^{\mathrm{r}}}{2}{\mathcal{L}}_{{\hat{\sigma }}_i}[{\hat{\rho }}^{\prime }],$$

$${\gamma }_{n,+1}^{\mathrm{r}}$$

$${\mu }_n=\mathrm{Re}\left\{\frac{{ϵ}_m({\omega }_{n,+1})+{ϵ}_d}{2{ϵ}_m({\omega }_{n,1})}\right\}\sqrt{\frac{3\pi {ϵ}_0\hslash {\gamma }_{n,+1}^{\mathrm{r}}{c}^3}{\sqrt{{ϵ}_d}{\omega }_{n,+1}^3}}.$$

$${\sigma }_{\mathrm{sca}}({\omega }_L)=⟨\hat{M}{⟩}_{\mathrm{SS}}^2=\mathrm{Tr}\{{\hat{\rho }}_{\mathrm{SS}}^{\prime }({\omega }_L)\hat{M}{\}}^2.$$

$${\varphi }_{\mu }^{\mathrm{S}}({\varrho }^{\prime })=$$

$${\varphi }_Q^{\mathrm{S}}({\varrho }^{\prime })=$$

$$\varrho^{\prime}_{E}=\ln\left(\Big{|}1+\tfrac{{2D\sqrt{\rho(1+\rho)}}(ix_{E}+(z_{E}-s))}{x_{E}^{2}+(z_{E}-s)^{2}}\Big{|}\right)++i\arctan\left({\tfrac{x_{E}}{\tfrac{x_{E}^{2}+(z_{E}-s)^{2}}{{2D\sqrt{\rho(1+\rho)}}}+(z_{E}-s)}}\right),$$

$${\mu }_{x^{\prime }}^{\prime }+i{\mu }_{z^{\prime }}^{\prime }$$

$$Q_{x^{\prime }{x}^{\prime }}^{\prime }+i{Q}_{x^{\prime }{z}^{\prime }}^{\prime }$$

$${\varphi }_{\mu }^{\mathrm{sc}}(\varrho ,\omega )$$

$${\varphi }_Q^{\mathrm{sc}}(\varrho ,\omega )$$

$$A_n^{\pm }(\varrho ,{\varrho }_E)$$

$$B_n^{\pm }(\varrho ,{\varrho }_E)$$

$$\Delta=\ln(\sqrt{\rho}+\sqrt{1+\rho})-\textrm{Re}\bigg{\{}\ln\bigg{(}1+\frac{2iD\sqrt{\rho(1+\rho)}}{\varrho_{E}-is}\bigg{)}\bigg{\}}$$

$${\varphi }_{\ell }^{(N)}(\varrho ,\omega )$$

$$A_n^0(\varrho )=\sum _{k=1}^N\mathrm{sign}(q_k)A_n^{+}(\varrho ,{\varrho }_k)\mathrm{and}{B}_n^0(\varrho )=\sum _{k=1}^N\mathrm{sign}(q_k)B_n^{+}(\varrho ,{\varrho }_k)$$

$$g_{\mu }^{n,\sigma }$$

Full text

I Introduction

Plasmonic nanostructures allow tailoring the emission characteristics of microscopic light sources Giannini2011 . The ability of surface plasmons (SPs) to modify the local density of photonic states and the near-to-far-field coupling efficiency of nanoemitters was firstly exploited for the conception, design and optimization of optical nanoantennas Bharadwaj2009 ; Biagioni2012 . In this context , the goal was amplifying both physical quantities to improve the inherent radiative properties of dye molecules and quantum dots, mainly. These nanoantenna-enhanced (faster, brighter or directional) emitters have found applications in areas such as photodetection, nonlinear optics or imaging Novotny2011 .

More recently, the scientific and technological focus in nanophotonics has shifted from the classical to the quantum optical regime Chang2006 ; Tame2013 . In the quest for ultra-compact nonclassical light sources, new design strategies for plasmonic devices are required Fernandez-Dominguez2018 . In order to develop functionalities operating with single photons Hartsfield2015 ; Gong2015 ; Hoang2016 , the interaction between quantum emitters (QEs) and the electromagnetic (EM) fields associated with SPs must be enhanced, and even pushed beyond the weak coupling regime Hummer2013 . This demands nanostructures yielding extremely large and highly structured spectral densities, in a similar way as nanoantennas do. On the contrary, in order to minimize radiation losses, which constraints the QE-SP interaction strength, the near-field of QEs must be effectively decoupled from their far-field. For this reason, plasmonic resonators perform as nanocavities for quantum optical applications Marquier2017 . Importantly, despite of their low quality factor (caused by metal absorption), the small effective volume of SPs give access to an unexplored parametric region of light-matter coupling, not accessible by other photonic technologies Schwartz2011 ; Zengin2015 ; DeGiorgi2018 .

In the strong coupling regime, light and matter excitations mix together, giving rise to hybrid states known as plasmon-exciton polaritons (PEPs). PEP characteristics can be controlled through the weight of their two constituents Chikkaraddy2016 . For quantum nanophotonics applications, this phenomenon makes it possible to tune the balance between the high coherence of SPs, and the high nonlinearities of optical transitions in QEs Torma2015 . Additionally, the small mode volume of SPs open the way to the formation of PEPs at the single QE level at room temperature. Thus, much research efforts have concentrated in this objective in the last years. Experimental evidence of strong coupling in ensembles of very few, or even single, QEs has been reported in various gap nanocavities Baumberg2019 : nanoparticle dimers Santhosh2016 ; Roller2016 , tip-substrate nanoresonators Gross2018 ; Park2019 or nanoparticle-on-mirror (NPoM) configurations Chikkaraddy2016 ; Kleeman2017 ; Leng2018 . Similarly, theoretical advances have revealed the geometric and material conditions most convenient for the realization of plasmonic strong coupling at the single emitter level Savasta2010 ; Manjavacas2011 ; Delga2014 ; Gonzalez-Tudela2014 ; Li2016 , as well as strategies to harness photon correlations in hybrid QE-SP systems Ridolfo2010 ; Gonzalez-Tudela2013 ; Alpeggiani2016 ; Saez-Blazquez2017 .

Apart from the strong light-matter interactions that they enable, SPs bring other opportunities to the emerging field of quantum nano-optics. Their deeply sub-wavelength nature unlocks attributes of QEs that are elusive to propagating EM fields, such as chiral phenomena associated to their polarization Lodahl2017 ; Downing2019 , mesoscopic effects Andersen2010 ; Neuman2018 , individual atomic bonds Benz2016 ; Carnegie2018 or light-forbidden excitons Zurita2002 ; Rivera2017 . Among the latter, special attention has been paid to quadrupolar excitons, as there are several theoretical predictions Filter2012 ; Yannopapas2015 ; Sanders2018 ; Cuartero-Gonzalez2018 that suggest that these can be brought up to time scales comparable to dipolar ones, even in the strong coupling regime.

In this Article, we investigate the interaction between the SPs supported by a NPoM cavity and single QEs of two types: those sustaining dipolar, and those supporting quadrupolar excitons. We present a transformation optics Pendry2012 ; Luo2013 (TO) approach which allows us to obtain analytical expressions for all the physical magnitudes characterizing the hybrid QE-SP system. In this context, TO allows the quantization of the plasmonic modes supported by the NPoM cavity and the parametrization of the QE-SP interaction Hamiltonian, in a similar way as other recently proposed methods Dzsotjan2016 ; Hughes2018 ; Franke2019 . Our approach sheds deep insights into two different plasmon-exciton phenomena: the near-field population dynamics in a spontaneous decay configuration and the far-field scattering spectra under dark-field laser illumination. Our model also accounts for finite-size effects associated with both excitonic charge distributions. Throughout the text, the differences and similarities in the light-matter strong coupling phenomenology for dipolar and quadrupolar transitions are discussed and analyzed. Despite the fact that a full three-dimensional (3D) TO framework for dipolar point-like sources is available Li2016 , we employ here its two-dimensional (2D) version Aubry2011 . There are three reasons justifying this choice. First, only the latter is fully analytical, which is instrumental for the description of quadrupolar transitions Cuartero-Gonzalez2018 and finite-sized QEs. Second, the 2D theory can be pushed beyond purely quasi-static approximation Demetriadou2017 , which allows computing the SP dipole moments and therefore, far-field spectra for the QE-SP system. Finally, the comparison between 2D and 3D results reveals that the 2D treatment reproduces all the phenomenology reported in 3D Li2018 .

The paper is organized as follows: Section II outlines the fundamental aspects of our TO approach. The spectral densities and SP coupling strengths for dipolar and quadrupolar QEs are presented in Sec. III. Section IV describes mesoscopic effects associated with the finite size of excitonic charge distributions. In Sec. V, we study the dynamics of the population exchange between QE and SPs in an spontaneous emission configuration. The calculation of the SP dipole moments and the scattering spectrum of the nanocavity under coherent pumping is described in Section VI. Finally, general conclusions are raised in Sec. VII

II Transformation Optics approach

The system under study is depicted in Figure 1(a): the interaction of a single QE with a NPoM plasmonic cavity characterized by the gap $\delta$ and nanoparticle diameter $D$ (related by the ratio $\rho=\delta/D$). The metal permittivity is given by a Drude fitting to silver, $\epsilon_{m}(\omega)=\epsilon_{\infty}-\omega_{p}^{2}/(\omega(\omega+i\gamma_{m}))$ with $\epsilon_{\infty}=9.7$, $\omega_{p}=8.91$ eV and $\gamma_{m}=0.06$ eV. The cavity is embedded in a dielectric medium with permittivity $\epsilon_{d}=4$, which models a DNA origami scaffolding Kaminska2018 . The QE, which can be located anywhere in the surroundings of the nanostructure, is modelled as a point-like EM source of dipolar or quadrupolar character. It is parameterized by its dipole, $\bm{\mu}$ or quadrupole, Q, moment. In order to describe light-matter interactions in this system, we benefit from the TO formalism previously developed to compute the absorption and scattering cross section of similar systems Aubry2011 . We consider the 2D version of the geometry in Figure 1(a), which has translational invariance along the out-of-plane direction. This simplifies greatly the calculation of the EM fields scattered by the NPoM geometry under point dipole and quadrupole excitations, which can be tackled fully analytically.

Let us introduce first the description of the QE sources in the 2D model. The dipolar exciton is characterized by its vectorial dipole moment $\bm{\mu}=(\mu_{x},\mu_{z})=\mu(\sin\alpha,\cos\alpha)$, where angle $\alpha$ (see Fig. 1(a)). The quadrupolar exciton is defined by its tensorial quadrupole moment

[TABLE]

Gauge invariance allows writing the quadrupolar tensor above in a traceless form, without modifying the light-matter coupling Hamiltonian Novotnybook . Therefore, by symmetry, there are only two independent entries, $Q_{xx}$ and $Q_{xz}$ in the tensor, which can be expressed in terms of the modulus $Q=\sqrt{\sum_{i,j}Q^{2}_{ij}}$ and the angle $\alpha$. The quadrupolar source can be interpreted as composed by 4 point-like charges (each of them of opposite sign to the adjacent ones) placed in the vertices of a square. For $\alpha=0$ ($\alpha=\pi/4$) the edges (diagonals) of the quadrupolar charge distribution are parallel to the $x$ and $z$-axes (see Sec. IV).

We employ the 2D model to obtain the Purcell factor, $P(\omega)$, experienced by each QE, which we will subsequently use to extract the corresponding spectral density and QE-SP coupling strengths. Taking advantage of the nanometric dimensions of the cavity, we can operate in the quasi-static regime to calculate the Purcell spectra for both QEs. These are given by the ratio of the power dissipated by the point-like EM source in the presence of the NPoM cavity and the embedding medium over the power radiated in free space

[TABLE]

where $\textbf{r}=(x,z)$ and $\textbf{ r}_{E}$ is the position of the emitter in the $xz$-plane.

The radiated powers in the denominator of Eqs. (2) and (3) can be computed through the free-space EM Dyadic Green’s function in 2D,

[TABLE]

where \textbf{G}^{(0)}(\textbf{r},\textbf{r}^{\prime},\omega)=\tfrac{1}{4i}\big{[}\textbf{I}+\left(\tfrac{\omega}{c}\right)^{2}\nabla\nabla\big{]}H_{0}^{(1)}(\omega|\textbf{r}-\textbf{r}^{\prime}|/c) and $H_{0}^{(1)}(\cdot)$ is the zero-order Hankel function of the First Kind. The expressions above yield $W_{\mu}^{(0)}(\omega)=\mu^{2}\omega^{3}/(16\epsilon_{0}c^{2})$ and $W_{Q}^{(0)}(\omega)=Q^{2}\omega^{5}/(64\epsilon_{0}c^{4})$. Note that in Eqs. (2)-(5) we have used that the power dissipated by dipole and quadrupole emitters is given by $W_{\mu}=\tfrac{\omega}{2}{\rm Im}\{\bm{\mu}\textbf{E}\}$ and $W_{\mu}=\tfrac{\omega}{2}{\rm Im}\{(\textbf{Q}\nabla)\textbf{E}\}$, respectively, where E stands for the electric field due the source in each case Novotnybook .

Equations (2) and (3) show that, neglecting radiative losses, we only need the quasi-static potentials $\phi_{\mu}(\textbf{r},\omega)$ and $\phi_{Q}(\textbf{r},\omega)$ in order to determine the Purcell spectra. This we do using TO, by applying the conformal map

[TABLE]

where $s=\delta+D\sqrt{\rho}/(\sqrt{1+\rho}+\sqrt{\rho})$, and $\varrho^{\prime}=x^{\prime}+iz^{\prime}$ and $\varrho=x+iz$ denote, respectively, the spatial coordinates (in complex notation) for the original and transformed frames. Under Eq. (6), the cavity maps into a metal-dielectric-metal waveguide of width $d=2\ln(\sqrt{\rho}+\sqrt{1+\rho})$, see Figure 1(b). Importantly the permittivities are not affected by the mapping. The QE sources transform into a periodic array of sources of the same character as the initial ones. They are distributed along $z^{\prime}$-direction with period $2\pi$. Thus, the Purcell spectra calculation reduces to solving Laplace’s equation in the geometry of Fig. 1(b). The solution in the original frame is easily obtained using $\phi_{\mu,Q}(\varrho,\omega)=\phi^{\prime}_{\mu,Q}(\varrho^{\prime}(\varrho),\omega)$. Details of the calculation as well as full analytical expressions for the quasi-static potentials are given in Appendix A.

Figure 2 shows the Purcell spectra for NPoM cavities of 30 nm diameter and different gap sizes: 0.9 nm (blue), 1.8 nm (red) and 2.7 nm (green). Two different QEs, placed at the gap center, $z_{E}=0.5\delta$, are considered: (a) a dipolar exciton oriented along $z$-direction ($\mu_{z}=\mu$) and (b) a quadrupolar exciton with a purely non-diagonal moment ($Q_{xz}=Q_{zx}=Q/\sqrt{2}$). In our parametrization, $\alpha=0$ for both QEs. Analytical predictions obtained from our TO approach (color solid lines) are compared against fully numerical 2D calculations (black dashed lines) using the Laplace’s Equation solver implemented in Comsol MultiphysicsTM. Both set of quasi-static spectra are in very good agreement. There are small discrepancies in the low frequency tail of Fig. 2(b), which we attribute to the failure of the assumption that the photonic environment is fully governed by the SP modes supported by the NPoM cavity. Note that our 2D model slightly underestimates Purcell factors when compared to full 3D calculations obtained for the same geometric parameters.

In Fig. 2(c), we plot the ratio between the quadrupolar and dipolar Purcell factors for the three cavity configurations in the panels above. We observe that $P_{Q}(\omega)$ is several orders of magnitude larger than $P_{\mu}(\omega)$ throughout the spectral window considered. However, the ratio between both magnitudes is largest in two regions, the low frequency tail and around $\omega_{\rm PS}\simeq 2.4$ eV. We anticipate here that the nanocavities do not support SPs in the former, and that the latter corresponds to the plasmonic pseudomode formed due to the spectral overlapping of high-order dark SP modes (see Refs. Delga2014 ; Li2016 ). In this spectral position, which lies in the vicinity of the surface plasmon asymptotic frequency (given by the condition $\epsilon_{m}(\omega)+1=0$ Li2018 ), the reduction of the gap size enhances quadrupolar transitions the most, yielding $P_{Q}/P_{\mu}>10^{5}$ for $\delta=0.9$ nm. In these conditions, the time scales of quadrupole and dipole exciton dynamics become similar, as their decay rates in free-space differ in 5-6 orders of magnitude Rivera2017 .

III Spectral Density and Coupling Strengths

The spectral density $J(\omega)$ contains information about the electromagnetic modes supported by the cavity as well as the coupling strength between each of them and the QE. Thus, it depends on the cavity geometry (diameter and gap size) and permittivity, the exciton characteristics (natural frequency and dipolar/quadrupolar moment) and its position and orientation. It can be expressed in terms of the Purcell factor as

[TABLE]

where $i=\mu,Q$ and $\gamma_{i}$ is the QE decay rate in free-space. Note that the equation above ensures that, in the weak-coupling regime, the general Wigner-Weisskopf theory for spontaneous emission (valid beyond the Markovian approximation) recovers an exciton population decaying monotonically in time with Purcell enhanced decay rate $P_{i}\gamma_{i}$ Petruccionebook . In order to calculate the spectral density experienced by dipolar and quadrupolar excitons in the vicinity of the NPoM cavity, we use the 2D Purcell spectra obtained in the previous section. The decay rates in free-space are taken from 3D calculations Klimov2005 , $\gamma_{\mu}=\omega^{3}\mu^{2}/(3\pi\epsilon_{0}\hbar c^{3})$ and $\gamma_{Q}=\omega^{5}Q^{2}/(360\pi\epsilon_{0}\hbar c^{5})$.

By simple inspection of Equations (25) and (25), we can see that the Purcell spectra and, therefore the spectral density for both dipolar and quadrupolar excitons are composed by a number of SP contributions (labelled with index $n$, that is related to the azimuthal order of the plasmonic mode) whose resonant condition reads

[TABLE]

which reproduces the quasi-static condition for absorption maxima obtained under plane-wave illumination Aubry2011 . Note that Eq. (8) is quadratic, which means that we can identify two different solutions for a given $n$. These correspond to SP modes with different parity (with respect to the NPoM gap cavity), which we label as $\sigma=1$ (even) and $\sigma=-1$ (odd). Note that the term even (odd) refers to the symmetric (antisymmetric) character of the SP electric fields across the gap.

Exploiting the Drude form of the metal permittivity, and using the high quality resonator limit Waks2010 , we can expand the spectral densities as a sum of Lorentzian SP terms of the form

[TABLE]

where $\gamma_{m}$ is the Drude damping rate (note that SP radiative damping is neglected in the quasi-static limit) and

[TABLE]

are the SP frequencies, with $\xi_{n,\sigma}=\frac{(\sqrt{\rho}+\sqrt{1+\rho})^{2n}+\sigma}{(\sqrt{\rho}+\sqrt{1+\rho})^{2n}-\sigma}$. Observe the remarkable similarity of the expression above with its 3D counterpart Li2018 .

Equation (10) shows that for both parities, SPs with increasing $n$ approach the pseudomode frequency, $\omega_{\rm SP}=\frac{\omega_{p}}{\sqrt{\epsilon_{\infty}+\epsilon_{d}}}$. Large $\rho$ provides faster convergence of the SPs frequencies to $\omega_{\rm SP}$. It also reveals that the frequency of even ($\sigma=+1$) modes increase towards this value, whereas it decreases for the odd ($\sigma=-1$) ones. As expected from the transformed geometry in Fig. 1, this phenomenology is equivalent to the dispersion of the SPs supported by metal-dielectric-metal waveguides, where both even and odd bands approach asymptotically the frequency for the single metal-dielectric interface (given by $\omega_{\rm PS}$ above).

The weight of each term in Equation (9) is given by $g_{i}^{n,\sigma}$, the QE-SP coupling strength between the dipolar or quadrupolar exciton ($i=\mu,Q$) and the SP mode of azimuthal order $n$ and parity $\sigma$. These constants contain all the information about the QEs and the SP mode spatial profile. Appendix B presents the analytical expressions for $g_{i}^{n,\sigma}$ that we obtain from our TO approach. Their dependence on the QE position and orientation is analyzed in detail below.

Figure 3 shows the spectral densities for a dipolar (a) and a quadrupolar (b) exciton in the same cavity configuration as the blue line in Figure 2. The QEs are placed at the gap center and their orientation is $\alpha=0$. In agreement with experimental values, we set $\mu=0.55$ e$\cdot$nm. We take $Q=0.75$ e$\cdot$nm2 for the quadrupole moment, which yields a ratio between free-space decay rates $\gamma_{Q}/\gamma_{\mu}\simeq 1\cdot 10^{-6}$ at the center of the frequency window considered, $\omega=2$ eV. Despite this inherent difference between both QEs, the spectral densities are equivalent at the pseudomode position, $J_{\mu}(\omega_{\rm PS})=J_{Q}(\omega_{\rm PS})\simeq 10^{-13}\,{\rm s}^{-1}$. Therefore, the interaction strength between the NPoM cavity and both QEs at $\omega_{\rm PS}\simeq 2.4$ eV is very similar. On the contrary, at lower frequencies, $J_{\mu}(\omega)\gg J_{Q}(\omega)$, which indicates that the quadrupole exciton couples more weakly than the dipolar one to low-order SPs (with small $n$). Note the large contrast at the dipolar SP, $J_{\mu}(\omega_{1,+1})\simeq 10^{2}J_{Q}(\omega\omega_{1,+1})$. In fact, the maximum at $\omega_{1,+1}=1.55$ eV barely stands out of the low-frequency tail of the pseudomode maximum in the quadrupolar spectral density, see Fig. 3(b).

In Figure 4, we study the dependence of the light-matter coupling strengths on the QE orientation for the different SPs supported by the NPoM cavity. Both dipolar ($\mu=0.55$ e$\cdot$nm) and quadrupolar ($Q=0.75$ e$\cdot$nm2) excitons are placed at the gap center $(x_{E},z_{E})=(0,0.5\delta)$. Fig. 4(a) and (b) display $g_{\mu}^{n,+1}$ and $g_{\mu}^{n,-1}$, respectively, versus $n$ and $\alpha$. These conturplots show that the maximum coupling takes place at $n<4$ for both $\sigma$, and that vertical (horizontal) dipolar QEs couple more efficiently to even (odd) SPs. This can be clearly seen in Fig. 4(c) which plot $g_{\mu}$ for the three orientations indicated by arrows in the previous panels. Only for $\alpha=\pi/4$ (green lines), the coupling strength to even (solid) and odd (dashed) SPs are comparable. Note that $g_{\mu}^{n,+1}>g_{\mu}^{n,-1}$ for very low azimuthal order, $n$, in this case. Fig. 4(d) and (e) display coupling strength maps for $g_{Q}^{n,+1}$ and $g_{Q}^{n,-1}$. They exhibit a similar dependence on $n$ and $\alpha$ as their dipolar counterparts. However, two main differences can be observed. First, although $g_{\mu}\simeq g_{Q}$ for large $n$, the maximum coupling is always lower for quadrupolar QEs within the range of geometric and material parameters considered. Second, the peak in $g_{i}^{n,\sigma}$ always takes place at lower $n$ for the dipolar excitons. These two circumstances are apparent in Fig. 4(f), which also shows that for a given $n$, $g_{Q}^{n,+1}\simeq g_{Q}^{n,-1}$ at $\alpha=\pi/8$ (green lines).

Once we have studied the orientation-dependence of QE-SP couplings, we investigate next the impact of the emitter position. We restrict our attention first to the symmetry axis of the cavity ($x_{E}=0$). Figure 5(a)-(b) display $g_{\mu}$ maps as a function of $n$ and $z_{E}/\delta$ between the gap center and the vicinity of the NP surface for both plasmonic parities. We can see that, in accordance with Fig. 4, the light-matter interaction is governed by low-order ($n<4$) even SPs. Note that the QE coupling to these modes barely depends on the emitter position. The associated electric field profile is constant along the NPoM gap. On the contrary, in accordance with the phenomenology reported for full 3D models Manjavacas2011 ; Li2016 , $g_{\mu}$ for both even and odd modes of higher $n$ increases as the QE approach the NP surface. This trend is visible in Fig. 5(c), which evaluates $g_{\mu}^{n,\sigma}$ for three $z_{E}$ values. A similar analysis is presented in Fig. 5(d)-(f) for quadrupolar excitons. The $g_{Q}$ maps reveal that, in contrast to dipolar QEs, the coupling vanishes for SPs with very low $n$, and the maximum takes place for $n>10$. As we have already discussed, the QE only interacts with $\sigma=+1$ modes at the gap center. The emitter displacement towards the NP surface increases both $g_{Q}^{n,+1}$ and $g_{Q}^{n,-1}$. The light-matter interaction is then fully governed by the plasmonic pseudomode. Indeed, the cuts at fixed QE position in Fig. 5(f) show that the coupling to even and odd SPs for large $n$ is maximum, and very similar, at $z_{E}=0.75\delta$. Importantly, the maximum coupling in this panel is higher than in Fig. 5(c). This indicates that, by displacing the emitter from the gap center, the plasmonic interaction for quadrupole QEs becomes larger than for dipolar ones.

We exploit the analytical power of our TO approach further and explore fully the spatial distribution of the QE-SP coupling, $g_{i}^{n,\sigma}(\varrho_{E},\omega)$, in the vicinity of the NPoM geometry. Figure 6 displays strength maps (in log scale) involving the dipolar (a,c) and quadrupolar (b,d) excitons and the lowest-frequency SP $(\omega_{1,+1})$ (a,b) and the plasmonic pseudomode $(\omega_{\rm PS})$ (c,d). We set all the parameters as in Fig. 5. The former mode corresponds to $n=1,\sigma=+1$, the coupling constant for the latter is calculated as Delga2014

[TABLE]

where the minimum order for even/odd parity is set by the condition $|\omega_{\rm PS}-\omega_{n_{\textrm{min}},\sigma}|\leq\gamma_{m}/2$. Notice that $n_{\textrm{min}}=7$ in Fig. 6(c)-(d) for both parities, which is in accordance with Fig. 3, which only shows five distinguishable peaks in $J_{i}(\omega))$ below $\omega_{\rm PS}$.

Figure 6(a) and (b) evidence that the coupling-strength maps associated to the lowest (dipolar) SP are focused within the gap of the NPoM geometry. However, the localization at the gap is significantly larger for the quadrupole QE ($\alpha=0$ for both excitons). Whereas the region of high $g_{\mu}$ spreads over the flat metal surface and the perimeter of the NP, $g_{Q}$ decays abruptly within a few nanometer range from the gap center. Let us remark again that all contourplots are in logarithmic scale. In contrast, Fig. 6(c) and (d) demonstrate that the pseudomode yields coupling maps insensitive to the cavity geometry. These are much more tightly bounded to the metal boundaries, within a sub-nm length scale, both at the NP and substrate surface. The gap does not seem to play any role in the spatial distribution of $g_{\mu}$ and $g_{Q}$, except from the overlapping of their tails across it. In accordance with the top panels, the quadrupole distribution is also more confined than the dipolar one. The remarkable contrast between Fig. 6(a) and (d) reveals that through the exploitation of higher order SP modes and multipolar excitons, spatial resolutions in the light-matter coupling well below the nanometer can be achieved Benz2016 .

IV Finite-size Effects

In this section, we extend our TO approach in order to address the emergence of mesoscopic effects Andersen2010 in the light-matter interactions due to the finite size of the QE Neuman2018 . The extreme confinement of the plasmonic coupling strength maps shown in Fig. 6 suggests that our NPoM cavity is an ideal platform to explore excitonic charge distributions beyond the point-like description of the EM source. As the spatial variation of $g_{\mu}$ or $g_{Q}$ approaches length-scales comparable to the QE dimensions, we can expect that this approximation breaks down. By inspection of Fig. 6, we can anticipate that these finite-size effects are higher for the plasmonic pseudomode than for SPs with low $n$.

A dipole EM source can be depicted as a pair of point-like charges of opposite sign and magnitude $|q|$. The vector between both charge positions is $\bm{\ell}=\mu/|q|(\sin\alpha,\cos\alpha)$ (note that we assume $\mu=|q|\ell$). The Purcell factor will no longer be given by Eq. (2). Instead, it reads now

[TABLE]

where, as the QE dimensions are much smaller than optical wavelengths ($\ell\ll 2\pi c/\omega$), $W_{\mu}^{(0)}(\omega)$ is given by Eq. (4). The quasi-static potential $\phi_{\ell}^{(2)}(\textbf{r})$ describes the EM fields scattered by the NPoM cavity excited by the two opposite charges separated a a distance $\ell$ (note that, for simplicity, we have dropped its frequency dependence above). The analytical expression for the potential generated by any neutral distribution of point-like charges, used to compute $\phi_{\ell}^{(2)}(\textbf{r})$, is provided in Appendix A.

An extended quadrupole source corresponds to a square-shaped distribution of four point-like charges with side vectors $\bm{\ell}=(Q/\sqrt{2}|q|)^{1/2}(\sin\alpha,\cos\alpha)$ and $\bm{\ell_{\bot}}=(Q/\sqrt{2}|q|)^{1/2}(\cos\alpha,-\sin\alpha)$ ($Q=\sqrt{2}|q|\ell^{2}$). The Purcell factor in this case is given by

[TABLE]

where $\bm{\ell}_{\pm}=\bm{\ell}\pm\bm{\ell_{\bot}}$, and $W_{Q}^{(0)}(\omega)$ is given by Eq. (5). The analytical expressions used to evaluate $\phi_{\ell}^{(4)}(\textbf{r})$ can be found in Appendix A.

Figure 7 reveals the complex phenomenology behind mesoscopic effects in QE-SP coupling, which depends very much on the emitter orientation (see sketches in all panels). Dipolar QEs with $\alpha=0$ (a) and $\alpha=\pi/2$ (c) are displayed in the top panels, whereas quadrupolar excitons with $\alpha=0$ (b) and $\alpha=\pi/4$ (d) are shown in the bottom panels. The geometric and material parameters are the same as in Fig. 3 (note that $z_{E}=0.5\delta$). Spectral densities calculated using the point-source approximation (blue) are compared against finite-size charge distributions for different $\ell$: 0.05 nm (brown), 0.4 nm (red) and 0.6 nm (orange). As expected, the former coincides with the point-source spectra in all cases, which proves the validity of Eqs. (LABEL:eq:Pmuext) and (13) in the limit $\ell\rightarrow 0$.

Spectral densities for vertical and horizontal dipoles in Fig. 7(a) and (c) show the opposite dependence on $\ell$, whereas $J_{\mu}(\omega)$ increases for $\alpha=0$, it decreases for $\alpha=\pi/2$. These deviations occur mainly at the pseudomode position, $\omega_{\rm PS}\simeq 2.4$ eV, whereas peaks in $J_{\mu}(\omega)$ at lower (a) and higher (c) frequencies, which are associated to low-order SPs with $\sigma=+1$ and $\sigma=-1$, respectively, are rather insensitive to $\ell$. This is evident in the insets of both panels, which plot the coupling strengths obtained from extended calculations normalized to the point-dipole predictions. Note that they are computed following the same procedure as descried in Sec. III. We can observe that $g^{\rm ext}/g\simeq 1$ for $n<6$, whereas the ratio increases (a) or decreases (c) significantly with $\ell$ for larger $n$. The contrast between both descriptions is maximum at the pseudomode, which allows us to gain insight into our findings through the map in Fig. 6(c) (evaluated for $\alpha=0$). Indeed, we can infer that the coupling enhancement in Fig. 7(a) is due to the fact that the exciton charges approach the metal boundaries as $\ell$ increases for $\alpha=0$, where $g^{PS}_{\mu}$ is maximum. On the contrary, they displace laterally, away from the gap center and towards regions of lower $g^{PS}_{\mu}$ for $\alpha=\pi/2$, yielding the coupling reduction in Fig. 7(b).

The bottom panels of Fig. 7 show that, for QEs located at the gap center ($z_{E}=0.5\delta$), the impact of finite-size effects are smaller for quadrupolar excitons than for dipolar ones. Fig. 6(d) shows that $g^{PS}_{Q}$ is more localized than its dipolar counterpart at the metal boundaries, which explains the insensitivity of both $J_{Q}(\omega)$ and $g^{\rm ext}/g\simeq 1$ to QE dimensions up to $\ell=0.4$ nm for both orientations. Only for $\ell=0.6$ nm (orange lines) deviations from the point-quadrupole approximation become apparent, which again, they take place mainly at the pseudomode frequency, due to the strongly confined character of the EM fields associated to high-order SPs. The spectral density and pseudomode coupling are only slightly lower than the point-quadrupole prediction for $\alpha=0$, while they are significantly higher for $\alpha=\pi/4$. This higher impact of mesoscopic effects in Fig. 7(d) can be attributed to two factors. First, the distance between the nearest point charges in the quadrupole distribution and the metal boundaries are smaller than in Fig. 7(c). Second, by increasing $\ell$, these charges (located along the vertical axis) interact more strongly with the odd ($\sigma=-1$) SPs supported by the cavity, while their counterparts remain along the $z=0.5\delta$ axis, where $g^{PS}_{Q}$ is minimum. Let us also stress that the coupling strength calculations, specially in Fig. 7(d), must be taken carefully. The fact that the pseudomode peak governs completely the spectral density means that the high-Q resonator limit Waks2010 , inherent to the modal decomposition of $J_{i}(\omega)$ in our approach, is not a valid assumption in this case.

V Exciton Population Dynamics

Once we have analyzed the dependence of the spectral density and coupling strengths on the various parameters of the system, we explore next the onset of strong coupling between dipolar and quadrupolar excitons and NPoM cavities. Using our TO approach, we can parameterize the Hamiltonian governing the coherent QE-SP interaction

[TABLE]

where $\hat{\sigma}_{i}$ and $\hat{a}_{n,\sigma}$ are the QE ($i=\mu,Q$) and SP annihilation operators (we take $\hbar=1$). The full density matrix of the system is then given by the master equation,

[TABLE]

where $\gamma_{m}$ is the Drude damping rate, and the Lindblad term $\mathcal{L}_{\hat{a}_{n,\sigma}}[\hat{\rho}]=2\hat{a}_{n,\sigma}\hat{\rho}\hat{a}^{{\dagger}}_{n,\sigma}-\{\hat{a}^{{\dagger}}_{n,\sigma}\hat{a}_{n,\sigma},\hat{\rho}\}$ accounts for the absorption losses experienced by the SP mode with indices $n$ and $\sigma$. Note that, for the moment, and in order to gain insight into the plasmon-exciton coupling phenomenology, we neglect the SP and QE radiative decay. These damping rates are much smaller than $\gamma_{m}$, and they do not affect the results presented below.

We study the temporal evolution of the exciton population, $n_{E}(t)=\langle e,\{0\}_{n,\sigma}|\hat{\rho}(t)|e,\{0\}_{n,\sigma}\rangle$, in a spontaneous emission configuration. Note that $|e,\{0\}_{n,\sigma}\rangle$ stands for the product of the QE excited state and the ground state of all SP modes. Thus, we set the initial density matrix for the system to $\hat{\rho}(t=0)=|e,\{0\}_{n,\sigma}\rangle\langle e,\{0\}_{n,\sigma}|$ and investigate the population dynamics. The exciton population can be calculated through Eq. (15) or, equivalently, through the Wigner-Weisskopf equation Petruccionebook (beyond the Markovian approximation Li2016 ; Li2018 ; Cuartero-Gonzalez2018 ) expressed in terms of the spectral density $J_{i}(\omega)$  Gonzalez-Tudela2014 ; Li2016 . We will pay special attention to the ocurrence of non-monotonic, reversible dynamics in $n_{E}(t)$, which we can relate to the onset of QE-SP strong coupling and the formation of PEPs at the single emitter level.

Figure 8 analyzes $n_{E}(t)$ for a vertically-oriented dipolar QE ($\mu=0.55$ e$\cdot$nm) placed at the gap of a NPoM cavity with $\delta=0.9$ nm and $D=30$ nm. Fig. 8(a) plots the spectral density at two different positions along the symmetry axis of the structure, $z_{E}=0.5\delta$ (solid line) and $z_{E}=0.9\delta$ (dashed line). It shows a significant enhancement in $J_{\mu}(\omega)$ as the emitter approaches the metal boundaries. Fig. 8(b) and (c) display the exciton population as a function of time and $z_{E}$ for two different QE frequencies (indicated by vertical arrows in the top panel): $\omega_{E}=\omega_{1,+1}=1.55$ eV and $\omega_{E}=\omega_{\rm PS}=2.4$ eV, respectively. If the QE is at resonance with the lowest-frequency SP (b), $n_{E}(t)$ undergoes a smooth monotonic decay. Importantly, this trend barely depends on the QE position. Taking into account the uniform $g_{\mu}^{1,+1}$ map in Fig. 6(a), we can conclude that the QE-SP interaction is governed by this mode in this case. On the contrary, when the QE is resonant with the plasmonic pseudomode, $n_{E}$ varies significantly, see  (8)(c). As expected from $g_{\mu}^{PS}$ distribution in Fig. 6(b), displacing the QE away from the gap center translates into a faster decay initially, and in the occurrence of Rabi oscillations in $n_{E}(t)$ for $z_{E}>0.7\delta$. Note that their pitch, the Rabi frequency, diminishes as $z_{E}$ increases further. They reveal the occurrence of QE-SP strong-coupling, and the formation of PEPs, the new eigenstates of the system.

Fig. 8(d) and (e) explore in a comprehensive manner the dependence of $n_{E}(t)$ on the QE natural frequency. The former corresponds to $z_{E}=0.5\delta$, and shows that reversible dynamics does not take place at any $\omega_{E}$ in this configuration. Notice though that the decay rate increases abruptly when the emitter is at resonance with a SP mode. This is particularly evident at frequencies approaching $\omega_{\rm PS}$. The latter is evaluated at $z_{E}=0.9\delta$ and unveils the emergence of reversible dynamics in the QE population. The Rabi oscillations become specially apparent in the vicinity of the plasmonic pseudomode, where the evolution of the QE population within the first 50 fs exhibits 5 well-defined maxima ($n_{E}>0.6$) and minima ($n_{E}\simeq 0$). Note that the emergence of these oscillations is accompanied by an underlying faster decay of the QE population, which can be linked to strong-coupling version of the phenomenon of quenching, widely investigated in the weak-coupling regime Delga2014 .

Figure 9 reproduces the study in Fig. 8 but for quadrupole QEs ($Q=0.75$ e$\cdot$nm2, $\alpha=0$). Fig. 9(a) evidences the higher sensitivity of the quadrupole spectral density on the emitter position. Whereas several SP maxima are apparent at the gap center (solid line), $J_{Q}(\omega)$ is completely governed by the pseudomode at $z_{E}=0.9\delta$ (dashed line). Fig. 9(b) and (c) reveal that, in agreement with the $g_{Q}$ contourplots in Fig. 6(c)-(d), $n_{E}(t)$ depends more strongly on $z_{E}$ than it does for dipolar QEs. The QE-SP interaction remains in the weak-coupling regime for $\omega_{E}=\omega_{1,+1}$ (b), although the decay rate experiences a strong reduction as $z_{E}$ increases (note the small oscillations at large $z_{E}$, which can be linked to the spectral detuning between the QE frequency and $\omega_{\rm PS}$). On the contrary, oscillations in $n_{E}(t)$ take place when the QE is only slightly displaced from $z_{E}=0.5\delta$ for $\omega_{E}=\omega_{\rm PS}$ (c). The system enters the strong-coupling regime in this case, yielding a clear reduction in the Rabi frequency as the emitter position approaches the metal surface. Fig. 9(d) shows that, regardless of $\omega_{E}$, the quadrupolar QE at the gap of the cavity always experience a monotonic decay (highly Purcell enhanced at the pseudomode). In contrast, Fig. 9(e) proves that in the gap boundaries, $n_{E}(t)$ develops Rabi oscillations for all QE frequencies. As expected, their pitch depends only moderately on $\omega_{E}$, as the QE-SP interaction is fully determined by the plasmonic pseudomode.

VI Scattering Spectrum

After exploring QE-SP strong-coupling through the temporal evolution of the exciton population, we turn our attention into the emergence of PEP signatures in far-field magnitudes, which are accessible experimentally. Specifically, we model a dark-field spectroscopy setup Lei2012 ; Ojambati2019 , in which the system is pumped coherently by a laser field with amplitude $E_{L}$ polarized along $z$-direction and with frequency $\omega_{L}$. The Hamiltonian describing such experimental configuration is

[TABLE]

where $\hat{H}_{\rm sys}$ is given by Eq. (14). The dipole moment operator of the QE-SP system is $\hat{M}=\sum_{n}\mu_{n}\hat{a}_{n,+1}+\mu\hat{\sigma}_{\mu}$, for dipolar QEs, and $\hat{M}=\sum_{n}\mu_{n}\hat{a}_{n,1}$, for quadrupolar ones. Note that only even SPs ($\sigma=+1$) contribute to the dipole moment of the NPoM cavity, and that we neglect the laser excitation of the quadrupole QE. Importantly, due to the pumping terms, we have $\hat{H}_{\rm exp}=\hat{H}_{\rm exp}(t)$. This temporal dependence can be removed under an unitary transformation (see for example Ref. Saez-Blazquez2018 for more details), obtaining $\hat{H^{\prime}}_{\rm exp}\neq\hat{H^{\prime}}_{\rm exp}(t)$ in the laser rotating frame.

In order to compute the dark-field scattering signal, we account for the radiative losses associated to both SPs and QEs in the master equation describing the dark-field setup,

[TABLE]

where $\hat{\rho^{\prime}}$ is the density matrix in the rotating frame. Note that Eq. (16) incorporates the radiative decay of both SPs and QE, which were absent in Eq. (14). The QE radiative decay rates above account for the effect of the embedding dielectric, having $\gamma_{i}^{\rm r}=\sqrt{\epsilon_{d}}\gamma_{i}$ ($i=\mu,Q$), where $\gamma_{i}$ are the decay rates in vacuum (see Sec. III). The SPs decay rate have a non-radiative and a radiative component, $\gamma_{n,\sigma}=\gamma_{m}+\gamma^{\rm r}_{n,\sigma}$. The latter is computed by introducing radiative corrections in our TO approach. Using a procedure very similar to the one presented in Ref. Demetriadou2017 for our NPoM geometry, we obtain

[TABLE]

and $\gamma^{\rm r}_{n,-1}=0$. Once $\gamma_{n,+1}^{\rm r}$ are known, the even SP dipole moments can be obtained by means of the method of images Aubry2011 , having Cuartero-Gonzalez2018

[TABLE]

Once Eqs. (16) and (17) are fully parameterized by means of TO analytical calculations, we can obtain the steady-state density matrix of the system by solving $\partial\hat{\rho}^{\prime}_{\rm SS}/\partial t=0$. Subsequently, we can compute the scattering cross section through the square of the expectation value of the total QE-SP dipole moment operator,

[TABLE]

In the Supplemental Material of Ref.Cuartero-Gonzalez2018 , we show that Eq. (20), restricted to the first excitation manifold and in the limit of low pumping ($E_{L}\rightarrow 0$), reproduces $\sigma_{\rm sca}$ for bare NPoM cavities.

After a brief description of our calculation of far-field spectra, we investigate next the scattering properties of the QE-SP hybrid systems considered in Sec. V. Mimicking the experimental configuration, we fix the QE frequency at resonance with the lowest, brightest SP mode, for which $\omega_{1,+1}=1.55$ eV and $\mu_{1}=46$ e$\cdot$nm, and focus in a narrow spectral window around it. Note that the exciton population dynamics reveal that the strongest QE-SP coupling takes place at the plasmonic pseudomode, but the dark and absorptive character of high-order SPs hampers its probing through the far-field spectrum around $\omega_{L}=\omega_{\rm PS}$. We anticipate though that, despite the significant detuning, strong-coupling signatures to high-order SPs can be recognized in $\sigma_{\rm sca}(\omega_{L})$ at lower frequencies.

Figure 10 shows normalized scattering spectra for a vertically-oriented dipolar QE. The normalization of the cross section is set so that $\sigma_{\rm sca}(\omega_{1,+1})=1$ for the bare NPoM structure (in absence of QEs). Grey dashed lines correspond to the bare cavity, and blue solid ones plot $\sigma_{\rm sca}(\omega_{L})$ when the QE is at the gap center (the spectra are the same in both panels). The former present a symmetric maximum centered at the SP frequency. The latter exhibit a well-defined Rabi doublet structure, with a central minimum at $\omega_{L}=\omega_{1,+1}$ and two maxima at the upper (U) and lower (L) PEP frequencies Savasta2010 ; Manjavacas2011 ; Saez-Blazquez2017 . This splitting is considered the fingerprint of QE-SP strong-coupling regime, and has been thoroughly analyzed in recent experimental reports on gap plasmonic cavities Zengin2015 ; Chikkaraddy2016 ; Santhosh2016 .

The gap-center spectrum (blue solid line) in Fig. 10 is clearly asymmetric, as the maximum below the SP frequency (LPEP) is significantly higher than the one above it (UPEP). In Fig. 10(a), we analyze the origin of this asymmetry, recently reported in different plasmonic systems Saez-Blazquez2018 ; NeumanOpt2018 . The cyan dashed line plots $\sigma_{\rm sca}(\omega_{L})$ obtained by considering only the lowest SP and the plasmonic pseudomode in the cavity field. We can observe that the doublet is symmetric in this case, which allows us to conclude that the height difference of the peaks in the full calculation originates from the interaction between the QEs and even SP modes with azimuthal indices between 2 and $n_{\rm min}$ (see Eq. (11)). The red solid line in this panel is evaluated for a QE displaced away from the gap center by $0.4D=12$ nm along $x$-direction. Fig. 6(a) shows that $g^{1,+1}_{\mu}$ is much lower in this position. The spectrum overlaps with the bare cavity, except in the vicinity of $\omega_{1,+1}$, where it develops a Fano-like profile, Ridolfo2010 , characteristic of the weak, or possibly the intermediate Leng2018 , coupling regime. The red vertical arrows indicate the PEP frequencies (eigenfrequencies of the Hamiltonian in Eq. (14)) in this configuration. Note that their separation is of the order of $\gamma^{\rm r}_{\mu}\simeq 10\,\mu$eV. The sharp dip in the spectrum is a consequence of the weak, coherent interaction between QE and cavity Ridolfo2010 . The right inset plots the square of the Hopfield coefficients for the LPEP (squares) and the UPEP (circles) as a function of $x_{E}/D$. These give the PEP content on the dipolar QE, $n_{E}=\langle e,\{0\}_{n,\sigma}|\hat{\rho}_{\rm SS}(\omega_{\rm PEP})|e,\{0\}_{n,\sigma}\rangle$, (blue dots) and lowest SP mode, $n_{1,+1}=\langle g,1_{1,+1}|\hat{\rho}_{\rm SS}(\omega_{\rm PEP})|g,1_{1,+1}\rangle$ (yellow dots). They show that, due to the reduction experienced by the QE-SP coupling, the lower (upper) PEP becomes more QE-like (SP-like) as $x_{E}/D$ increases.

Figure 10(b) explores the effect of moving the dipolar QE vertically. Red solid line plots $\sigma_{\rm sca}(\omega_{L})$ for emitter positions very close to the metal surface ($z_{E}=0.9\delta$). We can observe that both the Rabi splitting and the difference between LPEP and UPEP scattering maxima remain very similar to the ones at the gap center. On the contrary, the whole doublet structure has shifted significantly to lower frequencies (note that the scattering minima is no longer at $\omega_{1,+1}$). Orange dashed line plots the same spectrum but considering only the lowest SP and the pseudomode in the evaluation of Eq. (20). The position of the doublet is the same as in the full calculation but, once again, the asymmetry in the peaks height has vanished. This fact agrees with our interpretation, which links the differences in the scattering maxima with intermediate ($2\geq n\geq n_{\rm min}$) even SP modes. This approximate spectrum exhibits the same redshift as the exact one. Taking Fig. 6(a) and (c) into account, we can attribute this shifting of the Rabi doublet to the stronger coupling between the QE and the plasmonic pseudomode caused by the vertical displacement. The squared Hopfield coefficients in the right inset of this panel show that, similarly to Fig. 10(a), the balance between $n_{E}$ and $n_{1,+1}$ in both PEPs is lost as $z_{E}$ increases. Importantly, in contrast to the lateral displacement, this unbalance is accompanied here by an exponential growth of $n_{\rm PS}=\langle g,1_{\textrm{PS}}|\hat{\rho}_{\rm SS}(\omega_{\rm PEP})|g,1_{\textrm{PS}}\rangle$ (green dots). This verifies that, indeed, the redshift experienced by the scattering features originates from the stronger interaction between the QE and the plasmonic pseudomode. In fact, it can be interpreted as a result of the anticrossing between the UPEP and another, even higher frequency, PEP (not analyzed here) that is located at $\omega_{\rm PS}$ in the limit of low QE-SP coupling Cuartero-Gonzalez2018 .

The left insets in Fig. 10(a)-(b) plot spectra for the same system configurations as those rendered in solid lines in the main panels, but replacing the radiative decay rate $\gamma_{\mu}^{\rm r}$ in Eq.(17), which is of the order of $0.5\,\mu$eV (at $\omega_{E}=2$ eV), by a much larger rate $\gamma_{\mu}=20$ meV. This way, we account for both the non-radiative decay and dephasing experienced by QEs in a very simplified fashion Saez-Blazquez2017 . Note that the linewidth of J-aggregates, which exhibit sharp absorption and emission bands, are of the order of 10-30 meV at room temperature Valleau2012 . Therefore, by introducing a more realistic description of QEs, our approach yields smoother scattering spectra, which resemble those in recent experimental reports Chikkaraddy2016 ; Santhosh2016 ; Roller2016 ; Gross2018 ; Park2019 . Importantly, the physical discussion regarding the fingerprint of PEP formation in $\sigma_{\rm sca}(\omega_{L})$ above remains valid for these insets.

Figure 11 presents an analysis similar to the one in Fig. 10 but for quadrupolar QEs (note the narrower frequency window). Grey dashed and blue solid lines in both panels plot the bare cavity cross section and the spectrum for $z_{E}=0.5\delta$ (all system parameters are the same as in Fig. 9). In contrast to its dipolar counterpart, the gap-center spectrum does not exhibit a Rabi doublet, but a sharp Fano profile at $\omega_{L}=\omega_{1,+1}$. This is due to the fact that $g_{Q}^{1,+1}\ll g_{\mu}^{1,+1}$ at the gap center, see Fig. 6. As a result of the weaker coupling, the system eigenfrequencies (see vertical arrows) are close together in this case. Cyan dashed lines in Fig. 11(a) render $\sigma_{\rm sca}(\omega_{L})$ for the same configuration, but including only the lowest SP and the pseudomode in the calculation. The deviations from the exact result are apparent mainly at the scattering minimum, which reveals that intermediate SPs play a more minor role than in dipolar QEs. Red solid line plots the $\sigma_{\rm sca}(\omega_{L})$ for $x_{E}=0.4D=12$ nm. As expected from the tight localization of $g_{Q}^{1,+1}$ at the NPoM gap in Fig. 6(b), this spectrum coincides with the scattering cross section of the bare cavity, as QE-SP interactions vanish in this position. The square of the PEP Hopfield coefficients in the right inset shows that the system remains in the weak-coupling regime for all $x_{E}/D$ values. They demonstrate that the LPEP (UPEP) collapses rapidly into the quadrupole exciton (lowest SP mode) as the emitter moves away from the center of the gap.

The sensitivity of the scattering cross section to variations in the vertical position of the quadrupolar QE is analyzed in Fig. 11(b). Red solid and orange dashed lines plot $\sigma_{\rm sca}(\omega_{L})$ at $z_{E}=0.85\delta$ obtained from the full NPoM plasmonic spectrum and including only the lowest SP and pseudomode contributions in the calculation, respectively. The differences between them are even smaller than at the gap center, see Fig. 11(a). We can observe that, by approaching the emitter to the metal surface, the Fano-like profile in the blue solid line shifts to lower frequencies, but no Rabi doublet structure emerges in the spectrum. This means that the interaction between the QE and the lowest SP remains in the weak-coupling regime, despite the enhancement experienced by their coupling strength. Fig. 6(b)-(d) reveal that $g_{Q}^{\rm PS}$ grows much faster than $g_{Q}^{1,+1}$ with $z_{E}$, which explains why the main effect observed in $\sigma_{\rm sca}(\omega_{L})$ is the red-shift of the Fano feature. Again, this occurs due to the anticrossing with another PEP, whose initial content is mainly pseudomode Cuartero-Gonzalez2018 . The square of the LPEP Hopfield coefficients in the right inset shows that for larger $z_{E}/\delta$, $n_{1,+1}$ decreases, while $n_{\rm PS}$ increases, modifying the inherent character of this polariton, which now emerges from the hybridization of the QE exciton and the plasmonic pseudomode. On the contrary, the UPEP collapses into the lowest, bright SP in this process, decoupling completely from the quadrupole QE.

As in the exploration of dipolar QEs, we have introduced two left insets in Fig. 11(a)-(b), which render the scattering cross sections for the same system configurations as those plotted in solid lines in the main panels, but for a non-vanishing quadrupolar QE decay rate, $\gamma_{Q}=2$ meV. This allows us to show how the sharpness of the Fano-like spectral features in the main panels is reduced once a finite linewidth is introduced in the QE model.

VII Conclusions

We have presented a transformation optics approach that exploits two-dimensional conformal mapping to obtain a full analytical, insightful description of plasmon-exciton interactions in a nanoparticle-on-a-mirror cavity. Two different quantum emitters, supporting only dipolar or only quadrupolar transitions, have been thoroughly analyzed and compared. We have firstly computed the nanocavity spectral densities for both emitter families, which can be decomposed in terms of lorentzian contributions. This enables us to identify the plasmon-exciton coupling strengths for the full nanocavity electromagnetic spectrum, which becomes naturally quantized. Next, we have characterized in detail the dependence of plasmon-exciton coupling strengths on the emitter position and orientation. Special attention has been paid to mesoscopic effects taking place when the dimensions of the exciton charge distribution are comparable to the gap of the structure. Finally, the onset of the strong-coupling regime and the formation of plasmon-exciton polaritons has been investigated in two different, complementary, studies. First, we have revealed the occurrence of Rabi oscillations in the temporal evolution of the exciton population in a spontaneous emission configuration. Second, we have shown the emergence of a Rabi doublet structure in the dark-field scattering spectrum of the nanocavity-emitter system under laser illumination. We believe that our findings can serve as a guidance for the design and interpretation of experiments aiming to harness plasmon-exciton strong-coupling phenomena at the single emitter level.

Acknowledgement

This work has been funded by the Spanish MINECO under contracts FIS2015-64951-R, MDM-2014-0377-16-4 and through the “María de Maeztu” programme for Units of Excellence in R&D (MDM-2014-0377), as well as the EU Seventh Framework Programme under Grant Agreement FP7- PEOPLE-2013-CIG-630996.

Appendix A Quasi-static Potential Calculation

Here, we provide the analytical expressions behind the Purcell factor and spectral density calculations in Section II and  III, as well as the scattering potential for a neutral multiple point-charge distribution introduced in Section IV.

The source potentials describing the array of transformed point-sources in Fig. 1(b) can be written as

[TABLE]

where

[TABLE]

and the transformed dipolar and quadrupolar moments have the form (for small distances)

[TABLE]

Note that the conformal nature of the mapping in Eq. (6) preserves the character of the original excitation potential: a single dipole (quadrupole) source transforms into a periodic array of identical dipole (quadrupole) sources.

In order to solve Laplace’s Equation in the transformed frame and obtain the total quasi-static potentials, we Fourier transform Eqs. (21) and (22). Then, we impose that the scattered potentials have the same spatial dependence as the propagating SPs sustained by the metal-dielectric-metal geometry. Applying continuity conditions at the metal-dielectric interfaces, and performing an inverse Fourier transform, we obtain the scattered potentials in transformed space. Using $\phi^{\rm sc}_{\mu,Q}(\varrho,\omega)=\phi^{\prime\rm sc}_{\mu,Q}(\varrho^{\prime}(\varrho),\omega)$, we obtain their analytical expression in the NPoM frame

[TABLE]

and

[TABLE]

where

[TABLE]

and

[TABLE]

The total potentials can then be written as $\phi^{\rm tot}_{i}(\varrho)=\phi_{i}^{\rm S}(\varrho)+\phi^{\rm sc}_{i}(\varrho)\simeq\phi^{\rm sc}_{i}(\varrho)$ with $i=\mu,Q$.

Finally, the scattered potential for a neutral distribution of $N=2,4$ point charges $q_{k}$ (with $|q_{k}|=|q|$) located at positions $\varrho_{k}$, separated by distances given by the displacement $\ell$ reads

[TABLE]

with

[TABLE]

Appendix B Expressions for light-matter coupling strengths

By reshaping the quasi-static potentials in Eqs. (25), (LABEL:eq:pot_Q_sc) and  (31), we obtain the following analytical expressions for the light-matter coupling strengths that weight the various SP contributions to the spectral density in Equation (9)

[TABLE]

where $\alpha$ is the angle defining the QE orientation and

[TABLE]