Plasmonic quantum effects on single-emitter strong coupling

2.1 Theoretical framework

In a general hydrodynamic description, the linear nonlocal relation of the polarization vector P(r, ω) and the electric field E(r, ω) is described by the equation [47], [48], [49], [50]:

where e and m are, respectively, the magnitude of the charge and the mass of the electron, n(r, ω)=n0(r)+1e∇⋅P is the electron density with n₀ being the equilibrium charge density, γ is the phenomenological damping rate and ω_p(n₀) is the plasma frequency. The first term in Eq. (1) takes into account electron-electron interactions and represents the quantum pressure. In particular, g is the first-order term (with respect to perturbations of n) of the functional derivative of the energy functional G, i.e. g=(δG[n]δn)1. The complexity of the energy functional sets the level of approximation of the QHT. In this article we will take G[n] as the sum of the TF kinetic energy functional, the von Weizsäcker (vW) correction and the exchange-correlation (XC) energy term in the limit of the local density approximation, i.e. G[n, ∇n]=T_TF[n]+T_vW[n, ∇n]+E_XC[n]. The second term in Eq. (1) has the form of a viscoelastic tensor and takes into account the nonlocal damping [58], [59], which is especially relevant for geometries with small features [50]. Explicit expressions for these functionals can be found in Refs. [49] and [50].

If we set both g=σ≡0, Eq. (1) reduces to the usual Drude model in which the relation between E and P is purely local. We will refer to this approximation as LRA. In the TF hydrodynamic theory (TF-HT) only the TF kinetic energy functional is considered, i.e. G[n]=T_TF. In this approximation we will also set σ=0 and n₀=constant in the metal volume throughout the article. The principal characteristics for each approach are summarized in Figure 1.

Coupled to the continuity equation and Maxwell’s equations, Eq. (1) describes the optical response of metals. In its general form it takes into account nonlocal electron response, electron charge spill-out, size-dependent absorption and retardation effects. We have implemented this system of equations in a finite-element method (FEM) commercial platform using a 2.5D scheme [49], [50], [60] for axis symmetric structures. The single-molecule emission is approximated by a point-like dipole field and its interaction with the plasmonic system is implemented using a scattering field formulation that avoids numerical complications due to the singularity of the dipole field (see the Methods section for more details).

In order to study the strong interaction of the emitter with the metallic system described by Eq. (1) we employ the method developed by Alpeggiani et al. [61], which allows the study of non-perturbative effects by evaluating the perturbative photon emission decay rate. In this framework, the power spectrum (averaged over the solid angle) emitted to the far-field by a point-like dipole in proximity to a metal nanostructure is given by [61], [62]:

where the dipole spectrum S′ is defined as:

where ω and ω₀ are, respectively, the emitted frequency and the dipole transition frequency, γ_tot=γ_r+γ_nr is the perturbative dipole decay rate (including radiative and non-radiative contributions), δω the photonic (anomalous) Lamb shift, and q=γ_r/γ_tot is the quantum yield of the system, which gives the probability that a transition of the emitter results in photon emission to the far-field as opposed to being absorbed by the nearby environment.

S(ω) is the energy radiated by the dipole at frequency ω. Since no external excitation is present in our system, S(ω) represents the number of photons of energy ℏω emitted in the far-field per unit of time in all directions. Evaluating this quantity corresponds to measuring luminescence in experimental setups. The quantity S′(ω), by contrast, represents the total number of photons emitted by the molecule, including those that are lost (absorbed) in the plasmonic environment and do not reach the far-field.

All the quantities introduced can be calculated from the electromagnetic fields as detailed in the Methods section.

2.2 Spheres

Before delving into the strong coupling regime calculations, it is interesting to observe the impact of nonlocal and quantum effects on near-field and far-field optical properties by comparing the electromagnetic local density of states (LDOS) ρ(r₀, ω) and the extinction spectra, respectively. In the following, we normalize the extinction cross-section to the geometric cross-section (πR² for spheres and 2πR² for sphere dimers), which yields the extinction efficiency q_ext. Similarly, we normalize the LDOS to the free-space LDOS, ρ₀(ω)=ω²/(π²c³). Moreover, since the dipole parallel component is quenched at very small distances from metal surfaces [63], we only consider a dipole normally oriented with respect to the close-by metal surface.

