The interaction between excited states of a molecule and excited states of a metal nanostructure (e.g., plasmons) leads to hybrid states with modified optical properties. When plasmon resonance is swept through molecular transition frequency, an avoided crossing may be observed, which is often regarded as a signature of strong coupling between plasmons and molecules. Such strong coupling is expected to be realized when 2|⟨*U*⟩|/*ℏ*Γ > 1, where ⟨*U*⟩ and Γ are the molecule–plasmon coupling and the spectral width of the optical transition, respectively. Because both ⟨*U*⟩ and Γ strongly increase with decreasing distance between a molecule and a plasmonic structure, it is not obvious that this condition can be satisfied for any molecule–metal surface distance. In this work, we investigate the behavior of ⟨*U*⟩ and Γ for several geometries. Surprisingly, we find that if the only contributions to Γ are lifetime broadenings associated with the radiative and nonradiative relaxation of a single molecular vibronic transition, including effects on molecular radiative and nonradiative lifetimes induced by the metal, the criterion 2|⟨*U*⟩|/*ℏ*Γ > 1 is easily satisfied by many configurations irrespective of the metal–molecule distance. This implies that the Rabi splitting can be observed in such structures if other sources of broadening are suppressed. Additionally, when the molecule–metal surface distance is varied keeping all other molecular and metal parameters constant, this behavior is mitigated due to the spectral shift associated with the same molecule–plasmon interaction, making the observation of Rabi splitting more challenging.

## I. INTRODUCTION

Plasmonic excitations occurring in small nanoparticles provide an interesting avenue for light and energy manipulation below the diffraction limit. Collective excitations of the electrons in the nanoparticle can couple to nearby molecules, and applications in nano-optics,^{1–4} nanosensing,^{5–9} nanocatalysis,^{10–20} plasmonic photosynthesis,^{21–23} quantum information,^{24,25} and energy-harvesting technologies^{26–31} have emerged.

In standard molecular spectroscopy and photochemistry, the molecule–photon coupling weighted by the density of photon states is weak, as exemplified by the relatively slow spontaneous emission rates ($\u223c109s\u22121$ for allowed electronic transitions), although using a high-intensity radiation field can lead to a strong and often non-linear molecular response. Different situations are encountered when the molecule interacts with one or a few localized photon mode(s), as encountered in optical cavities.^{32–37} In such situations, the system response is best described in terms of hybrid radiation-matter states (polaritons) that diagonalize the strongly coupled part of the Hamiltonian, which include the molecule, these localized modes, and their mutual interaction. The most prominent manifestation of this hybridization is the avoided crossing (Rabi splitting, Ω_{R}), observed in the optical response of such systems when the molecular transition frequency approaches and crosses a cavity mode frequency. The coupling strength reflected in the observed splitting also depends on the number of molecules *N* and is given by $\Omega R=2gN$ if the molecular aggregate is much smaller than the cavity mode wavelength. Here, *g* is the coupling parameter of a dipole **D** (representing the transition dipole of a single molecule) to the radiative mode, given for a rectangular cavity of volume *V* and dielectric constant *ɛ* = *ɛ*_{r}*ɛ*_{0} by *g* = *D*(*ℏω*/2*ɛV*)^{1/2}. By this experimental measure, “strong coupling” is often defined by the observability of this spectral structure, namely, by the condition Ω_{R}/Γ > 1, where Γ is the width of the polaritonic peaks. The possible consequences of strong coupling for transport phenomena^{38–46} and chemical rates^{29,35,47–70} have been subject to intense discussion over the last decade.

Along with Fabry–Pérot (FP)-type cavities, where the localization length is limited by the optical mode wavelength and is bounded by the cavity size, strong coupling phenomena are also observed in interacting plasmon–exciton systems.^{71–74} In particular, metallic nanostructures offer alternative setups for localization of electromagnetic modes, where considerably higher dissipation rates due to losses in the metal constituents are compensated by the considerably smaller, sub-wavelength localization volumes. Intermetallic gaps between plasmonic structures are known to support particularly strong localization. Such localities have been identified as “hotspots” in studies of surface-enhanced spectroscopies^{75–78} and are now referred to as plasmonic cavities.^{79–84} A cavity-like structure is not, however, an absolute requirement. It was recently pointed out that strong coupling, as defined above, may be realized even in the vicinity of single plasmonic particles.^{85–88}

In addition to the usual spectroscopic implications (i.e., the observed Rabi splitting), other dynamical processes may be modified in such strong coupling situations. Photochemical processes can be enhanced or induced at plasmonic interfaces due to the focused nature of the local EM field^{89,90} or as a consequence of plasmon-induced generation of hot electrons.^{91,92} Moreover, strong coupling induced modification of inter- and intra-molecular interactions may lead to new or modified chemistry even without incident light.^{65,93–95}

Compared to FP cavities where many emitters are needed to achieve strong coupling, the much smaller volume of plasmonic cavities makes it possible to observe this phenomenon down to a single nanodot^{80,88,96–98} and even a single molecule.^{85,99–101} In this case, averaging over molecular position and orientation is not an option and the dependence on these geometrical parameters should be addressed explicitly. Furthermore, the geometry of the plasmonic structure can be used to tune the plasmonic resonance. It is important to note, however, that a practical manifestation of “strong coupling” is a relative concept: the possibility to observe Rabi splitting in the optical response of a metal-plasmonic structure composite reflects not only the magnitude of the molecule–radiation field coupling, an intrinsic property of this composite, but also the linewidth, which can have many origins, including spectral connection within the molecule and its interaction with the thermal environment. Limiting ourselves to lifetime broadening, it is often the case that, at close proximity to a metal structure, this attribute of the optical response is also dominated by the molecule–metal interaction, namely, by the effect of the metal structure on the radiative and non-radiative relaxation of the excited molecule. It is therefore useful to examine the relative magnitudes of the molecule–plasmon coupling strength and the metal induced lifetime broadening as providing an intrinsic bound to the presence of strong coupling as manifested by the observability of strong coupling. A recent observation of Rabi splitting in the zero-phonon line of a single molecule located in a plasmonic cavity at low temperature^{102} suggests that “strong coupling” may be more pervasive in such systems and can be observed if other sources of spectral line broadening are eliminated.

