Strong plasmon-molecule coupling at the nanoscale revealed by first-principles modeling

TL;DR

This paper demonstrates that first-principles TDDFT calculations can accurately model strong plasmon-molecule interactions at the nanoscale, predicting significant vacuum Rabi splitting and enabling detailed studies of polaritonic systems.

Contribution

It introduces a parameter-free ab initio approach using TDDFT to analyze nanoscale plasmon-molecule hybrids, bridging quantum optics and chemistry.

Findings

Vacuum Rabi splitting predicted in small plasmon-molecule systems

Single-molecule coupling strength exceeds 200 meV

Cavity quantum electrodynamics applies at nanometer scales

Abstract

Strong light-matter interactions in both the single-emitter and collective strong coupling regimes attract significant attention due to emerging quantum and nonlinear optics applications, as well as opportunities for modifying material-related properties. Further exploration of these phenomena requires an appropriate theoretical methodology, which is demanding since polaritons are at the intersection between quantum optics, solid state physics and quantum chemistry. Fortunately, however, nanoscale polaritons can be realized in small plasmon-molecule systems, which in principle allows treating them using ab initio methods, although this has not been demonstrated to date. Here, we show that time-dependent density-functional theory (TDDFT) calculations can access the physics of nanoscale plasmon-molecule hybrids and predict vacuum Rabi splitting in a system comprising a few-hundred-atom aluminum nanoparticle interacting with one or several benzene molecules. We show that the cavity quantum electrodynamics approach holds down to resonators on the order of a few cubic nanometers, yielding a single-molecule coupling strength exceeding 200 meV due to a massive vacuum field value of 4.5 V/nm. In a broader perspective, our approach enables parameter-free in-depth studies of polaritonic systems, including ground state, chemical and thermodynamic modifications of the molecules in the strong-coupling regime, which may find important use in emerging applications such as cavity enhanced catalysis.

Equations

$$g(N)=\sqrt{N}{\mu }_1{E}_{\mathrm{vac}},$$

Full text

Light-matter interaction between an optical mode and a quantum emitter is accurately described in terms of the cavity quantum electrodynamics (cQED) formalism Scully and Zubairy (1997). The interaction can be weak or strong, depending on circumstances, leading to drastically different behavior. The regime of strong light-matter coupling, unlike its weak counterpart, leads to formation of hybrid cavity-emitter eigenmodes manifested in coherent energy exchange between the system subsystems occurring on timescales that are (much) faster than the corresponding damping rates. Thus, the emitter and the cavity form a unified light-matter hybrid polariton, whose properties, including spontaneous emission and chemical potential, can be tuned Törmä and Barnes (2015); Baranov et al. (2018); Kockum et al. (2019). Because of their compositional nature, polaritons are useful for photon-photon interactions, which leads to remarkable nonlinear and quantum optical phenomena, including photon blockade and Bose-Einstein condensation of exciton polaritons Yoshie et al. (2004); Birnbaum et al. (2005); Kasprzak et al. (2006). On the other hand, strong light-matter coupling may lead to changes in emitter properties, including photochemical rates Hutchison et al. (2012); Herrera and Spano (2016); Galego et al. (2016); Martínez-Martínez et al. (2018); Munkhbat et al. (2018) and exciton transport Feist and Garcia-Vidal (2015); Zhong et al. (2016).

Traditionally cavity-exciton polaritons were realized in atomic or solid state systems utilizing high quality factor optical cavities at low temperatures. These efforts resulted in remarkable effects such as polariton Bose-Einstein condensates and superfluidity Khitrova et al. (2006). Polaritonic behavior associated with “dressing” of the emitter by a cavity field is usually captured by traditional quantum optical approaches such as Jaynes-Cummings or Dicke models Jaynes and Cummings (1963); Garraway (2011); Törmä and Barnes (2015); Baranov et al. (2018); Kockum et al. (2019). However, these quantum optical formalisms treat matter in an extremely simplified manner, that is, as a two-level system, leading to oversimplifications and inconsistencies in the description of material subpart.