Figure 2 depicts the extinction efficiency and the LDOS enhancement spectra for Drude-like Na spheres of radius R varying from 2.5 nm to 20 nm. The LDOS refers to the states that can be excited by a point-like dipole placed at d=1 nm from the nanoparticle, normally oriented (as shown in the inset of Figure 2B). The QHT results are compared against the LRA and the TF-HT calculations. It is interesting to notice that for the smallest particle (R=2.5 nm) both the extinction and the LDOS are strongly affected by the model used to describe the metal particle. The LRA extinction shows the usual dipole resonance at ω=ωp/3, which appears shifted toward the blue in the TF-HT spectra and toward the red in the QHT spectra [48]. The TF-HT shows a set of resonances above the plasma frequency due to longitudinal modes [37]. These modes are not visible in the QHT spectra, due to additional nonlocal losses introduced by the viscoelastic tensor [50], which visibly broadens the dipole resonance. As the particle size increases, the extinction spectra for the different models converge to the same curve, as would be expected. Surprisingly, this is not the case for the LDOS spectra, which remain neatly different even for the largest nanoparticle (R=20 nm). The LRA spectra show a set of new resonances in addition to the dipole resonance that tend to sum up into a large super-resonance peak around ω=ωp/2 for bigger particles. Note that this is the surface plasmon resonance condition for a flat surface. The limit of a dipole infinitesimally close to the sphere is, in fact, equivalent to a dipole nearby an infinitely large sphere. The TF-HT, however, displays resonances above the LRA limit ωp/2, as already pointed out by Christensen and co-authors [54]. More interesting is the fact that as the particle size is increased, the TF-HT spectrum seems to build up similarly to the LRA case, although this time into a much broader resonance above ωp/2. For the QHT spectra, with respect to the extinction, everything is blurred by the presence of nonlocal losses, which are always large due to the fact that the dipole is extremely close to the metal surfaces. Although broad, a resonance can be seen converging below ωp/2 for larger distances. Similar results have been recently observed by Gonçalves et al. [64] using the Feibelman d-parameters approach, in which the surface effects are embedded in an effective frequency-dependent parameter. The QHT, however, goes beyond purely surface effects and offers a full nonlocal (k-dependent) description of the free-electron dynamics, which is particularly relevant when the emitter is in very close proximity to the metal surface. In the QHT spectra, a new resonance can be seen right below the bulk plasma frequency ω_p, persisting even for the largest particle. It is, however, hard to identify the nature of the mode in question since the field maps are mostly characterized by the strong evanescent fields of the dipole.

Figure 2:

Far-field versus near-field properties for metallic spheres of different radii R, for the local response approximation (LRA), Thomas-Fermi hydrodynamic theory (TF-HT) and quantum hydrodynamic theory (QHT) models, respectively.

(A) The extinction efficiency spectra. (B) The normalized local density of states (LDOS) for states excited by a normally oriented dipole placed at a fixed distance d=1 nm. In our calculations we assumed a Drude-like metal with bulk parameters ℏω_p=5.89 eV and ℏγ=0.066 eV corresponding to sodium (Na).

Next, we fix the particle radius at R=20 nm and we vary the distance d of the emitter from the nanoparticle. Figure 3 shows the spontaneous (total) and radiative emission rate enhancement, γ_tot/γ₀ and γ_r/γ₀, respectively, for all considered models. As expected, at large distances (d=16 nm), all models converge to an identical spectrum, in which the dipole resonance is clearly distinguishable from higher order modes. As the distance is reduced, we obtain, for the total rate, similar results as for the LDOS in Figure 2 for the largest particles. Interestingly, the radiative decay rate enhancement remains similar for all distances, throughout the different models. This is not unexpected since γ_r takes into account only those modes that radiate to the far-field. It is worth noting in all the spectra the presence of an antiresonance associated to a minimum in the radiative emission rate. This antiresonance behavior corresponds to the situation in which the induced dipole moment has an equal amplitude and opposite phase to the emitting dipole, resulting in an almost-zero radiation to the far-field.

Figure 3:

Normalized total and radiative decay rates for an emitter placed near a metallic sphere (as in Figure 2) with a fixed radius of R=20 nm; the emitter is normally oriented with respect to the metal surface and placed at distances d (shown in nanometers).