In this paper, we study the emergence of light–matter strong coupling near several prototype plasmonic structures, aiming to map the phenomenon as a function of nanostructure geometry and molecular position and orientation. We examine the dependence of strong coupling as experienced by a single molecule, represented by a point-dipole and positioned near metallic spheres, ellipsoids, and bispherical dimers, on these structural parameters. While light–matter interaction has important quantum ramifications, the information needed to establish the occurrence of strong coupling can be obtained from classical considerations and we invoke this approach in the present study. Furthermore, our focus on nanoparticles and their close vicinity makes it possible to ignore retardation effects and work in the electrostatic limit.

The models used in this study and our method of calculation are described in Sec. II. Details of the calculations performed for the different models are provided in Sec. III and in the supplementary material. Some representative results are presented and discussed in Sec. IV. Section V concludes.

## II. MODELS AND METHOD

Figure 1 displays the models used in this study. The plasmonic system comprises of either a metallic ellipsoid (that includes a sphere as a limiting case), shown in Fig. 1(a), or a metallic bispherical dimer, shown in Fig. 1(b)—systems for which analytical solutions to the relevant electrostatic boundary value problem can be obtained. The metal is represented as a continuum dielectric whose dielectric response function is given in the Drude form

where *ɛ*_{b} is the background permittivity from the bound electrons in the metal, *ω*_{p} is the bulk plasma frequency, and *γ*_{p} is a phenomenological damping constant. In the calculations presented below, the following parameters were used: for gold,^{103}^{,}*ɛ*_{b} = 9.84, *ω*_{p} = 2.872 eV = 4.363 × 10^{15} rad/s, and *γ*_{p} = 0.072 eV = 1.094 × 10^{14} rad/s and for silver,^{104}^{,}*ɛ*_{b} = 5, *ω*_{p} = 3.98 eV = 6.05 × 10^{15} rad/s, and *γ*_{p} = 0.039 eV = 5.93 × 10^{13} rad/s. These metal particles are embedded in a uniform dielectric environment with a dielectric constant *ɛ*_{s}, taken to be 1. Other model dielectric functions can, of course, be used. In particular, to include effects of interband transition in the metal or the motion of the ionic core, one often considers Drude–Lorentz-type dielectric response functions.^{105} The molecule is taken as a point-dipole characterized by its magnitude and its position and orientation relative to the plasmonic structure. Finite molecular size effects can be considered as described in Ref. 105 but are disregarded in the present calculation.

The consequences of coupling between a dipole emitter and a plasmonic nanostructure can be calculated in terms of the dyadic Green’s function $G\u0332(r,r\u2032,\omega )$ that describes the response at **r** to a point-current positioned at a given location **r**′ and oriented relative to this plasmonic structure.^{106} In particular, for a dipole oscillating at frequency *ω*, this Green’s function yields the field associated with the oscillating dipole, leading to the interaction between the dipole and the plasmonic system in the form $U\u0304=\u2212D\u22c5E$, where **D** = **D**_{ge} is the transition dipole between the relevant states *g* and *e* of the molecular emitter and its direction reflects the molecular orientation in space and **E** is the electric field associated with the response of the plasmonic structure to this dipole, calculated at the dipole position.^{107} In the quantum case, **D** and **E** are replaced by their quantum operator counterparts, but classical considerations are known to be sufficient for the present purpose.^{109–113} The relaxation induced by the plasmonic nanostructure can be obtained from the same calculation: at steady state, the energy dissipation rate can be obtained from the average work per unit time that this field does on the oscillating dipole. Furthermore, for a dipole oscillating in an environment characterized by a real dielectric response function, the spontaneous radiative relaxation rate can be described by the golden rule, $\gamma R=2\omega 3\u210f\epsilon 0|D|2\rho (r,\omega )$, in terms of the local density of states $\rho (r,\omega )=2\omega \pi c2ImTrG\u0332(r,r,\omega )$, providing an elegant framework for analyzing the effect of an inhomogeneous dielectric environment on the emission rate (Purcell effect).

In the electrostatic limit and point-dipole emitter model considered here, this formulation can be made significantly simpler and proceeds along the following steps:

- Following Gersten and Nitzan,
^{110}the point-dipole, positioned at**r**_{1}, is assumed to oscillate with a frequency*ω*and a constant amplitude*D*according to^{114}The Laplace equation is solved for the given dipole-nanostructure configuration using(2)$D(t)=De\u2212i\omega t+D*ei\omega t.$*ɛ*_{s}= 1 and the given*ɛ*_{ns}(*ω*) [Eq. (1)] with suitable parameters together with standard electrostatic boundary conditions at the surface(s) of the dielectric nanostructure to yield the electrostatic potential Φ(**r**,*ω*) and the electric field**E**(**r**,*ω*) =**E***(**r**, −*ω*) = −$\u2207$Φ(**r**,*ω*) anywhere in the system (excluding the point occupied by the dipole). The reaction field—the part**E**_{r}of the field**E**that arises from the polarization induced in the dielectric structure—is obtained as**E**_{r}(**r**,*ω*) =**E**(**r**,*ω*) −**E**_{d}(**r**,*ω*), where**E**_{d}is the field of the bare dipole. Note that these response functions depend on*ω*through the frequency dependence of*ɛ*_{ns}(*ω*) and that retardation effects are insignificant for the assumed small system. Also note that the linear response nature of standard electrostatics implies that these responses are linear in the driving dipole**D**. In particular, the component*E*_{r}(*ω*) of the field at the dipole position and in the dipole direction is used to define the needed component*G*(*ω*) of the Green’s function,This is essentially the electrostatic limit of the relevant component of the dyadic Green’s function mentioned above (calculated at the dipole position) except that the term associated with the bare dipole field is ignored.(3)$Er(\omega )=G(\omega )D.$ - This information is used to evaluate the dipole–metal nanostructure coupling by the following approximate procedure. The plasmonic response of the metal nanostructure is assumed to be describable as a linear combination of damped harmonic oscillators with frequencies fitted to the plasmonic peaks. In what follows, we assume that the dipole frequency is close to one of these peaks, of frequency
*ω*_{pl}, and ignore the contributions of other plasmons. Introducing a dimensionless coordinate*y*for the plasmonic oscillator, we write the dipole–plasmon interaction in the (assumed bilinear) formThe coordinate(4)$U=\u2212MyD.$*y*is proportional to the polarization induced in the plasmonic particle by the presence of the dipole, and*M*is the coupling parameter. The structure of Eq. (4) implies that*My*is the field associated with this polarization at the position and direction of the dipole, i.e.,In the quantum analog of Eq. (4), $y\u0302$ and $D\u0302$ are operators,(5)$My=GD.$*My*is the classical analog of $iM(a\u0302\u2212a\u0302\u2020)$, and in a vacuum system of volume*V*,*M*is given byUsing the operator equivalent of Eq. (4) in the form $U\u0302=\u2212My\u0302D\u0302$, the magnitude of the coupling matrix element between the lower molecular state(6)$M=\u210f\omega 2\epsilon 0V.$*g*with one photon and the upper molecular state*e*without a photonis a measure of the molecule–radiation field coupling. Our task is to evaluate(7)$\u3008U\u3009=\u3008g,1|U\u0302|e,0\u3009=MDge$*M*for a given molecule-nanostructure configuration. - To evaluate
*M*, we assume that the underlying dynamics associated with the plasmonic peak at*ω*_{pl}is that of a damped harmonic oscillator of a fictitious mass*m*and damping*γ*_{pl}(the plasmon width), which is driven by the dipole*D*via interaction (4),The long time solution of Eq. (8) under driving (2) is $y\u0304(\omega )e\u2212i\omega t+y\u0304(\u2212\omega )ei\omega t$ with(8)$my\u0304\u0308+m\gamma ply\u0304\u0307+m\omega pl2y\u0304=M\u0304D.$Transforming $y\u0304$ to a dimensionless coordinate(9)$y\u0304(\omega )=1mM\u0304D(\omega )\omega pl2\u2212\omega 2\u2212i\omega \gamma pl.$*y*according to $y(\omega )=(2m\omega pl/\u210f)1/2y\u0304(\omega )$ [thereby redefining $M\u0304$ as $M=(\u210f/2m\omega pl)1/2M\u0304$] and using Eq. (5), this leads (see Sec. S1.1 of the supplementary material) to a relationship between*G*(*ω*) and*M*,The field-dipole coupling parameter(10)$G(\omega )2=2\omega pl\u210f2M4(\omega pl2\u2212\omega 2)2+(\omega \gamma pl)2.$*M*and the parameters*ω*_{pl}and*γ*_{pl}are estimated by using the Laplace equation to obtain the reaction field component**E**_{r}(*ω*) at the position of a driving dipole**D**(*ω*) near the metal nanostructure and fitting $G(\omega )2=E(\omega )\u22c5D(\omega )/D(\omega )22$ to the form given by Eq. (10). In Sec. IV, we present the so-obtained coupling parameter in terms of an effective volume*V*_{eff}defined by(11)$M=\u210f\omega pl2\epsilon 0Veff.$ - The total dipole moment of the composite system
**D**_{tot}(*ω*) is calculated as**D**_{tot}(*ω*) =**D**(*ω*) +**D**_{ns}(*ω*), where**D**_{ns}(*ω*) is the dipole induced on the nanostructure, which is obtained as the following integral over the nanoparticle volume*V*_{ns}:(12)$Dns(\omega )=\epsilon ns(\omega )\u2212\epsilon s\epsilon 0\u222bVnsE(r,\omega )dr.$ - The radiative decay rate, modified by the presence of the nanostructure, is obtained from Refs. 110 and 115,Obviously, the ratio |(13)$\Gamma R(\omega )=\epsilon s\omega 33\u210fc34\pi \epsilon 0Dtot(\omega )2.$
**D**_{tot}|^{2}/|**D**|^{2}represents the Purcell lifetime modification factor for the present situation. - The non-radiative relaxation affected by the presence of the dielectric nanostructure is calculated as the rate of Ohmic heat generation on the nanostructure,
^{115}given by the following integral over the nanostructure volume:In addition to the non-radiative relaxation rate associated with dissipation in the dielectric nanostructure, an intrinsic molecular component, denoted as $\Gamma NR0$, results from the interaction of the molecule with its thermal environment and, for large molecules, also with intramolecular relaxation processes. This relaxation is molecule specific and cannot be quantitatively accounted for in our generic model. In the calculations reported below, unless otherwise stated, we have taken this intrinsic molecular relaxation rate to be equal to the radiative decay rate of the free molecule, implying that the emission yield of the free molecule is 0.5.(14)$\Gamma NRns(\omega )=\epsilon 0Im[\epsilon ns(\omega )]2\u210f\u222bVnsE(r,\omega )2dr.$ - The emission quantum yield of the molecular emitter–plasmonic nanostructure system is given by(15)$Y(\omega )=\Gamma R(\omega )\Gamma R(\omega )+\Gamma NRns(\omega )+\Gamma NR0(\omega ).$
- The strong coupling condition on this level of treatment is given bywith ⟨(16)$2\u27e8U\u27e9>\u210f\Gamma =\u210f\Gamma R+\Gamma NRns+\Gamma NR0$
*U*⟩ given by Eq. (7).

The above calculation scheme provides a rather straightforward, albeit approximate, procedure for mapping strong coupling regions in a system of coupled nanostructure and molecular emitter, where the driving frequency *ω* is naturally taken to be the relevant molecular transition frequency *ω*_{0}. It suffers, however, from a drawback, originating from the fact that since the dipole emitter is assumed to drive the system, such calculation cannot account for any spectral shift that is another manifestation of the molecule–radiation field coupling. The criterion (16), in which the coupling and relaxation rates are calculated for the driving frequency *ω*, may thus be viewed as an approximate criterion for the onset of strong coupling, but may not fully account for the spectral structure and dynamics beyond it.

