Strong Coupling between Localized Surface Plasmons and Molecules by Coupled Cluster Theory

Plasmonic nanocavities enable the confinement of molecules and electromagnetic fields within nanometric volumes. As a consequence, the molecules experience a remarkably strong interaction with the electromagnetic field to such an extent that the quantum states of the system become hybrids between light and matter: polaritons. Here, we present a nonperturbative method to simulate the emerging properties of such polaritons: it combines a high-level quantum chemical description of the molecule with a quantized description of the localized surface plasmons in the nanocavity. We apply the method to molecules of realistic complexity in a typical plasmonic nanocavity, featuring also a subnanometric asperity (picocavity). Our results disclose the effects of the mutual polarization and correlation of plasmons and molecular excitations, disregarded so far. They also quantify to what extent the molecular charge density can be manipulated by nanocavities and stand as benchmarks to guide the development of methods for molecular polaritonics.

S trong coupling between molecules and quantum plasmons 1 in nanocavities 2 leads to the formation of hybrid plasmon−molecule states: polaritons, or more specifically, plexcitons. 3 These new states manifest distinct features compared to the original states, 4−7 potentially resulting in modified chemical/photochemical reactivity 8−10 and relaxation dynamics, 11,12 along with other coherent processes. 13,14 Modeling accurately the molecules, nanostructures and their coupling is of the utmost importance to support the experimental advances. Coupling descriptions proper for resonant optical cavities, such as using the molecular dipole only, 15,16 should be amended 17−19 when the coupling affects the molecule on a submolecular level. 18,20,21 There are currently no approaches that treat the coupled plasmon−molecule system nonperturbatively, that is, that include the relaxation of the molecular electronic density upon polariton formation together with plasmon−molecule correlation. In this work we extend the QED-CC method, 22,23 already applied to resonant optical cavities, to realistic nanoplasmonic cavities. We quantize the plasmonic excitations starting from the classical dielectric description of nanostructures with arbitrary shapes. Each quantized plasmonic mode is associated with a surface charge density (discretized into point charges), analogous of a transition charge density in a molecular system. A similar quantization approach for plasmons in Drude metals, derived from macroscopic QED, 19,24,25 was adopted in ref 17. The quantization framework we are presenting accounts for the geometrical features of the nanoparticles setups without disregarding the molecular complexity. Indeed, coupled cluster (singles and doubles excitations) is recognized as a highly accurate method in quantum chemistry: 26 we extend it to include self-consistent and correlated molecule−plasmon hybridization. With this method, we compute the interaction between a nanocavity realized after ref 18 and two realistic molecules: porphyrin and para-nitroaniline (PNA). For porphyrin, we highlight the role of the geometrical features of the system and the role of electron−plasmon self-consistency and correlation in polaritons, beyond the standard Jaynes-Cumming picture. 27 We use PNA to quantify the same effects, already predicted on the ground and excited states electron densities 22 for a resonant cavity, in the case of a nanocavity.

