Dipolar coupling of nanoparticle-molecule assemblies: An efficient approach for studying strong coupling

The coupling of light and matter in the strong coupling (SC) regime leads to the emergence of excited states of mixed nature,¹ which are characterized by a coherent energy exchange between the subsystems at a rate that is much faster than the respective damping rates. Emitter and the cavity thus form a light–matter hybrid (polariton) with modified and tunable properties,^2–4 including nonlinear and quantum optical phenomena,^5–7 photochemical rates,^8–14 thermally activated ground-state chemical reactions under vibrational SC,^15–17 and exciton transport.^18,19

Theoretical analysis of these phenomena is non-trivial as polaritons reside at the intersection between quantum optics, quantum chemistry, and solid state physics. While quantum optical approaches such as Jaynes–Cummings or Dicke models have been used extensively,^2–4,20,21 they are ill-suited for describing the material aspects. This has spurred the development of advanced theoretical techniques in recent years based on various quantum optical and quantum chemistry methods.^{9–11,13,14,22–28} While those methods that provide an accurate account of quantum chemical effects are computationally very demanding, methods with a more approximate treatment of the materials’ degrees of freedom are limited with regard to chemical specificity. At the same time, ensemble effects and the collective interaction of molecules and/or nanoparticles (NPs) are of immediate experimental interest,^15,29 calling for methods that allow one to bridge between system specific predictions and computational efficiency.

We have recently demonstrated the usefulness of density functional theory (DFT)^30,31 and time-dependent DFT (TDDFT)³² calculations for studying polariton physics in NP–molecule systems.²⁷ Here, building on this study, we demonstrate that the observed effects can be reproduced nearly quantitatively using a subsystem approach where the different units, specifically NPs and molecules, only interact via dipolar coupling (DC). This approach allows one to very quickly evaluate the coupled response for a wide range of geometries at a computational cost that is orders of magnitudes smaller than a full TDDFT approach. In addition, it enables one to combine response function calculations from a variety of sources, including, e.g., different first-principles techniques, exchange-correlation functionals, and classical electrodynamic calculations.

In Sec. II, we provide an overview of the methodology used in this study as well as computational details. We then consider different scenarios, such as (i) coupling as a function of NP–molecule distance (Sec. III B), (ii) coupling as a function of the number of molecules and the size of the NP (Sec. III C), and (iii) mixing dynamic polarizabilities from different sources and/or types of calculations (Sec. III D). Finally, we summarize the conclusions and provide an outlook with respect to anticipated improvements and extensions.

In the dipolar coupling limit, the response of a system to an electrical field can be described by its dynamic polarizability tensor α_μν(ω). In Secs. II A and II B, we review the dynamic polarizability and the coupled response of a set of polarizable units within this framework.

The dynamic polarizability α(t) is a linear response function that defines the induced dipole moment d(t) in response to the external field E_ext(t),

Here and in the following, μ and ν refer to the Cartesian coordinates of the response and external electric field, respectively. We emphasize that α(t) is a causal response function that is zero in the negative time domain.

In frequency space, the dipole response is obtained via Fourier transformation and convolution theorem as

Hence, the α_μν component of the polarizability tensor describes the strength of the response along μ given a field along ν.

In practice, a typical calculation of the polarizability tensor within real-time TDDFT (RT-TDDFT) consists of recording the induced time-dependent dipole moment following a weak electric field impulse.³³ Then, the response in frequency space is obtained by Fourier transformation,

By choosing the external electric field $Eνext(t)=K0δ(t)$ to be aligned along ν, the μν component of the polarizability tensor is obtained as

The full polarizability tensor can thus be calculated by carrying out at most three calculations with external fields along each Cartesian direction; this number can be reduced further in the presence of symmetry. Here, the strength of the impulse, K₀, is set to be sufficiently weak so that the full response given by RT-TDDFT is dominated by the linear response regime.

The time-propagation approach allows calculating metallic nanoparticles with hundreds of atoms³⁴ in particular when combined with exchange correlation (XC) functionals that can accurately account for the d-band position in noble metal NPs.