To overcome the latter limitation, the calculation should be done self-consistently (see, e.g., Refs. 117–118). Here, the system, comprising the molecular emitter and the nanostructure, is driven by an external oscillating electric field, **E**(*t*) = **E**(*ω*) cos(*ωt*) = **E**(*ω*)Re *e*^{−iωt} with a real amplitude **E**(*ω*), which couples to the molecular emitter. The dynamics of the latter can be represented as a driven damped harmonic oscillator, which in Fourier space takes the form

where ** α** is the static polarizability tensor and $G\u0332(r,\omega )$ is defined in Eq. (3). Here,

*ω*

_{0}is the transition frequency of the free emitter, while Γ

^{0}represents the combined radiative, $\Gamma R0$, and non-radiative, $\Gamma NR0$, relaxation rates associated with this free emitter. For simplicity, we assume that the molecular polarizability tensor is diagonal

**=**

*α**α*

_{x}

**e**

_{x}+

*α*

_{y}

**e**

_{y}+

*α*

_{z}

**e**

_{z}with

**e**

_{j},

*j*=

*x*,

*y*,

*z*being unit vectors. Equation (17) can then be written for each component of

**D**in the form

Here, the static polarizability of the two-level molecular model is given by

where $Dge(j)$ is the *j*th component of the molecular transition matrix element **D**_{ge}. Making another simplifying assumption that $G\u0332$ is diagonal (this will be true for geometries considered in our calculations below) leads to

The absorption line shape *A*(*ω*) may be calculated as the work done on the molecular dipole per unit time, averaged over the driving period, $Ej(t)D\u0307j(t)\u0304$. This yields (see Sec. S1.2 of the supplementary material)

The *αG* term in Eq. (20) represents the effect of the nanostructure on the dynamics of the molecular emitter. Its real and imaginary parts lead, respectively, to a spectral shift and to modifications in the radiative and nonradiative relaxation rates. When *G*(*ω*) is characterized by a resonance structure, e.g., a distinct plasmonic peak as seen in Eq. (10),^{119} a split peak structure may appear in *A*(*ω*) and will indicate the onset of strong coupling. The following points should be noted:

The classically calculated absorption line shape, Eq. (21), can be used to calculate the optical response of the molecular dipole-nanostructure composite as a function of frequency, thereby directly looking at the spectral manifestation of strong coupling, not just at the approximate condition Eq. (16) for its onset. The needed response function

*G*is obtained as a function of*ω*from the solution of the Laplace equation as described above Eq. (3). The same calculation also yields the radiative and nonradiative decay rates and the quantum yield as functions of*ω*, as described above.In setting up the model of Eq. (17), we aimed to examine the effect of proximity of the molecular emitter to a plasmonic nanostructure. We have therefore implemented the driving by an external field as acting only on the molecular emitter. In reality, both the molecule(s) and nanostructure comprise the system that responds to the incident field, and the effect of directly driving the nanostructure by this field should be included. In this case,

**E**(*ω*) on the RHS of Eq. (17) will be replaced by the vector sum of the incident field and the field created by the polarization induced on the nanostructure by the incident field—an important ingredient in the electromagnetic theory of surface enhanced spectroscopies. Consequently, after making the same simplifying assumptions as above,*E*_{j}(*ω*) in Eq. (20) will be augmented by the corresponding additive field from that polarization. This will not change the consequences regarding strong coupling as determined by the denominator in Eq. (20).The simple criterion [Eq. (16)] for strong coupling, namely, the appearance of split peaks at the crossing of the molecular transition and a plasmonic resonance frequency, is derived from a model in which the plasmon is represented, like the molecule in Eq. (17), as a damped harmonic oscillator with bilinear coupling between these oscillators.

^{105}This is a reasonable approximation where the distance between the molecular emitter and the plasmonic structure is large enough to represent the response of the latter by its dipolar plasmon alone. In the general case, while the condition for strong coupling is again contained in the properties of the denominator in Eq. (20), its actual manifestation depends on the properties of the response function*G*and can be determined only numerically. Note that*G*accounts not only for plasmonic effects but also for other electrostatic aspects of the response, such as the lightning rod effect.^{120}

Whichever level of analysis we aim at, a prerequisite of the calculation is the evaluation of the response function *G* and the radiative and non-radiative relaxation rates of the molecular emitter when oscillating in close proximity to the given nanostructure. The models depicted in Figs. 1(a) and 1(b) admit analytical solutions for these functions. These are described in Sec. III and in the supplementary material. Using these solutions, we will analyze in Sec. IV, for silver and gold particles (as represented by their dielectric response function) in some select configurations, the emergence of strong coupling and the molecular emission properties as functions of geometrical parameters.

## III. RESPONSE FUNCTION, RELAXATION RATES, AND EMISSION YIELD

In this section, we present analytical results for the response function *G* and the radiative and non-radiative relaxation rates for several geometries associated with the structures displayed in Figs. 1(a) and 1(b) (the case of a molecule near a spherical particle is obviously a limit of the spheroid results presented below). The solutions of the Laplace equations are obtained as infinite sums that are computed numerically by truncating the series while ensuring convergence beyond a desired accuracy. For quick estimates, we focus on the procedure described in Eqs. (3)–(16). In Sec. III C, we provide some details on the calculation based on Eq. (21).

### A. A molecular (point) dipole near a prolate spheroidal nanoparticle

The solution of the Laplace equation for this geometry is best done in the prolate spheroidal coordinate system (supplementary material, Sec. S2.1). In this coordinate system, the nanoparticle surface is defined by *ξ* = *ξ*_{0}. The potentials inside and outside the spheroidal particle (solutions of the Laplace equation) are written in the form