This can be easily seen in the dipole-dipole approximation where the radiative decay rate enhancement reduces to [65]:

Assuming a Drude-like dielectric constant for the nanoparticle, i.e. ε=1−ωp2ω2+iγω, we obtain the following antiresonance condition:

As the distance d→0, the antiresonance ω_ar→ω_p. At frequency ω_ar, the radiative decay rate enhancement is minimum and its value is:

From Eq. (6) it can be seen that the minimum value decreases monotonically with d. This is quite evident in Figure 3 for the LRA and TF-HF, but it is not the case for the QHT, where γrγ0|min increases going from d=2 nm to 1 nm. This fact can be related to spill-out effects, where classical formulas break down. In Figure 4 we plot the equilibrium and induced charge densities, n₀ and n₁, respectively, along the z-axis as a function of the distance from the emitter for the QHT. Contrarily to the other models (LRA and TF-HT, for which there is no electron spill-out), in this case the emitter is overlapped with the electron density tail that decays exponentially away from the particle’s ion edge.

Figure 4:

Equilibrium and induced charge densities, n₀ and n₁, respectively, along the z-axis as a function of the distance from the emitter calculated using the quantum hydrodynamic theory (QHT) at ℏω=4 eV.

The dipole emitter is shown in red. The shaded area represents the nanoparticle interior delimitated on the left by the ion edge. The inset shows the absolute value of the induced density in log scale.

It is worth pointing out that although exact electrodynamic solutions exist for both the LRA [66] and the TF-HT [54], [67] cases, their implementation involves infinite sums of higher-order multipoles. While higher-order multipole contribution to radiative rates is negligible, this is not the case for non-radiative rates [68], especially if small distances between the dipole and the metal surface are considered. The numerical evaluation of these large analytical sums may limit the minimum distance that can be practically computed to 1 nm. Our FEM implementation based on the scattering field formalism (see Methods section) allows calculation of the non-radiative decay rates also for sub-nanometer distances.

The quantities analyzed so far are independent (besides a scaling factor) from the dipole moment p of the molecule. The non-perturbative quantities defined in Eqs. (2) and (3), by contrast, strongly depend on these parameters. Figure 5A shows the power and dipole spectra (S and S′, respectively) as a function of the emitter frequency ω₀ for a dipole moment of p=|p|=25.7 D. As the emitter frequency is swept, anti-crossing appears around the main peak in the total emission rate spectrum. Interestingly, this is a dark mode that can be clearly seen in the S maps as a dark horizontal line around ω=ωp/2. The possibility of strong coupling with dark modes was already pointed out by Delga et al. [69] where a collection of quantum emitters interacting with a metal nanoparticle was theoretically analyzed. In contrast, the dipolar and the quadrupolar modes appear as light horizontal lines, although they do not show any strong coupling. Another interesting feature clearly visible in the S maps of Figure 5A is the deep horizontal line associated to the antiresonance in the radiative decay rate spectra. Because this is a far-field effect, there is no strong interaction with the emitter. The strong coupling between the emitter and the nanoparticle is clearly visible for all models in the dipole spectrum S′ maps (second row of Figure 5A).

Figure 5:

Fluorescent emission modification induced by a metallic nanosphere.

(A) Power and dipole strong coupling spectra as a function of the emitter frequency ω₀ for a dipole amplitude of p=25.7 D. (B) Dipole spectra as a function of the oscillator strength f for a dipole placed at distance d=1 nm from a metallic sphere. The emitter frequency is ℏω₀=4.13 eV for the local response approximation (LRA), ℏω₀=4.33 eV for the Thomas-Fermi hydrodynamic theory (TF-HT) and ℏω₀=4.00 eV for the quantum hydrodynamic theory (QHT) calculations, respectively.

The analysis of Rabi splitting is more easily done in terms of the dimensionless oscillator strength f=23mω0p2e2ℏ of the emitting dipole rather than the dipole moment. Figure 5B shows the dipole spectrum S′ as a function of the oscillator strength for ω₀ tuned at the peak of γ_tot. From these maps it is possible to roughly identify the oscillator strength necessary to observe Rabi splitting, namely fthLRA≃0.05, fthTF≃1 and fthQHT≃2. The oscillator strength threshold for the TF-HT and the QHT is more than one order of magnitude bigger than the LRA case. While this is somehow expected in the QHT case because of the extra losses due to viscosity, it is a surprising result for the TF-HT. In fact, although similar systems have already been investigated by some of the authors of this article [70], this result has not been observed before. On the contrary, Tserkezis et al. [56], [71] have demonstrated robustness of strong coupling against nonlocal corrections, when collective emitters are considered. As we will see later, however, this behavior is hindered in the vicinity of noble metal systems, since photon emission spectra are strongly affected by d-band transitions.