For completeness, we note that for molecules it is usually more convenient to evaluate the excitation spectrum directly in the frequency domain using the Casida approach for linear-response TDDFT (LR-TDDFT).^35,36 In this case, the response can be formally obtained by solving a matrix eigenvalue equation, which yields the excitation energies ω_I and the corresponding transition dipole moments ⟨Ψ_I|r_μ|Ψ₀⟩, and the dynamic polarizability can then be obtained as

To deal with the divergence along the real frequency axis and to obtain (artificial) broadening of the spectrum, it is common to set ω → ω + iη in both Eqs. (3) and (5), resulting in a Lorentzian line shape with a peak width that is determined by the η parameter.

Below, we also consider coupling to a homogeneous sphere described by a local dielectric function ε(ω) according to Mie theory.³⁷ To this end, we note that in the quasi-static limit the dynamic polarizability tensor of a sphere with volume V is diagonal and takes the form³⁸ 

Given the dynamic polarizability of two or more units, we can evaluate the response of the coupled system. This is a widely known approach for molecular assemblies (see, for example, Refs. 39–41), but we present it here in whole for completeness.

Let $α0(i)(ω)$ be the irreducible polarizability tensor of the individual units. For a system composed of N units, the induced dipole moments, given the total electric field $Etot(i)(ω)$ at each unit, are obtained as

The Coulomb interaction between units i and j, located at $R→i$ and $R→j$⁠, is given in atomic units by $1/|R→ij|$⁠, where $R→ij=Ri−R→j$⁠. The dipole–dipole interaction between point dipoles at $R→i$ and $R→j$ is given by a tensor

Therefore, the total electric field at each polarizable unit is obtained as

By substituting Eq. (9) into Eq. (7) and solving for the induced dipole moment, we obtain

where the reducible unit-wise polarizability tensor is identified,

Each unit contributes to the total dipole moment of the coupled system. Assuming a uniform external electric field throughout the coupled system [$E→ext(1)=Eext(2)=⋯=Eext(N)$], the total polarizability tensor for a coupled system comprising N units is obtained by carrying out a double summation over all units of the system,

To present the results, we use the dipole strength function given by the imaginary part of the dynamic polarizability,

The dipole strength function equals the electronic photoabsorption spectrum safe for a constant multiplier and satisfies the Thomas–Reiche–Kuhn sum rule analogously to the oscillator strength.

The coupled systems considered in this study comprise one Al NP with 201, 586, or 1289 atoms and one or more benzene molecules, which have been investigated in whole using RT-TDDFT calculations in Ref. 27. The formalism outlined in Sec. II B requires material specific inputs in the form of the dynamic polarizability of the individual subsystems, i.e., isolated Al NPs of different sizes and an isolated benzene molecule, respectively.

The data for these systems have been calculated in Refs. 27 and 42 using the PBE XC functional⁴³ in the adiabatic limit and RT-TDDFT via the δ-kick technique³³ (Sec. II A), as implemented using linear combination of atomic orbital (LCAO) basis sets³⁴ in the GPAW code.⁴⁴ The projector augmented-wave⁴⁵ method was employed with double-ζ polarized (dzp) basis sets, as provided in GPAW. The wave functions were propagated up to 30 fs using a time step of 15 as. Further details of these calculations are given in Refs. 27 and 42.

For benzene, we also evaluated the first 16 roots of the excitation spectrum using LR-TDDFT within the Casida approach,^35,36 as implemented in the NWChem suite.⁴⁶ Calculations were carried out using the B3LYP functional^47,48 and the 6-311G^* basis set.^49,50 The excitation spectra were subsequently transformed into dynamic polarizabilities via Eq. (5). The convergence with respect to the number of roots is demonstrated in Figs. 1 and 2 of the supplementary material.

All spectra obtained in this work or taken from Ref. 42 are broadened using η = 0.1 eV. For the purpose of extracting coupling strength parameters, a coupled oscillator model⁵¹ was fitted to the obtained spectra. The details of the fitting scheme are outlined in Note 1 of the supplementary material.