$Plm(\xi )$ and $Qlm(\xi )$ are the associated Legendre functions of the first and second kinds, respectively. The two foci of the nanoparticle overlap with that of the coordinate system itself at $\xb1a2$. Here and below, the dielectric response functions of the nanoparticle and the environment are taken as *ɛ*_{d}(*ω*) and *ɛ*_{s} (assumed frequency independent), respectively. The dipole **D** is situated at **r**_{1}(*ξ*_{1}, *η*_{1}, 0) (implying that it is on the *xz* plane) with an arbitrary orientation relative to the nanoparticle surface, $D=D\xi \xi \u0302+D\eta \eta \u0302$, where $\xi \u0302$ and $\eta \u0302$ are the unit vectors along the coordinates *ξ*, *η*, respectively. Following Ref. 110, we impose continuity of the potential Φ and of the component of the displacement field normal to the spheroidal surface, $\epsilon \u2202\Phi (\xi ,\eta ,\varphi )\u2202\xi $, on this surface. To this end, the dipole potential [last term of Eq. (23)] is written in prolate spheroidal coordinates [Eq. (S12) in the supplementary material] and imposing these boundary conditions leads to the coefficients *A*_{l,m} and *B*_{l,m} as follows:

and

where

and

The solution Φ^{out}(*ξ*, *η*, *ϕ*) of Eq. (23) is related to the response dyadic $G\u0332$ according to $\u2212\u2207\Phi (r)=G\u0332\u22c5D$. In particular, the component of $G\u0332$ in the dipole direction is evaluated here. To exemplify, for the configuration where the dipole is placed along the major axis of the nanoparticle and is oriented perpendicular to the nanoparticle surface, i.e., *η*_{1} = 1, $D=D\xi \xi \u0302$, the component $Gsphdmaj,\u22a5$ is obtained in the following form (maj ≡ *major*, sphd ≡ *spheroid*):

For the same configuration, the non-radiative decay rate associated with heat production in the spheroid is calculated from Eq. (14), which leads to

In addition, from Eq. (12), we calculate the total dipole in the molecule–dielectric spheroid system. The result is

This can be used with Eq. (13) to calculate the radiative decay rate. Finally, the emission yield is calculated from Eq. (15). The detailed derivation of Eqs. (24)–(30) is provided in the supplementary material. Similar results for several other configurations are also detailed in the supplementary material.

### B. A molecular (point) dipole between two nanospheres

As in Ref. 121, the Laplace equation for this problem is solved in a bispherical coordinate system (*μ*, *η*, *ϕ*) (supplementary material, Sec. S4.1). The nanospheres of radii *r*_{R} and *r*_{L} (their surfaces defined by *μ*_{R} and *μ*_{L}, respectively) are characterized by the dielectric response functions *ɛ*_{R}(*ω*) and *ɛ*_{L}(*ω*) with their centers positioned on the *z* axis at *z* = *a* coth *μ*_{i}, *i* = L,R, *x* = *y* = 0 (in Cartesian coordinates) where the poles of the coordinate system are at *z* = ±*a* on the *z* axis. The dipole is also positioned on this line at (*μ*_{1}, *π*, 0) (*x*_{1} = *y*_{1} = 0 in Cartesian coordinates).

#### 1. Parallel configuration

For a dipole oriented along the line connecting the sphere centers, $D=D\mu \u0302$, where $\mu \u0302$ is the unit vector corresponding to the coordinate *μ*. The azimuthal symmetry of this configuration simplifies the solution. The solution of the Laplace equation for *μ* > *μ*_{R} > 0, namely, inside the right sphere, is written in the form

For *μ* < *μ*_{L} < 0, namely, inside the left sphere, we have

The potential outside both the nanoparticles takes the form

where the dipole potential is written explicitly [last term in Eq. (33)]. The boundary conditions, given by

yield an infinite linear system of equations for the coefficients (see Sec. S4 of the supplementary material) that can be solved numerically after truncating at a desired order *n*_{max} to give the needed coefficients up this order. Expressions for the component of $G\u0332$ in the dipole direction and the radiative and nonradiative rates in terms of these coefficients are provided in Sec. S4.2 of the supplementary material.

#### 2. Perpendicular configuration

As before, the molecule is positioned at (*μ*_{1}, *π*, 0) on the line connecting centers of the two nanoparticles with parameters defined as in the parallel configuration case. The environment is described by a frequency independent dielectric constant *ɛ*_{s} as usual. The molecular transition dipole is taken to be oriented perpendicular to the line connecting the sphere centers, given by $D=D\eta \u0302$, where $\eta \u0302$ is the unit vector of the coordinate *η*. As this configuration no longer has an azimuthal symmetry, the solutions of the Laplace equation for *μ* > *μ*_{R} > 0 (inside the right sphere) and *μ* < *μ*_{L} < 0 (inside the left sphere) are written as

and

respectively. The potential outside both the nanoparticles is given by

The same boundary conditions given by Eqs. (34)–(37) yield an infinite linear system of equations for the coefficients (supplementary material, Sec. S4.3) and are solved numerically as done for the parallel configuration. The radiative and nonradiative rates and the relevant component of $G\u0332$ in terms of these coefficients are provided in Sec. S4.3 of the supplementary material.

### C. The self-consistent procedure

For the configurations discussed above, once the needed component of the Green’s dyadic has been evaluated, we can use Eqs. (21) and (19) and the given free molecule parameters Γ^{0} and *ω*_{0} to calculate the absorption line shape for the given dipole-metal nanostructure configuration. We emphasize again that Eq. (21) is an expression for the absorption line shape for a model in which the incident light interacts with the molecular dipole only. In such a model, a distinct plasmon-dominated peak accompanying the molecular response is a manifestation of “strong coupling.” Generalization to the case where the incident radiation drives both the molecular dipole and the metal nanostructure is straightforward, but more costly, and the resulting line shape will also show the signature of interference between the dipole and the metal responses.

## IV. RESULTS AND DISCUSSIONS