■ RESULTS AND DISCUSSION
Theory. The first subsystem we focus on is the nanostructure. We start from a classical dielectric perspective, namely computing the plasmon modes from a continuous medium model. As such, it is not an atomistic description and does not include dielectric nonlocal effects, electron spillout, or chemical bonds. However, the description we adopt is surprisingly accurate: 28 it reproduces the inhomogeneous electric field due to atomistic irregularities in the nanostructures 29 even at subnanometric level. 30 Our goal is to describe the nanostructure at a quantum level as a material system, which differs from other approaches that quantize the electromagnetic field in the presence of dielectric media. 24,31 By building on polarizable continuum model applied to nanoparticles (PCM-NP) approach, 32 we shall achieve quantization (Q-PCM-NP) in the quasi-electrostatic framework, where the retardation effects are disregarded. 33 Similar quantization frameworks for the plasmonic modes were used in different contexts not involving a molecule 34 or using a model description of the molecule. 35 To classically compute the nanoparticle dielectric response properties, we solve Maxwell's equations for a continuous, frequency-dependent dielectric (nanostructure) under an external electromagnetic perturbation. In the quasi-static limit, the nanoparticle experiences an electrostatic potential at a given frequency, V(ω), which induces polarization charges q(ω) on the nanoparticle surface. The equation defining such response charges is where Q IEF (ω) is the frequency-dependent response to the external perturbation V(ω). We refer to IEF as to the specific integral equation formalism 36 formulation of the PCM problem. In the present work, V(ω) is the potential produced by the molecular transition density. The frequency-dependent linear charge density response function of the nanostructure is taken in a diagonal form 37 where K is the diagonal linear-charge-response-matrix derived from eigenvalues Λ p of the appropriate IEF matrix, 37 ε(ω) is the frequency-dependent dielectric function, and S is the matrix storing the electrostatic potential between discrete points of the dielectric. This is the starting point of the Q-PCM-NP quantization procedure. Quantization of the Plasmonic Modes: Q-PCM-NP. We assume that the metal nanoparticle can be characterized by a Drude-Lorentz dielectric function where Ω P 2 is the squared plasma frequency of the bulk metal, ω 0 is the natural frequency of the bound oscillators (Lorentz model), and γ is the damping rate. By defining Each matrix element of Q IEF (ω) is evaluated on representative points of the tesserae, k and j. In eq 5, we recognize the spectral form of a linear response function, 38 Q quant (ω). Our goal now is to identify the quantities characterizing the excited states of the nanostructure (that is, plasmons) relevant to devise the coupling with the molecule. For the quantum description of the nanostructure, the response function Q quant (ω) is formally written in terms of the surface charge operator q. 39 We conveniently write Q quant (ω) matrix elements in its spectral representation 40 as By comparing eq 5 and eq 6, we verify that ω p represents the excitation frequencies of the plasmonic system. γ p ′ = γ/2 for each plasmon state p is a decay rate, while we can identify ⟨0|q̂k|p⟩ with The label q p,k represents the transition charge sitting on the kth tessera associated with the mode p and they are the analogous of transition densities in molecules. Expressed differently, the ensemble of the q p,k set of charges represents the normal mode of the plasmonic system at frequency ω p (see Figure 1a−c as an example).
We are now in the position to write the quantum plasmonic Hamiltonian as where ω p is the p-mode frequency and b̂p † , b̂p are the corresponding bosonic creation and annihilation operators. 17,25 Assuming the localized surface plasmons to be bosons is not a result of the present derivation (at the linear response level, a Fermionic or bosonic excitation would yield the same Q quant (ω)). It is rather based on the general behavior of plasmons in extended systems. 1 In the Supporting Information, we report numerical checks of the present approach against analytical results for a point-dipole molecule interacting with a Drude metal nanosphere. 41 We can now write the general quantized plasmon−molecule Hamiltonian aŝ=̂+̂+Ĥ where Ĥm ol is the standard molecular Hamiltonian (which we reduce to the electronic Hamiltonian). Making use of the response charges for the quantized plasmon, the interaction term Ĥi nt reads 39 where V̂j is the molecular electrostatic potential operator, evaluated at each tessera representative point j. Our formulation of Ĥi nt presents two major advances with respect to simplified models (see Supporting Information for the formal comparison): we obtain an intuitive representation of the inhomogeneous electromagnetic field in terms of polarization charges and the interaction is general on the molecular side, meaning that it can be interfaced to advanced electronic structure models. Through this interface, our method can describe phenomena occurring in both the weak and strong coupling regime. 42,43 Radiative corrections to the quasi-static description, related to radiation reaction, 44 can also be included in this model, at the quantum level. The plasmon radiative decay rate can be calculated from the transition dipole moment associated with the charges q p,j and it can then be added to the nonradiative lifetime γ p . For the systems treated in the numerical section, this correction is negligible. An estimate of the molecule−plasmon coupling g n between the plasmonic mode p and the molecular electronic transition from the ground state S 0 to the excited state S n is given by Here 0 and 1 are the occupation numbers of the plasmonic mode p. The coupling g n is the central quantity in simplified approaches to strong coupling, such as the Jaynes-Cummings model (JC).
Extending QED Coupled Cluster (QED-CC) to Plasmon−Molecule Systems. Determining the eigenfunctions of the plasmon−molecule Hamiltonian in eq 9 requires an accurate description of the molecule. In this section, we outline the extension of the QED-CC treatment (already applied to resonant cavities 22,23 ) to the case of plasmons. The detailed algorithm is available in the QED-CC dedicated section of the Supporting Information. The simplest description of the molecule without plasmon interactions is a single Slater determinant, |HF⟩ (Hartree−Fock). To include electron correlation we rely on coupled cluster theory. 26 More explicitly, we apply the exponential of the cluster operator to the |HF⟩ determinant The cluster operator T̂is expressed aŝ where T 1 are linear combinations of single excitations, T 2 are linear combinations of double excitations and so on. Inspired by the formal similarities of the interaction in molecule− plasmon and molecule−photon systems, we extend the newly developed quantum electrodynamics coupled cluster theory (QED-CC) 22 from photons to plasmons. Here the plasmon− molecule interaction is described nonperturbatively, allowing also the ground state to couple to the electromagnetic field. 23 In QED-CC, the cluster operator T̂also includes the bosonic creation and annihilation operators bp † , bp of the quantized plasmon mode p̂=̂+̂+̂+̂+ Here we identify three contributions: electronic (T1, T2), plasmonic (Γ̂1) and mixed plasmon−molecule excitations (S1 1 , S2 1 ). By the effect of such excitations, QED-CCSD-1 (SD denotes single and double excitation of the electronic subspace, whereas 1 is the excitation order of the plasmonic modes) includes the correlation induced on the electronic states by the plasmonic transition. We note that, in the limit where the interaction ĝ→ 0, CCSD and QED-CCSD-1 are equivalent.
Calculations with Molecules. Free-Base Porphyrin. In this section, we exploit Q-PCM-NP coupled to QED-CCSD-1 and EOM-CCSD to simulate plasmon−molecule interactions. Here we consider a free-base porphyrin coupled to the plasmonic nanotip, as sketched in Figure 2a. A similar setup was adopted in refs 17 and 18. According to eq 11, we define the couplings between the ground state |S 0 , 1⟩ of the free-base porphyrin and the |S 1 , 0⟩, |S 2 , 0⟩ states as g 1 and g 2 respectively. In Figure 2b, we show the associated interaction maps, where the couplings g 1 and g 2 are displayed as a function of the molecule−nanotip displacement. The two quasi-degenerate dipolar transitions of free-base porphyrin at 4.05 eV (S 1 ) and 4.07 eV (S 2 ) have different symmetries. Depending on the tip displacement, the interaction element g varies between 0 meV at the center to 36 meV at 6 Å displacement. By displacing the tip accordingly, a polariton is created with either of the two states. When combining the two quasi-degenerate transitions, = + g g g eff 1 2 2 2 , a doughnut-like shape of the transition density is obtained, in agreement with previous results on zinc phtalocyanine. 17 We stress that this doughnut-like shape is not observable without including the geometry of the nanostructure. In Figure 2c, we show the absorption spectrum of the relaxed tip−molecule complex (QED-CCSD-1) and compare it with the polaritons obtained by the Jaynes-Cummings model (see the Supporting Information for details). The point at which we perform the calculations is highlighted as the black Nano Letters pubs.acs.org/NanoLett Letter square in panel g eff of Figure 2b. In QED-CCSD-1, the light− matter system is allowed to correlate both the ground and the excited states, guaranteeing accurate results with respect to the truncated basis. Even more, the interactions between all singly and doubly excited determinants with and without plasmons are included. Conversely, the JC case is treated by correlating the electronic states only, hence neglecting the mixed molecule−plasmon correlation contributions. As a further difference between JC and QED-CCSD-1, the interaction in the former is a simple norm of the g 1 and g 2 interactions, and it does not allow for the S 1 and S 2 states to interact through the plasmon. For QED-CCSD-1, the presence of the polaritonic correlation and several electronic states reveals that the plasmon contribution is distributed over other electronic states displayed as the small peak at 3.88 eV in Figure 2c. At the same time, the Rabi splitting between the two main polaritonic peaks is reduced (46 meV) with respect to JC (70 meV). While the small peak is due to the inclusion of different states, the difference in the Rabi splitting is due to the inclusion of the mixed plasmon−molecule correlation. To corroborate this view, we compare in the Supporting Information an extended JC model using three states with the QED-CCSD-1.
Although the description improves with the extended JC model, it still misses the intensities and the positions of the polaritonic peaks, confirming the role of the mixed plasmon− molecule correlation in the description. Para-nitroaniline. Para-nitroaniline (PNA) presents an intense absorption peak in the UV range. We hence expect large Rabi splittings when the PNA is appropriately aligned with the electromagnetic field of the nanotip. The transition density associated with the bright transition of PNA at ∼4.5 eV is parallel to the long axis of the molecule, making it an excellent candidate for single-molecule strong coupling.
In Figure 3, we investigate the coupling of the PNA molecule with the nanotip. We identify the maximum (parallel) and the minimum (perpendicular) coupling conditions, as sketched in Figure 3a, and plot the Rabi splitting of the plasmon−molecule coupled system in Figure  3b. When the molecular transition dipole is perpendicular to the plasmon mode, there is little-to-no effect on the spectrum. Instead, when the molecule is gradually rotated to be parallel to the mode, the Rabi splitting increases until the maximum of 179 meV (parallel). By gradually rotating the molecule in the opposite direction (antiparallel) we observe a slightly smaller splitting (175 meV). The 4 meV difference between the parallel and antiparallel orientations is a signature of the realistic molecular description. Indeed, a dipolar interaction would result in a symmetric spectrum with a Rabi splitting independent of the parallel/antiparallel orientation. In addition, we observe a red-shift of the polaritonic peaks when including the mixed plasmon−molecule correlation energy via QED-CCSD-1.
The Jaynes-Cummings model fails at describing groundstate properties, as no change in the molecular ground state is predicted by the model. By instead using QED-CCSD-1, we show in Figure 3c how the ground-state electronic density changes due to the molecule−plasmon interaction. Since the coupling in the present case is smaller than the one considered in ref 22 the charge localization effects as discussed in the same reference are reduced. We quantify this effect in Figure 4a. We note the opposite trend to what was described in ref 22, namely that here the plasmon induces a small charge transfer from the donor to the acceptor. This inversion is a signature that the molecules interact differently with plasmonic nanotips compared to optical cavities, whereas the Jaynes-Cummings model treats them equivalently. The charge transfer properties of the polariton states are also substantially changed, reminiscent of what is shown in ref 22. However, in the present case the sum of the polaritons' charge transfer adds up to the charge transfer of the singlet state without the plasmons, as displayed in Figure 4b. This again highlights that interactions with nanoparticle plasmons and optical cavities are not equivalent. The plasmon-induced changes in the electronic density also imply a modification of other groundstate properties as well as the potential energy surfaces. One example is the minor induced change of the dipole moment (− 0.03 D). One important limitation in this ground-state study is the inclusion of only a single plasmon mode: higher energy modes are expected to contribute to the ground-state coupling and enhance the minor effects seen here, as the ground-state correlation is not a resonant property. 23

■ CONCLUSIONS
In this work, we have introduced a new methodology, based on PCM-NP 45,46 and QED-CC, 22 to treat the quantum coupling between nanoparticles of arbitrary shapes and molecules. We compared our description of polaritons in free-base porphyrin with the simplified Jaynes-Cummings model, highlighting the role of the mixed molecule−plasmon correlation. Finally, we have analyzed the effects on the ground state of PNA in a realistic plasmonic nanocavity setup, showing that the modification of ground states density previously investigated by QED-CC is still present. Our representation of quantum plasmonic modes as apparent surface charges, in synergy with coupled cluster, provides a simple-yet-accurate interface to tackle mixed nanoparticle−molecules systems. In conclusion, we believe our method provides solid ground to support and guide experiments investigating and manipulating the chemical properties at the submolecular level. Accurate calculations of plexciton-affected reaction-barriers are already possible, as well as exploring effects related to chiral nanotips. The method is also prone to various extensions on the nanostructure side: the inclusion of multiple EM modes or the extension to classical atomistic models. 47