In the following, we consider Al NPs with 201, 586, and 1289 atoms, which exhibit a plasmon resonance at about 7.7 eV [Fig. 1(a)]. Below, we analyze their coupling to two or more benzene molecules. The response of the latter is obtained either at the PBE or the B3LYP level yielding the first bright excitation at 7.1 eV and 7.2 eV, respectively [Fig. 1(b)].

Using the dynamic polarizabilities of the isolated units in Eq. (12), we first inspect the optical spectrum of a system comprising an Al₂₀₁ NP and two benzene molecules as a function of the NP–molecule distance. Reference data from RT-TDDFT calculations obtained using the PBE XC functional for the full system are available from Refs. 27 and 42. These full TDDFT calculations take into account not only coupling to all orders but also charge transfer and renormalization of the underlying states and excitations due to NP–molecule interactions.

The spectra obtained in the DC approximation and from full TDDFT calculations are in good agreement over the entire range of distances considered here [Fig. 2(a)]. Since the DC calculations only include dipolar interactions, the good agreement suggests that for the system under study, higher-order terms, charge transfer, and orbital hybridization play a rather small role in the part of the configuration space considered here. This need not hold for all systems, especially for NPs in direct contact, for which charge transfer and higher-order multipoles are more likely to become prominent. The importance of these effects should hence be assessed on a case-by-case basis.

For the shortest distance of 3 Å, a slight underestimation becomes apparent of both the position and width of the lower polariton, which we attribute to the absence of orbital overlap and to a lesser extent the neglect of higher-order multipole interactions. In the full TDDFT calculation, the former contribution, which is the more dominant one owing to the fact that higher order multipoles should be even more blue detuned from the benzene transition than the dipolar one, results in a broadening of the benzene transition, which couples with rate g to the Al NP plasmon. An increase in the decay rate (γ) of one of the strongly coupled elements causes the Rabi splitting Ω to increase as $Ω=4g2−(γpl−γex)2$ (for zero detuning),³ causing the TDDFT-computed polaritons to be broader and farther apart than in the calculations using the DC approximation. In addition, since the Al-NP plasmon is blue-detuned from the benzene transition, the two subsystems do not contribute equally to both polaritons. The lower polariton has a more benzene-like character and is more influenced by the width of the benzene transition. Conversely, the upper polariton is more Al-plasmon-like and is not influenced significantly by the decay rate of the molecular transition.

Despite the limitations discussed in the previous paragraph, the clear strength and benefit of the DC calculations lies in their very rapid computation. Indeed, once the response functions of the individual units have been computed, it is straightforward (and orders of magnitude faster than with full TDDFT modeling) to map out the configuration space. The DC calculations of hybrid systems can even have certain advantages compared to full TDDFT calculations using incomplete basis sets, as numerical factors that can affect the accuracy of the latter are effectively avoided.

Following the evolution of the spectrum as a practically continuous function of the NP–molecule distance reveals the emergence of a clear lower polariton state starting at distances of about 10 Å [Fig. 2(b)]. The coupling strengths g extracted from the DC calculations agree well with the full TDDFT calculations, both with respect to magnitude and distance dependence [Fig. 2(c)]. Interestingly, one observes the difference in g between the two types of calculations to increase with distance, which might appear unintuitive at first. This increasing difference stems from the different sources of error in the TDDFT and DC calculations. In the former, the use of localized basis sets that shift in space between calculations induces numerical errors, which exaggerate the blue-shift of the lower polariton at large distances (Fig. 3 of the supplementary material). In the latter, the missing effect of state hybridization leads to overly pronounced features, i.e., the polariton peaks are higher and the valley between them is deeper. In the coupled oscillator model, an increased value of g makes the lower polariton peak more pronounced and the separation between lower and upper polaritons larger.