Here, we present numerical results that map single molecule strong coupling regions about spherical, spheroidal, and bispherical structures. Three comments are in order. First, because strong coupling is defined by Eq. (16) or by the appearance of split peaks at the overlap of the molecular and plasmon resonances as calculated from Eq. (21), it necessarily depends on properties of the molecular resonance (non-radiative relaxation aside from that due to coupling to the metal nanostructure as well pure dephasing) that are not addressed in the present analysis. The results shown below are obtained by assuming that the quantum yield for emission by a free molecule is 0.5 and that pure dephasing is absent.

Second, our classical calculations disregard quantum mechanical effects, mainly tunneling and non-locality of the dielectric response, which are expected to become important at small molecule–metal distances.^{105} Classical estimates should be regarded reliable only at distances larger than, say, 0.3 nm where both effects are small.

Finally, as discussed in Sec. II, estimates of strong coupling based on Eq. (16) are made at a given frequency, usually taken (as in Figs. 2–7) as the transition frequency *ω*_{0} of the free molecule. The metal-induced spectral shift, satisfying the inequality (16) at this frequency, while indicating strong coupling, does not therefore guarantee the appearance of a split peak in the actual plasmon–molecule system. The actual spectral shape is obtained from Eq. (21), as demonstrated in Figs. 8 and 9.

Figures 2–7 show the coupling strength ⟨*U*⟩, the lifetime broadening Γ, and the coupling-broadening ratio (CBR) 2|⟨*U*⟩|/*ℏ*Γ calculated at the resonance frequency of the free molecule for several geometries. A value larger than 1 of the CBR indicates strong coupling by the criterion (16). For definiteness, in the calculation shown below, the molecular transition dipole is taken to be 10 D. However, the free molecule transition frequency *ω*_{0} is taken to be equal to the dipolar plasmon frequency of the metal corresponding to geometry under consideration. For example, when considering a molecular dipole near small gold and silver spheres, the molecular frequencies are taken equal to corresponding dipolar plasmon frequencies, 2.62 and 3.37 eV for gold and silver spheres, respectively. However, the coupling strength, which is calculated by fitting the calculated *G* to the plasmon peak, is, by this definition, frequency independent.

The calculations of the radiative and metal-induced nonradiative contributions to Γ [Eqs. (13) and (14)] have been extensively discussed before,^{110} and typical results showing the dependence of these rates on the molecular orientation and its distance from the metal structure and on the size and shape of this structure are shown in the supplementary material (Sec. S5). We emphasize again that the only contributions to the linewidth Γ that are considered here are the lifetime (radiative and non-radiative) of the free molecule and the metal induced modifications of these lifetimes and that in the calculations reported below, we have assumed that a free molecule has an intrinsic non-radiative relaxation rate that is equal to its radiative rate (implying an emission yield of 0.5 for the free molecule). Note that if the intrinsic non-radiative relaxation rate of the free molecule is disregarded, the remaining contributions to the linewidth are proportional to the square of the molecular transition dipole *D*_{ge}, whereas the coupling ⟨*U*⟩ [Eq. (7)] is linear in this parameter. Consequently, and perhaps counter-intuitively, the CBR satisfies 2|⟨*U*⟩|/*ℏ*Γ ∝ 1/*D*_{ge} and decreases with the increasing molecular transition dipole. The relaxation rates (and magnitude of coupling) are normalized by the radiative relaxation rate $\Gamma R0$ of a free molecule, which is calculated from $\epsilon s\omega 03Dge2/(3\u210fc34\pi \epsilon 0)$ at the respective molecular frequency.

Figure 2 depicts the CBR for a molecule near spherical gold and silver particles of radius 10 nm, in normal and parallel orientations relative to the sphere surface, as a function of the molecule–surface distance. Also shown is the combined relaxation rate Γ. It is seen that “strong coupling” as defined by the CBR > 1 criterion is not necessarily associated with close proximity. Both and ⟨*U*⟩ and Γ go quickly to zero when the molecule–surface distance increases; however, their ratio varies relatively slowly and remains larger than 1 in all the distance range studied. It is interesting to note that in both the normal and parallel configurations, this ratio goes through maximum for silver and decreases beyond a distance $\u223c1$ nm, while for gold, it rises monotonically with increasing distance (results for an extended distance range are shown in Fig. S5.1 of the supplementary material, but should be regarded cautiously because of the use of the electrostatic approximation). The metal induced nonradiative relaxation, being the major contribution to the combined relaxation rates, is very large when close to the nanoparticle surface. As a consequence, for silver, there is an optimum distance of the molecule from the surface where the CBR is maximum. For gold, on the other hand, because the extent of the effect of nonradiative decay rate is even more, for the distance considered, CBR is larger as the molecule is distant from the surface. Qualitatively, for both perpendicular and parallel orientation, these observations are common.

These estimates of strong coupling show considerable dependence on the size and shape of the dielectric nanostructure. Figure 3 shows the CBR parameter and the lifetime broadening Γ plotted against the molecule–(spherical gold) nanoparticle surface distance for different particle sizes in the perpendicular dipole orientation, while Figs. 4 and 5 show similar results for different particle shapes by considering metal spheroids of constant volume (= volume of a sphere of radius 10 nm) and different aspect ratios *p*/*q* when the molecule is positioned on the long axis with a normal orientation relative to the nanoparticle surface. In Figs. 4 and 5, we show the magnitude of coupling ⟨*U*⟩ (panel b of Figs. 4 and 5) and total relaxation rates Γ (panel c of Figs. 4 and 5). More details on the radiative and non-radiative contributions to Γ for a molecular dipole positioned near gold and silver nanospheroids are provided in Sec. S5.2 of the supplementary material.

Figures 6 and 7 show similar quantities for a molecular dipole positioned between two metal spheres. The (point) dipole is placed on the axis connecting the sphere centers (called the intersphere axis) and oriented parallel (and normal) to this axis. The molecular transition frequency is taken to be in resonance with the stronger plasmon peak associated with the corresponding configuration (see Fig. S8 for the optical response of gold and silver bispherical structures).