More advanced theoretical techniques developed recently allow for more sophisticated effects including, e.g., multiple electronic resonances, accounting for atomic vibrations, and light-matter interactions beyond the point dipole approximation. Significant progress along these lines has been achieved by several groups using various quantum optical and quantum chemistry methods Herrera and Spano (2016); Galego et al. (2016); Martínez-Martínez et al. (2018); Galego et al. (2015); Ruggenthaler et al. (2014); Ling Luk et al. (2017); Neuman et al. (2018). However, typically either molecules or electromagnetic fields in these approaches are treated in a simplified manner. Here, we model the entire plasmon-exciton system by time-dependent density-functional theory (TDDFT) Runge and Gross (1984).

Plasmon-molecule interactions can be modeled computationally using density-functional theory (DFT) Hohenberg and Kohn (1964); Kohn and Sham (1965) and/or TDDFT approaches, thanks to the relatively small number of atoms involved in these interactions. TDDFT allows the whole system to be treated on the same footing, enabling one to track effects related to modification of the matter subpart, which are inaccessible by purely quantum optical or classical electromagnetism methods. In fact, analytical and numerical quantum chemistry approaches, including DFT and TDDFT, have been successfully applied to study plasmon-molecule interactions, charge transfer, chemical enhancement and electromagnetic effects in surface-enhanced Raman scattering (SERS) experiments Arenas et al. (1996, 2002); Bjerneld et al. (2004); Zhao et al. (2006); Jensen et al. (2008); Morton and Jensen (2011); Corni and Tomasi (2001); Lombardi and Birke (2008); Shegai et al. (2009). However, vacuum Rabi splitting and strong plasmon-molecule coupling have not been in the focus.

Related to the above, quantum electrodynamics density-functional theory (QED-DFT) has been developed recently Ruggenthaler et al. (2014), enabling the modeling of strong coupling between electromagnetic cavity modes and electronic excitations Flick et al. (2017, 2018). While the commonly-used TDDFT approaches with classical description of electric fields may not be able to describe the full range of quantum optical phenomena that QED-DFT aims at, we demonstrate in this paper that the cavity mode created by a localized surface plasmon resonance as well as its strong coupling with excitons can be described already within the standard TDDFT.

Additionally, many recent experimental observations cannot be easily explained with current theories Hutchison et al. (2013); Orgiu et al. (2015); Thomas et al. (2016). Here, we take a step towards helping to elucidate such observations by demonstrating the possibilities of TDDFT for studying polariton physics. We argue that theoretical predictions obtained on this small scale may be useful for understanding and modeling of more extended systems.

We model light-matter interactions by employing the real-time-propagation (RT) TDDFT approach Yabana and Bertsch (1996) based on the localized basis sets Kuisma et al. (2015); Larsen et al. (2009) as implemented in the open-source GPAW package Mortensen et al. (2005); Enkovaara et al. (2010). This RT-TDDFT code is combined with extensive analysis tools Rossi et al. (2017), that are utilized for analyzing the electron-hole transition contributions to resonances and visualizing them as transition contribution maps (TCM) Malola et al. (2013); Rossi et al. (2017). See Methods for detailed description.

We focus on proof-of-principle plasmonic systems based on idealized aluminum (Al) nanoparticles. The free-electron-like electronic structure of Al greatly simplifies the analysis in contrast to noble metals with $d$-electron-screened plasmons. In addition, Al has attracted recent interest as an alternative plasmonic material Knight et al. (2014).

Our model systems consist of Al nanoparticles (regular truncated octahedra with 201, 586, and 1289 atoms) and benzene molecules (see Fig. S4). The \ceAl201 nanoparticle exhibits a plasmon resonance at $\hbar\omega_{\mathrm{Al}}=7.7\,\mathrm{eV}$ (Fig. 1a), whose collective nature is recognizable in the TCM (Fig. 1b) Rossi et al. (2017). The benzene molecule has a doubly-degenerate resonance at $\hbar\omega_{\mathrm{B}}=7.1\,\mathrm{eV}$ with a transition dipole moment $\mu_{1}=4.45\,\mathrm{D}$ (Figs. 1a, S1). When the nanoparticle and benzene molecules are placed in proximity (bottom part of Fig. 1a) their resonances couple and form two polaritons. The lower (at 6.92 eV) and upper (at 7.68 eV) polaritons are, respectively, mixed symmetric and antisymmetric plasmon-molecule states as apparent in the TCMs of the coupled system (Fig. 1de). A comparison with the TCMs of the uncoupled constituents (Fig. 1bc) reveals the individual contributions to the polaritons. Transitions from just below the Fermi energy to just above correspond to the plasmon, while the contribution from about $-2$ to 3 eV originates from the molecular exciton. Crucially, the two polaritonic states differ with respect to the sign of the molecular contributions. In case of the lower polariton (LP) plasmonic and molecular transitions are in-phase, while for the upper polariton (UP) they are out-of-phase. This respectively symmetric and antisymmetric combination is clearly visible in the induced densities. At the LP the dipoles of both particle and molecules are parallel, while they are antiparallel at the UP. Such an arrangement is archetypal for a strongly coupled system. This observation is valid also for all other calculated spectra, indicating the presence of strong or near-strong coupling already in the case of a single molecule placed 3 Å away from the \ceAl201 particle. We note here, that the density of states of the molecular elements in the coupled system does not significantly differ from their uncoupled states up to and including the LUMO (Fig. S2).

