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

Strong light-matter interactions facilitate not only emerging applications in quantum and non-linear optics but also modifications of materials properties. In particular the latter possibility has spurred the development of advanced theoretical techniques that can accurately capture both quantum optical and quantum chemical degrees of freedom. These methods are, however, computationally very demanding, which limits their application range. Here, we demonstrate that the optical spectra of nanoparticle-molecule assemblies, including strong coupling effects, can be predicted with good accuracy using a subsystem approach, in which the response functions of the different units are coupled only at the dipolar level. We demonstrate this approach by comparison with previous time-dependent density functional theory calculations for fully coupled systems of Al nanoparticles and benzene molecules. While the present study only considers few-particle systems, the approach can be readily extended to much larger systems and to include explicit optical-cavity modes.


I. INTRODUCTION
The coupling of light and matter in the strong coupling (SC) regime leads to the emergence of excited states of mixed nature, 1 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][9][10][11][12][13][14] thermally-activated ground-state chemical reactions under vibrational SC [15][16][17] as well as exciton transport. 18,19 Theoretical analysis of these phenomena is nontrivial 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][3][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][10][11]13,14,[22][23][24][25][26][27][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, e.g., 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) 32 calculations for studying polariton a) Corresponding author: erhart@chalmers.se physics in NP-molecule systems. 27 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 the following section, we provide an overview of the methodology used in this study as well as computational details. We then consider different scenarios, including (i ) coupling as a function of NP-molecule distance (Sect. III B), (ii ) coupling as a function of the number of molecules and the size of the NP (Sect. III C), and (iii ) mixing dynamic polarizabilities from different sources and/or types of calculations (Sect. III D). Finally, we summarize the conclusions and provide an outlook with respect to anticipated improvements and extensions.

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

A. Dynamic polarizability
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. 33 Then, the response in frequency space is obtained by Fourier transformation By choosing the external electric field E ext ν (t) = K 0 δ(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 0 , 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 34 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 µ |Ψ 0 , 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η both in 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. 37 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 38

B. Dipolar coupling
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 α (i) 0 (ω) 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 E (i) tot (ω) 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 = R i − 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 ), the total polarizability tensor for 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.

C. Computational details
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 the previous section requires material specific input 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 43 in the adiabatic limit and RT-TDDFT via the δ-kick technique 33 (Sect. II A) as implemented using linear combination of atomic orbitals (LCAO) basis sets 34 in the gpaw code. 44 The projector augmented-wave 45 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. 46 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).
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 51 was fitted to the obtained spectra. The details of the fitting scheme are outlined in Supplementary Note 1.

A. Dynamic polarizability of individual NPs and molecules
In the following we consider Al NPs with 201, 586 and 1289 atoms, which exhibit a plasmon resonance at about 7.7 eV (Fig. 1a). 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.3 eV, respectively (Fig. 1b).

B. Coupling as a function of distance
Using the dynamic polarizabilities of the isolated units in Eq. (12), we first inspect the optical spectrum of a system comprising an Al 201 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. 2a). 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.
For the shortest distance of 3Å a slight underestimation becomes apparent of both 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 of the decay rate (γ) of one of the strongly coupled elements causes the Rabi splitting Ω to increase as Ω = 4g 2 − (γ pl − γ ex ) 2 (for zero detuning), 3 causing the TDDFT-computed polaritons to be broader and farther apart than in the calculations using the DC approximation. Also, 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 benzenelike 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.
In spite of the limitations discussed in the previous paragraph, the clear strength and benefit of the DC cal-culations 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. 2b). 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. 2c). 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 ( Supplementary Fig. 1). 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 deeper. In the coupled oscillator model, an increased value g both makes the lower polariton peak more pronounced and the separation between lower and upper polariton larger.

C. Coupling as a function of the number of molecules and NP size
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 . 27 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 201 case with several benzene molecules, one observes a gradual shift of spectral weight from the up- per 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 character (benzene molecule) becomes invalid as the system forms a fully-mixed state (Fig. 3a). 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. 27 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 201 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 586 and Al 1289 NPs one observes a systematic underestimation of both position and width of the lower polariton (Fig. 3b,c). As noted above (Sect. 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. Optical response of a system comprising Al201 and two benzene molecules at a distance of 3Å, where the dynamic polarizability of the former is alternatively either by using the PBE XC functional and RT-TDDFT within the gpaw code or by using the Mie approximation in the quasi-static limit with an experimental dielectric function. The dynamic polarizability of the latter is alternatively obtained by using the PBE XC functional and RT-TDDFT within the gpaw code or by using the B3LYP XC functional and Casida LR-TDDFT within the NWChem code. The response of the isolated systems is shown for reference at the bottom of the figure using shaded lines.

D. Response functions from multiple sources
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 201 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 polari-ton. 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. 53 By matching the diameter of the Mie sphere to the distance between opposing {100} facets of the Al 201 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 Rakic 54 is considerably larger than the deviation between the present calculation and the spectrum based on the data from McPeak et al.

IV. CONCLUSIONS AND OUTLOOK
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. 27 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 firstprinciples 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 dis-tances, where higher-order multipoles and charge transfer effects become important. One should note, however, 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 55 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.

SUPPLEMENTARY MATERIAL
See supplementary material for details on the fitting of the spectra to the coupled oscillator model.

SOFTWARE USED
The gpaw package 44,57 with LCAO basis sets 58 and the LCAO-RT-TDDFT implementation 34 was used for the RT-TDDFT calculations. The NWChem suite 46 was used for Casida LR-TDDFT calculations. The ase library 59 was used for constructing and manipulating atomic structures.
The NumPy, 60 SciPy 61 and Matplotlib 62 Python packages and the VMD software 63,64 were used for processing and plotting data. The emcee 65 Python package was used to fit spectra to the coupled oscillator model.