Next, we consider the effect of the gap size for the likes of the configurations shown in Fig. 6, i.e., molecule in the gap in between nanosphere dimers, specifically for the parallel molecular orientation relative to the intersphere axis. In particular, in Fig. 7, the CBR parameter, the coupling magnitude, and the combined relaxation rate are displayed for both gold and silver nanoparticle dimers. In all calculations, the molecule is assumed to be positioned in the middle of the gap. To note, the molecular frequency *ω*_{0} is taken to be in resonance with the higher plasmon peak for the configuration with two 10 nm metal nanoparticles with a 2 nm gap. For this reason, the calculated combined relaxation rate Γ in Fig. 7(c) shows a clear maximum for 2 nm gap size for silver. This also explains the slight deviation from the trend for gold.

The configuration considered here, with the molecule positioned in the space between two nanoparticles, is often referred to as a plasmonic cavity. The effective cavity volume can be determined from the calculated M using Eq. (11). Obviously, this “effective volume” can also be defined for configurations that cannot be perceived as cavities, such as used in Figs. 2–5 (dipole–sphere/spheroid systems), and it depends on the particular geometry considered, including the dipole position and orientation. It is therefore of limited value for molecules interacting with plasmonic structures, but because of its prominence in many discussions, we show several examples of this parameter in the supplementary material (Figs. S6 and S7).

The interplay between coupling and broadening that determine “strong coupling” according to Eq. (16) strongly depends on geometry. While this fact is widely appreciated, the particulars of these dependence are sometimes surprising as exemplified in points (C) and (D) below.

In the plasmonic structures studied, the intrinsic coupling-broadening ratio (CBR) parameter is usually large and by itself satisfies the strong coupling criterion. Here, “intrinsic” implies that we consider only spectral broadening effects that reflect (a) the molecular radiative lifetime and fluorescence yield associated with the transition considered and (b) the effect of the electromagnetic molecule–metal interaction on the molecule radiative and non-radiative relaxation rates. This intrinsic broadening disregards other sources of spectral broadening, such as congestion of many overlapping transitions and environmentally induced thermal relaxation.

While both the molecule–plasmon coupling and the intrinsic molecular lifetime broadening (by “intrinsic,” we mean broadening associated only with the isolated molecule, the plasmonic nanostructure, and their interaction) strongly decrease with the increasing molecule–metal surface distance,

^{122}the CBR remains large and the corresponding criterion for “strong coupling” persists at large distances. This observation by itself is meaningless as other sources of broadening and limited resolution will usually mask the Rabi structure. This does indicates, however, that if other sources of broadening and relaxation are eliminated (such as with zero-phonon transitions at cryogenic temperatures), signature of strong coupling would be observed in single molecule plasmonic cavities, in agreement with a recent observation.^{102}Because of the interplay between local electromagnetic coupling and particle induced relaxation, “hotspots” characterized by a large enhancement of the local electromagnetic fields do not necessarily stand out with regard to the CBR. To illustrate this, consider the perpendicular and parallel orientation of a molecular dipole of 10 D transition dipole moment positioned 1 nm apart from a silver nanosphere of 10 nm radius. For the perpendicular orientation, the coupling strength is larger in magnitude than the parallel orientation: −⟨