To determine the fidelity of TDDFT for modeling strong coupling, we calculated the photoabsorption spectra of \ceAl201, \ceAl586, and \ceAl1289 coupled to $N$ benzene molecules, whose dipole moment is aligned collinearly with that of the plasmon vacuum field. The molecules are positioned at the corners of two opposing {100} facets of the particles (Fig. S4), for simplicity neglecting relaxation. The coupled systems exhibit the emergence of the LP and UP modes, the splitting of which increases with the number of molecules, see Fig. 2 (TCMs show qualitatively identical symmetric/antisymmetric mixtures of plasmonic/molecular states, see Fig. S3). The Rabi splitting $2\Omega$ of the absorption spectra for \ceAl201 and one benzene molecule equals 730 meV, a value comparable to the plasmon width. A non-negligible part of the UP/LP separation originates, however, from the plasmon-exciton detuning and the coupling strength $g$ is 200 meV for one benzene, a value smaller than the geometrical mean of the plasmon and exciton widths (300 meV). The strong coupling condition is, however, fulfilled for $N\geq 3$ for \ceAl201 and for still larger $N$ for larger particles.

One of the signatures of strong coupling is the square-root dependence of $g$ on the number $N$ of identical excitons with transition dipole moment $\mu_{1}$ interacting with an optical cavity in the standard QED expression,

[TABLE]

where $E_{\mathrm{vac}}=\sqrt{\hbar\omega_{\mathrm{Al}}/2\varepsilon\varepsilon_{0}V}$ is the vacuum field and $V$ is the volume of the electromagnetic mode Khitrova et al. (2006). Using the \ceAl201 particle as an example, the expected scaling is shown in Fig. 3 by the black dotted line using $\mu_{1}$ of a benzene molecule and $E_{\mathrm{vac}}^{(201)}=4.45\,\nicefrac{\mathrm{V}}{\mathrm{nm}}$ with $\hbar\omega_{\mathrm{Al}}=7.7\,\mathrm{eV}$ and $V=3.3\,\mathrm{nm}^{3}$ corresponding to the volume of the \ceAl201 particle, which is an adequate estimate for a nanometer-sized particle Koenderink (2010). To estimate the coupling strengths $g$ from our first-principles data, we fit the absorption spectra with a coupled harmonic oscillator model Wu et al. (2010) (see Sec. S4 in the SI and Figs. S5-S8), yielding the $g$-values shown by squares in Fig. 3. We notice that the coupling strength values obtained using the standard QED expression in Eq. (1), which is insensitive to spatial variations in the electric field, are larger than the ones derived from coupled oscillator fitting. A similar observation holds for the larger particles, except that due to a larger volume vacuum field and coupling strengths are reduced.

Equation (1) is accurate under the condition that all excitonic modes are coherently coupled with the same, maximal rate to the cavity. This is true for structures such as Fabry-Pérot cavities, photonic crystal slabs, or micropillars, in which the anti-node of the mode is accessible to excitonic modes of molecules or semiconductors by precise placement via trapping or doping Khitrova et al. (2006); Baranov et al. (2018). For isolated plasmonic particles the maximum of its mode is, however, typically inside the particle Koenderink (2010), such as in the present case. Consequently, the molecular coherent dipole moment $\mu_{\mathrm{coh}}$ interacting with the cavity is smaller than its maximum value of 4.45 D. It is in the range 2 to 4 D and increases with the number of Al atoms in the considered range of particle volumes. The reduction of efficiency is typically expressed in terms of an efficiency factor $\eta$ determined by the ratio of the mode energy density at the position of the exciton to its maximum value Yang et al. (2016).