2.3 Sphere dimers

Metallic structures that are characterized by few-nanometer or even sub-nanometer inter-particle distances achieve the largest optical field confinement [25], [29] and smallest optical cavity volumes, and are ideal candidates for achieving single-molecule room temperature strong coupling [56], [62], [70], [71], [72]. We consider dimers of Drude-like Na spheres of radius R and gap size g=2d. As in the previous section, it is instructive to first analyze far-field against near-field properties. Figure 6 shows the extinction efficiency and normalized LDOS spectra for dimers characterized by a fixed gap g=2 nm for different sphere radii. The extinction spectra refer to an incident plane wave the electric field of which is polarized along the dimer axis, and the LDOS spectra refer to the states that can be excited by a point-like dipole placed at the center of the gap and oriented along the dimer axis, as shown in the insets of Figure 6A and B, respectively. The QHT results are compared against the LRA and the TF-HT calculations. Contrarily to the single sphere of the previous section, in the dimer system, it is possible to clearly observe higher order resonances as R increases. As expected, the dipole resonance clearly shifts toward lower energies. In particular, in the LDOS spectra, the dipole resonance is always clearly distinguishable (for all models), even for the larger dimers, contrarily to what we observed for the single spheres. As we will show later, this allows a strong coupling of the molecule and the dipolar mode. Differences between the LRA, TF-HT and QHT are analogous to the single sphere case.

Figure 6:

Far-field versus near-field properties for metallic dimers constituted by identical spheres of various radii R, for the local response approximation (LRA), Thomas-Fermi hydrodynamic theory (TF-HT) and quantum hydrodynamic theory (QHT) models, respectively.

(A) The extinction efficiency spectra for a fixed gap size g=2 nm. (B) The normalized local density of states (LDOS) for a dipole parallel to the dimer axis and placed at the center of the gap, i.e. d=1 nm. Material parameters as in Figure 2.

In Figure 7A we plot the spontaneous (total) and radiative emission rate enhancement for different gap sizes and a fixed particle radius for all considered models. Similarly to the single sphere case, at large gaps (and larger emitter distances), all models converge to the same spectrum for both γ_tot and γ_r. At such large gap sizes (g=32 nm), nonlocal and quantum effects are indeed expected to be negligible. As the gap closes, differences between models start to emerge, although, they show qualitatively similar resonances below the plasma frequency. The total decay rate enhancement reaches similar peak values at the dipolar resonance, but it tends to increase at lower energies in the QHT spectra because of the stronger losses due to the viscosity term in Eq. (1). Also in this case, it is possible to identify the antiresonance behavior in the radiative emission rates, although no simple analytical formulas are known for dimers.

Figure 7:

Emission properties of a single molecule placed at the center of a plasmonic dimer.

(A) Normalized total and radiative decay rates for an emitter placed at the center of a metallic sphere dimer as shown in the inset for a fixed radius of R=20 nm; the emitter is parallel to the dimer axis and placed at distances d=g/2 (indicated in nanometers) from the particles. (B) Power and dipole strong coupling spectra as a function of the emitter frequency ω₀ for a dipole amplitude of p=25.7 D.