*U*⟩_{⊥}/*ℏ*= 1.272 × 10^{14}s^{−1}(perpendicular) and −⟨*U*⟩_{∥}/*ℏ*= 8.935 × 10^{13}s^{−1}(parallel). However, the associated nonradiative relaxation rate is also larger for perpendicular orientation: $\Gamma NR\u22a5=6.437\xd71012s\u22121$ (perpendicular) and $\Gamma NR\u2225=2.091\xd71011s\u22121$ (parallel) (the radiative decay rate is smaller than these number by $\u223c1$ order of magnitude at this small distance from metal surface, so does not have significant contribution to the CBR). As a result, the value of CBR for parallel orientation is more than twice the value for perpendicular case (CBR^{⊥}= 36.733, CBR^{∥}= 80.744).

Finally, consider the self-consistent treatment of Eq. (21). We reiterate that the results shown in Figs. 2–7 are based on a calculation of ⟨*U*⟩, Γ_{R}, Γ_{NR} computed in a model in which the molecular dipole drives the system and do not reflect the actual line shape. In particular, such a calculation does not account for spectral line shifts induced by the proximity of the molecule to the metal nanostructure and for the relative peaks intensities, even if Rabi splitting occurs. Therefore, the results shown in Figs. 2–7 only indicate the fulfillment of the strong coupling criterion [Eq. (16)] but not the observability of an actual Rabi splitting. Figures 8 and 9 show results based on the self-consistent calculation, Eqs. (17)–(21), which yields the actual absorption line shape for the process in which a dipolar emitter, coupled to the metal nanostructure, is driven by an incident external field.

Figure 8 shows this line shape, calculated for a dipole emitter near a spherical gold nanosphere (radius 10 nm), oriented perpendicular to the sphere surface and parallel to the incident field. The transition frequency *ω*_{0} is taken 2.62 eV, in resonance with the dipolar plasmon frequency of the gold nanosphere. As before, the relaxation rate Γ^{0} of the free molecule is taken to be twice the radiative relaxation rate, corresponding to an emission yield of 0.5 for the free molecule. Figure 8(a) depicts the absorption line shape for a dipole situated at 0.5 nm away from the sphere surface for different values of the molecular transition dipole moment in the range 2–10 D. Figure 8(b) displays the absorption line shape for an emitter with transition dipole moment 10 D placed at different distances (0.5–5 nm) from the sphere surface. Similar results for silver are shown in Fig. S5.6. Figure 9 shows similar results for a molecular dipole in the middle of the gap between two identical gold spheres (radii 10 nm) oriented parallel to the intersphere axis, where the incident electric field is parallel to this axis as well. For a bispherical structure, the plasmonic response for an incident plane field is characterized by a double-peak structure, which results from the interacting dipolar plasmonic responses of the individual spheres and depends on the inter-sphere distance. The emitter transition frequency is taken to be in resonance with the lower frequency and more intense one of these peaks (see Fig. S8a). The results shown in Fig. 9(a) are for an intersphere gap size 1 nm and different values of the molecular transition dipole moment, while Fig. 9(b) shows the line shapes for a transition dipole 10 D at varying gap sizes.

Consider first the case of an emitter near a spherical gold nanoparticle (Fig. 8). Keeping in mind that in this model calculation the incident field is taken to be coupled only to the emitter, it is the emitter absorption peak that is seen for weak coupling [orange line in Fig. 8(a) and black line in Fig. 8(b)]. The main effect of strong coupling upon increasing the molecular transition dipole is seen to be a strong red shift of the molecule-dominated peak. The particle’s optical signature, resulting in the molecule–particle interaction that is seen on the high energy side of the spectrum, is small, leading us to conclude that a pronounced Rabi splitting is not expected in this configuration. In contrast, the optical response of a system comprising a molecule positioned between two gold spheres does show a clear evidence of such splitting, as seen in Fig. 9. This observation is consistent with that of Ref. 123 that for a dimer-like structure of gold nanoparticles; a strong coupling situation is likely to result when the interparticle gap is less than 2 nm.

While “strong coupling” between a single molecule and a metal-plasmon excitation, manifested in the split-peak structure of the absorption line shape, is seen in the calculated spectra in Fig. 9, these results also indicate that the simpler model calculation of the CBR shown in Figs. 2–7 is not by itself sufficient to predict the observability of this spectral feature. The CBR, Eq. (16), does provide an estimate of the relationship between the molecule–metal nanostructure coupling and the associated broadening but does not necessarily lead to an observed Rabi splitting, which is usually assumed to be a prime manifestation of such strong coupling. We have identified the strong spectral shift, another manifestation of this coupling, as the main reason for this apparent discrepancy. This may also lead to apparent counter-intuitive observations: for example, comparing the line shapes in Fig. 9(a), a clear Rabi splitting is seen for an emitter with transition dipole 4 D, while a single peak is seen for a 2 D emitter. In the stronger coupling case of an emitter with transition dipole moment $\u22656$ D, a larger Rabi splitting is indicated by the calculation, but the dominance of the lower energy peak may appear as a single peak in a realistic observations. It should be noted that this shift will also affect the distance dependence of the radiative and radiationless relaxation rates of a molecule approaching a metal surface and may potentially play a role in the observed “quenching of quench,” where the emission yield of an emission appears to increase, rather than decrease, as the molecule–surface distance decreases (see Ref. 100 and references therein). We emphasize however that no such trend in the emission yield *Y*, Eq. (15), was seen in our calculations.

## V. CONCLUSIONS

Strong radiation matter coupling is often characterized in the literature as a relative concept, by comparing the absolute magnitudes of the coupling ⟨*U*⟩ and the broadening Γ, and the strong coupling criterion 2|⟨*U*⟩|/*ℏ*Γ, an indication that Rabi splitting may be observed when the molecular and plasmon optical transition come into resonance, is often used. While usually observed when an optical mode interacts with many molecules, it has been noted, as outlined in the Introduction, that this criterion can be satisfied even for single molecules interacting with metallic plasmonic structures. A necessary condition for observing this phenomenon in such a setting is that the coupling-broadening ratio (CBR), 2|⟨*U*⟩|/*ℏ*Γ, satisfies the strong coupling criterion when Γ and ⟨*U*⟩ are derived from the interaction between a single molecular transition when the molecule interacts only with the plasmonic structure, in the absence of other sources of broadening. Using simple analytically solvable model structures, representing the molecule as a point dipole and treating radiation–matter coupling in the long wavelength approximation, whereas all distances are assumed small relative to the radiation wavelength, we have found that this condition is easily satisfied. Surprisingly, we observe that this remains true when the molecule–metal distance increases so that the absolute coupling itself becomes small. This suggest that Rabi splitting can be observed in systems involving a single molecule interacting with small metal particles if other sources of broadening can be suppressed, as was recently observed using the low temperature molecular zero-phonon transition.^{102}

It should be kept in mind that even if the strong coupling criterion, Eq. (16), is satisfied, such an observation by itself is not an indication of strong radiation–matter coupling when this coupling is compared to other molecular energetic parameters. Furthermore, by considering the results of a self-consistent calculation that addresses directly the absorption line shape in a system of interacting a molecule and the plasmonic structure, we find that the intrinsic CBR satisfying the strong coupling criterion does not guarantee the observation of Rabi splitting in the absorption spectrum because the coupling induced spectral shift together with the strong frequency dependence of the imaginary part of plasmonic response can result in a spectrum dominated by a single, albeit shifted, molecular peak.

## SUPPLEMENTARY MATERIAL

See the supplementary material for short derivations of Eqs. (10) and (21), details of the calculations described in Sec. III, and results for some additional configurations similar to the ones described in Sec. IV.

## ACKNOWLEDGMENTS

M.S. acknowledges the financial support from the Air Force Office of Scientific Research under Grant No. FA9550-19-1-0009. A.N. acknowledges the support of the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available within the article and its supplementary material.

## REFERENCES

_{2}O

_{6}nanosheets

_{2}reduction on Au nanoparticles

_{2}O nanocages

We work in the dipole gauge and disregard the dipole self-interaction term in the Hamiltonian. While strong coupling [by the criterion of Eq. (16)] is found in our systems, we never encounter situations of ultra-strong coupling that will require taking this term into account.

In a more general formulation, this would be the Fourier component **D**(*ω*)*e*^{−iωt} + **D**(−*ω*)*e*^{iωt} of a general time dependent dipole **D**(*t*). In accordance, Eq. (3) will take the form *E*_{r}(*ω*) = *G*(*ω*)*D*(*ω*).

Such a peak appears when a pole of the response function is characterized by an imaginary part that is smaller than its distance from other poles. Such behavior is manifested by the dipolar plasmon of some metals and semiconductors. In the electrostatic limit under consideration, these poles depend on the shape, but not the size of the dielectric nanostructure. An example is shown in the supplementary material, Sec. S5.7, for the dipolar plasmons of a prolate spheroid.

For the metal–molecule distances considered in our calculations, the intrinsic lifetime broadening is dominated by the molecule–metal interaction.