Before extracting efficiency factors $\eta$ from the TDDFT data, we estimate their values on the basis of classical electromagnetic calculations, which are acceptable at a semi-quantitative level for picocavities Benz et al. (2016), with a local Drude-permittivity tailored to match the TDDFT absorption spectrum for \ceAl201. We find that in the spatial region occupied by the benzene molecule $\eta$ ranges from 0.2 to 0.4 (Fig. 4a), due to the rapid decrease of the plasmon-induced electric field. However, such a low value of $\eta$ predicted by calculating the mode profile of a bare cavity does not hold in practice. It is known that the presence of a material with a refractive index greater than the background causes additional field localization, increasing the coupling strength Yang et al. (2016). While a single benzene molecule cannot be meaningfully described in terms of a refractive index, it has the identical effect. This is confirmed by the electric field enhancements from TDDFT (Fig. 4e-g). In comparison to the induced field of the bare plasmon, in the LP/UP of the coupled system benzene focuses the electric field, modifying the cavity and, consequently, its vacuum field Yang et al. (2016).

We now turn to the extraction of the coupling efficiencies $\eta$ from the TDDFT-derived coupling strengths. The efficiencies modify the slope predicted by Eq. (1) and the resulting dependencies are shown as dashed lines in Fig. 3. $\eta_{201}=0.47$ gives very good agreement for $N\leq 3$ (note the dashed black line marking the theoretical dependence $g=\eta_{201}\sqrt{N}\mu_{1}E_{\mathrm{vac}}$ in Fig. 3). For $N>3$ the calculated coupling strengths fall, however, below the slope set by $\eta_{201}$. Such a deviation for large $N$ is not observed for the two larger particles, which, respectively, have efficiency factors of $\eta_{586}=0.61$ and $\eta_{1289}=0.78$. These larger efficiencies imply that the plasmonic modes of the larger particles are used more efficiently than in the case of \ceAl201. Minor deviations for the two larger nanoparticles are seen for $N<3$. However, they are probably caused by inaccurate fitting of the absorption spectra, as the LP/UP are not very clearly defined.

The efficiency factors obtained from classical electromagnetic and TDDFT simulations are consistent. We note that in both scenarios the addition of molecules leads to a systematic red shift of the plasmon (Fig. S9). This shift is the largest for \ceAl201 and scales inversely with particle volume in agreement with predictions Antosiewicz et al. (2012), further supporting the validity of our calculations. This effect goes hand-in-hand with strong coupling and complicates the analysis. The calculations yield a $\sqrt{N}$ dependence of the coupling strength for the larger particles (\ceAl586, \ceAl1289). For the smallest particle (\ceAl201) the efficiency, however, varies with the number of benzene molecules for $N\gtrsim 3$. In this case, the volume occupied by the molecules is non-negligible relative to the particle volume and every additional molecule changes the properties of the nanoscale cavity and subsequently modifies the coupling.

We will now interpret these results from the perspective of classical electromagnetism. The molecules can be represented as a polarizable background with excitons modelled by Lorentzians Antosiewicz et al. (2014). This background can be treated as part of the cavity and consequently modifies its local optical density of states (LDoS) Yang et al. (2016); Zhang et al. (2015). This modification manifests itself as focusing of the electric field around the molecules, clearly visible for \ceAl201 coupled with one benzene molecule (Fig. 4cd), in comparison to the resonance of \ceAl201 alone (Fig. 4b). These LDoS modifications cause an increase of the mode volume, decreasing (or at least not enhancing) the coupling strength per molecule. Such changes of the cavity induced by adding molecules are not visible for larger particles, as for \ceAl586 and \ceAl1289 the coupling strengths follow the $\sqrt{N}$ dependence.

A final point we discuss is the effect of a diminishing overlap between the vacuum field of the cavity and the molecular exciton. It is cleanly demonstrated via the distance dependence of the response, which we performed for the \ceAl201 particle coupled to two benzene molecules (Fig. 5a). The Rabi splitting of LP/UP decreases with increasing separation while the TCMs for LP/UP show qualitatively identical symmetric/antisymmetric mixtures of plasmonic/molecular states. Furthermore, the reduction of $\eta$ Yang et al. (2016) leads to a decrease of the coupling strength (Fig. 5b).