Having confirmed the ability of the DC approximation to reproduce the distance dependence for an Al NP with two benzene molecules and the concomitant emergence of strong coupling, it is now instructive to analyze the spectrum as a function of the number N of benzene molecules, which provides another means for tuning the coupling. As their number increases, the effective mass of benzene excitations in the NP–molecule hybrid system increases, enhancing the coupling rate proportionally to the square root of N. The resulting large coupling strengths induce much larger changes of the polariton spectra than when the NP–molecule distance is varied. In particular, one notes that both the lower and upper polaritons exhibit significant changes with the Rabi splitting growing monotonically with N.²⁷ 

We consider Al NPs with three sizes containing 201, 586, and 1289 atoms coupled to up to eight benzene molecules. Again the DC approximation works well overall, although the agreement becomes worse as NP size and/or the number of benzene molecules increases (Fig. 3).

For the Al₂₀₁ case with several benzene molecules, one observes a gradual shift of spectral weight from the upper to the lower polariton. For small N, it is sensible to attribute the lower (upper) polariton primarily benzene (Al-NP) character. For larger N, however, the distinction between plasmonic (Al-NP) and excitonic characters (benzene molecule) becomes invalid as the system forms a fully mixed state [Fig. 3(a)]. This marked change between the balance of the upper and lower polaritons with N is the result of a gradual shift from a red-detuned molecular transition of one benzene with respect to the NP–plasmon to a blue-detuned one for N = 8.²⁷ This change is mostly ascribed to the red-shift of the NP–plasmon with increasing N, which may not be captured quantitatively in the DC approximation. The position of the lower polariton in this Al₂₀₁ sequence is slightly blue-shifted relative to the full TDDFT data. On the other hand, the DC calculations yield rather sensible results for the width of both polaritons.

By comparison, for the larger Al₅₈₆ and Al₁₂₈₉ NPs, one observes a systematic underestimation of both the position and width of the lower polariton [Figs. 3(b) and 3(c)]. As noted above (Sec. III B), the differences between DC and TDDFT calculations at a distance of 3 Å are likely related to the absence of orbital hybridization in the DC approximation. The orbital hybridization in TDDFT may also be the cause of other, potentially small, effects that escape the DC approximation. In addition to the above-mentioned red-shift of the NP–plasmon, the presence of the molecule may itself modify the mode volume of the cavity (beyond the explicitly studied plasmon-transition coupling) and increase the coupling strength.⁵² 

While the full TDDFT calculations can capture a variety of contributions that are missed by DC calculations, they are effectively limited by the need to describe the system using one common XC functional. Yet, functionals that perform well for metals are often poorly suited for molecular systems and vice versa. A distinct advantage of the DC approximation is its ability to combine response functions from multiple different sources, including but not limited to different types of electronic structure calculations as well as classical electrodynamics simulations. Here, we specifically exploit this feature to analyze the effect of the XC functional with regard to the description of the excitation spectrum of benzene. Simultaneously, we highlight the possibility for inter-code coupling when using the DC approximation.

We consider a system comprised of a Al₂₀₁ NP and two benzene molecules at a distance of 3 Å using dynamic polarizabilities calculated using different XC functionals (PBE and B3LYP) and different electronic structure codes (GPAW and NWChem) (Fig. 4). In comparison to PBE, the B3LYP XC functional gives the first bright excitation higher in energy, closer to both the experimental value and the plasmon resonance of the Al NP. The larger spectral overlap leads to less detuning, which primarily manifests itself in a more pronounced transfer of oscillator strength from the upper to the lower polariton. As proof-of-concept, we also demonstrate the DC approximation using the polarizability of a Mie sphere with the dielectric function of bulk Al by McPeak et al.⁵³ By matching the diameter of the Mie sphere to the distance between opposing {100} facets of the Al₂₀₁ NP, 16.2 Å, the spectra are similar. We note that the deviation between the Mie spectra based on the dielectric functions measured experimentally by McPeak et al. and Rakić⁵⁴ is considerably larger than the deviation between the present calculation and the spectrum based on the data from McPeak et al.