Let us now study the strong coupling regime. Figure 7B shows the power and dipole spectra as a function of the emitter frequency ω₀ for g=2d=2 nm and a molecule with a dipole moment of p=25.7 D. As anticipated earlier, it is possible in this case to clearly see anti-crossing behavior around the dipolar mode as well other higher order radiative modes, on top of the coupling with the dark mode observed for the spheres. The dipole spectra as a function of the oscillator strength in Figure 8 show that all models predict a similar threshold f_th~0.5 for a gap g=2 nm (or d=1 nm). For a smaller gap, g=1 nm (or d=0.5 nm), the oscillator strength threshold reduces to f_th~0.1 for the LRA and TF-HT, while the QHT spectra become too broad to identify any threshold. Interestingly, when nonlocal broadening is considered, there appears to be an optimal gap distance for being able to observe strong coupling. Figure 9A shows the dipole spectrum for different values of the dimer gap size, i.e. the distance from the emitter. The LRA and the TF-HT predict a similar behavior: as the gap shrinks, the coupling strength increases generating a Rabi splitting that increases indefinitely. In the case of the QHT, the Rabi splitting appears for a sufficiently large coupling but disappears for the smallest gap size, as mentioned above. This peculiar behavior is due to the increase of losses as the emitter overlaps with the electron density tail, as shown in Figure 9B. This is clearly visible in the log-scale plot, where the self-consistent QHT density accounts for a non-zero charge density value in the gap region.

Figure 8:

Dipole spectra as a function of the oscillator strength f for a dipole placed at the center of the sphere dimer.

In (A) g=2d=2 nm and the emitter frequency is ℏω₀=2.42 eV for local response approximation (LRA), ℏℏω₀=2.46 eV for Thomas-Fermi hydrodynamic theory (TF-HT), and ℏω₀=2.38 eV for quantum hydrodynamic theory (QHT) calculations, respectively. In (B) g=2d=1 nm and the emitter frequency is ℏω₀=2.15 eV for LRA, ℏω₀=2.22 eV for TF-HT, and ℏω₀=2.08 eV for QHT calculations, respectively.

Figure 9:

Influence of the electron density tail on the emission properties.

(A) Dipole spectra for different dimer gap sizes and fixed oscillator strength f=2. The oscillator frequency ω₀ is as in Figure 8 for the different models. In (B) the equilibrium electron density is along the dimer axis in the linear (top) and log (bottom) scale. The emitter position is depicted in red.

Although in this work we approximate the emitter as a point-dipole source, in this overlap regime, the effects due to the finite size [21], [73] of the light source might not be negligible, and would require a more sophisticated description of the emitter. One possibility would be to use the QHT to describe both the metal structure and the emitter. In this case, the charge transfer between the emitter and the structure and the nonlocality of the emitter would also be considered.

2.4 Gold structures

So far we have considered Drude-like metal nanoparticles and even if we considered parameters for sodium, our previous results still remain nearly ideal. Although this was useful to understand the underlying behavior of the systems, it is also important to characterize the impact of more practical materials. Any plasmonic experiment is in fact performed using noble metal nanoparticles.

Figure 10A shows decay rate enhancements for a gold sphere dimer characterized by a gap g=2 nm. The spectra are clearly affected by the larger losses due to interband transitions. This also strongly impacts the strong coupling regime as shown in Figure 10B and C. We cannot distinguish an anti-crossing or a Rabi splitting behavior. Instead, as we increase the dipole strength (Figure 10C) we can only observe a shift of the emitter frequency towards lower energies. This is due to the dielectric component (ε_core) of the gold permittivity in the high-end part of the visible spectrum.

Figure 10:

Emission properties of a dipole emitting near the gold systems shown in the insets, for a gap g=2 nm for different models.

(A) Total decay rates and (B) dipole spectra as a function of the emitter frequency ω₀ and the oscillator strength f (C) for a sphere dimer. Total decay rates (D), dipole spectra as a function of the oscillator strength f and local scattering field map (E) for a cone dimer.

In order to avoid this issue, we consider a dimer constituted by two sharp tips or cones, as depicted in the inset of Figure 10D. We limit, however, the sharpness of the structure to a radius of curvature of the tip of r_tip=0.4 nm. The cone dimer allows the dipolar resonance to shift towards lower energies, increase its strength and at the same time reduce the mode volume. These properties allow observation of a Rabi splitting for oscillator strength of the order of f_th~0.1 for the LRA and TF-HT and f_th~0.2 for the QHT, as shown in Figure 10E. Figure 10F shows the norm of the scattered field |E_s|=|E–E_dipole|, E_dipole being the electric field associated to the dipole. The impact of the electron density spill-out is clearly visible for the QHT maps.

Although the impact of d-band electrons is detrimental for observing strong coupling, atomistic features [42], [74] in real systems might provide some extra field confinement, allowing strong coupling observation in real systems [72]. Another important role could be played by the spatial extension of real molecules [21], which goes beyond the point dipole approximation used in this work.