In conclusion, we demonstrated the suitability of first-principles TDDFT for studying strong coupling between plasmons and molecular excitons. This approach allows us to capture relevant interactions at the atomic level. This is important for studying coupling of molecules to picocavities slightly larger than the molecules themselves Benz et al. (2016); Chikkaraddy et al. (2016)) as well as ultrastrong coupling that leads to modification of the molecular ground state Galego et al. (2015); George et al. (2015). Furthermore, we have shown the degree to which the simple cQED description holds for small systems. For small single-particle cavities, such as the \ceAl201 particle considered here, the presence of molecules modifies the cavity and, consequently, its mode volume. The mode volume, in turn, affects the coupling strength, whose value does not increase as quickly as predicted by QED. This slower increase of $g$ is noticeable for the \ceAl201 particle, but for larger sizes any deviations from Eq. (1) become negligible. Despite the deviations between a simple cQED formula and TDDFT calculations, the order of magnitude agreement between them is quite remarkable. In particular, it is rather surprising that such a simple formula at all holds at such small scales. Based on the obtained coupling strengths, calculated mode volumes and field enhancements, reaching single/few-molecule ultrastrong coupling for single-particle cavities may prove challenging.

We acknowledge financial support from the Swedish Research Council (VR), the Knut and Alice Wallenberg Foundation, and the Polish National Science Center via the project 2017/25/B/ST3/00744. Computational resources were kindly provided by the ICM-UW (Grant #G55-6) and by the Swedish National Infrastructure for Computing at PDC, Stockholm. TJA also acknowledges support from the HPC-Europa3 program via project HPC17XPO18.

Methods

Computational details. The DFT and TDDFT calculations were carried out using the Perdew-Burke-Ernzerhof (PBE) Perdew et al. (1996); *1997_PRL_78_1396_perdew exchange-correlation functional in the adiabatic limit. The spectra are calculated using the $\delta$-kick technique Yabana and Bertsch (1996) in the linear-response regime and employing the dipole approximation for light-matter interaction. The default projector augmented-wave (PAW) Blöchl (1994) data sets and double-$\zeta$ polarized (dzp) basis sets provided in GPAW were used. The dzp basis set of Al includes diffuse 3p functions, which are important for describing plasmon resonances Rossi et al. (2015). To minimize spurious effects due to the basis-set superposition error, the so-called ghost-atom approach was used separately for each nanoparticle size to keep the total system basis set as intact as possible despite the changing number of surrounding molecules. In general, while the used basis sets might not be adequate for yielding numerical values at the complete-basis-set limit, they are expected to be sufficient for the purposes of the present work. A grid spacing parameter of 0.3 Å was chosen to represent densities and potentials, and the molecules/particles were surrounded by a vacuum region of at least 6 Å. The Hartree potential was evaluated on a larger grid with at least 120 Å vacuum around the system and a coarser grid spacing of 1.2 Å, and subsequently refined to the original grid. For the time propagation, we used a time step of $\Delta t=15\,\mathrm{as}$ and total propagation time of at least $T=30\,\mathrm{fs}$, which is sufficient for the used Lorentzian spectral broadening with $\eta=0.1\,\mathrm{eV}$ corresponding to a full width at half-maximum of 0.2 eV.

Transition contribution map (TCM). A TCM is used for visualizing the Kohn-Sham (KS) electron-hole transition contributions to photoabsorption. Briefly, the photoabsorption cross-section $S(\omega)$ is expressed in the basis of occupied ($i$) and unoccupied ($a$) KS states as $S(\omega)=\sum_{ia}S_{ia}(\omega)$ Rossi et al. (2017). In TCM, the elements of the matrix $S_{ia}(\omega)$ at a chosen resonance energy are plotted on a Gaussian-broadened two-dimensional plane spanned by the energy axes for occupied ($\varepsilon_{o}$) and unoccupied ($\varepsilon_{u}$) KS states. See Ref. Rossi et al., 2017 for a detailed description of TCM construction. The maps are augmented with the corresponding densities of states, and a diagonal line is drawn to indicate the KS eigenvalue difference corresponding to $\omega$.