In this study, we have analyzed the efficacy of the DC approximation for capturing the emergence of SC in Al NP–benzene hybrid systems. To this end, we compared optical spectra obtained from DC calculations with the results from an earlier study that applied TDDFT to the entire system.²⁷ We find that overall the DC approximation is able to reproduce the full TDDFT results for this system well over a large size range of considered NPs. Deviations become more pronounced at short NP–molecule distances as orbital overlap and higher-order multipole interactions start to play a role.

After computation of the dynamic polarizabilities of the isolated components, the DC approximation is many orders of magnitude faster than a full scale TDDFT. This enables not only rapid screening of configuration space but computations for much larger systems than those accessible by fully fledged TDDFT calculations on complete systems (or other first-principles approaches) while still maintaining the underlying accuracy of the individual coupled elements. It thereby becomes possible to study, e.g., ensembles of multiple NPs, mixtures of multiple NPs interacting with multiple molecules or aggregates, and—by extension of the formalism—also the properties of such systems inside of cavities.

The approach employed here has a long history primarily in the context of molecular aggregates.^39,40 Here, we have demonstrated that it is equally applicable to NP and NP–molecule assemblies. More importantly, it is shown that by using dynamic polarizabilities from first-principles calculations, near-quantitative predictions become possible also for strongly coupled systems.

We note that calculations via the DC approximation are complementary to electrodynamics simulations using either the finite-difference time-domain method or the discrete dipole approximation, where the latter bears some technical similarities to the present approach. Specifically, DC calculations enable one to address systems with units in the nanometer size range that are commonly difficult or impossible to access using classical electrodynamics.

As noted, the DC approximation as used here has shortcomings that need to be addressed in the future. In particular, it will break at short NP–NP or NP–molecule distances, where higher-order multipoles and charge transfer effects become important. In particular, the latter require methods such as TDDFT calculations of the full system. One should also note that particles in solution are usually protected by a ligand shell, which prevents very short particle–particle distances. In the future, the approach could therefore be extended to include higher-order multipoles and to account for ligand shells around particles. Since systems of interest often are composed of particles in solution, the dielectric screening due to the solvation medium could also be included in the approach, which can be accomplished via polarizable continuous medium⁵⁵ approaches commonly used in quantum chemistry methods. As a final note, it is also of interest to include retardation effects, to accurately describe large systems, and to allow for periodic systems.

See the supplementary material for a convergence study with respect to the number of excitations in Eq. (5) as well as details regarding the fitting of the spectra to the coupled oscillator model.

We gratefully acknowledge the Knut and Alice Wallenberg Foundation (Grant No. 2019.0140, J.F., P.E.), the Swedish Research Council (Grant No. 2015-04153, J.F., P.E.), the Academy of Finland (Grant Nos. 332429, T.P.R; 295602, M.K.), and the Polish National Science Center (Grant No. 2019/34/E/ST3/00359, T.J.A.). The computations were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at NSC, C3SE, and PDC partially funded by the Swedish Research Council through Grant Agreement No. 2018-05973 as well as by the CSC—IT Center for Science, Finland, by the Aalto Science-IT project, Aalto University School of Science, and by the Interdisciplinary Center for Mathematical and Computational Modeling, University of Warsaw (Grant No. G55-6).

This study used the TDDFT data published in Ref. 42. The new data generated in this study are openly available via Zenodo at http://doi.org/10.5281/zenodo.4501057, Ref. 56.

The GPAW package^44,57 with LCAO basis sets⁵⁸ and the LCAO-RT-TDDFT implementation³⁴ was used for the RT-TDDFT calculations. The NWChem suite⁴⁶ was used for Casida LR-TDDFT calculations. The ASE library⁵⁹ was used for constructing and manipulating atomic structures. The NumPy,⁶⁰ SciPy,⁶¹ and Matplotlib⁶² Python packages and the VMD software^63,64 were used for processing and plotting data. The emcee⁶⁵ Python package was used to fit spectra to the coupled oscillator model.