Dipolar coupling of nanoparticle-molecule assemblies: An efficient approach for studying strong coupling
Jakub Fojt, Tuomas P. Rossi, Tomasz J. Antosiewicz, Mikael Kuisma, Paul Erhart
DDipolar coupling of nanoparticle-molecule assemblies: An efficient approach for studying strong coupling
Jakub Fojt, Tuomas P. Rossi, Tomasz J. Antosiewicz, Mikael Kuisma, and Paul Erhart a)1) Department of Physics, Chalmers University of Technology, SE-412 96 Gothenburg, Sweden Department of Applied Physics, Aalto University, 00076 Aalto, Finland Faculty of Physics, University of Warsaw, 02-093 Warsaw, Poland Department of Chemistry, University of Jyv¨askyl¨a, 40014 Jyv¨askyl¨a, Finland
Strong light-matter interactions facilitate not only emerging applications in quantum and non-linear opticsbut also modifications of materials properties. In particular the latter possibility has spurred the developmentof advanced theoretical techniques that can accurately capture both quantum optical and quantum chemicaldegrees of freedom. These methods are, however, computationally very demanding, which limits their appli-cation range. Here, we demonstrate that the optical spectra of nanoparticle-molecule assemblies, includingstrong coupling effects, can be predicted with good accuracy using a subsystem approach, in which the re-sponse functions of the different units are coupled only at the dipolar level. We demonstrate this approachby comparison with previous time-dependent density functional theory calculations for fully coupled systemsof Al nanoparticles and benzene molecules. While the present study only considers few-particle systems, theapproach 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 ofmixed nature, which are characterized by a coherent en-ergy exchange between the subsystems at a rate that ismuch faster than the respective damping rates. Emit-ter and the cavity thus form a light-matter hybrid (po-lariton) with modified and tunable properties, includ-ing nonlinear and quantum optical phenomena, pho-tochemical rates, thermally-activated ground-statechemical reactions under vibrational SC as well asexciton transport. Theoretical analysis of these phenomena is non-trivial as polaritons reside at the intersection betweenquantum optics, quantum chemistry, and solid statephysics. While quantum optical approaches such asJaynes–Cummings or Dicke models have been usedextensively, they are ill-suited for describing thematerial aspects. This has spurred the development ofadvanced theoretical techniques in recent years basedon various quantum optical and quantum chemistrymethods.
While those methods that pro-vide an accurate account of quantum chemical effects arecomputationally very demanding, methods with a moreapproximate treatment of the materials degrees of free-dom are limited with regard to chemical specificity. Atthe same time, ensemble effects and the collective interac-tion of, e.g., molecules and/or nanoparticles (NPs) are ofimmediate experimental interest, calling for methodsthat allow one to bridge between system specific predic-tions and computational efficiency.We have recently demonstrated the usefulness of den-sity functional theory (DFT) and time-dependentDFT (TDDFT) calculations for studying polariton a) Corresponding author: [email protected] physics in NP-molecule systems. Here, building on thisstudy, we demonstrate that the observed effects can bereproduced nearly quantitatively using a subsystem ap-proach where the different units, specifically NPs andmolecules, only interact via dipolar coupling (DC). Thisapproach allows one to very quickly evaluate the cou-pled response for a wide range of geometries at a compu-tational cost that is orders of magnitudes smaller thana full TDDFT approach. In addition, it enables oneto combine response function calculations from a vari-ety of sources, including, e.g., different first-principlestechniques, exchange-correlation functionals, and classi-cal electrodynamic calculations.In the following section, we provide an overview ofthe methodology used in this study as well as compu-tational details. We then consider different scenarios,including ( i ) coupling as a function of NP-molecule dis-tance (Sect. III B), ( ii ) coupling as a function of the num-ber of molecules and the size of the NP (Sect. III C),and ( iii ) mixing dynamic polarizabilities from differentsources and/or types of calculations (Sect. III D). Fi-nally, we summarize the conclusions and provide an out-look with respect to anticipated improvements and ex-tensions. II. METHODOLOGY
In the dipolar coupling limit, the response of a sys-tem to an electrical field can be described by its dynamicpolarizability tensor α µν ( ω ). In the following sections,we review the dynamic polarizability and the coupled re-sponse of a set of polarizable units within this framework. A. Dynamic polarizability
The dynamic polarizability α ( t ) is a linear responsefunction that defines the induced dipole moment d ( t ) in1 a r X i v : . [ c ond - m a t . m e s - h a ll ] J a n esponse to the external field E ext ( t ), d µ ( t ) = (cid:90) ∞−∞ (cid:88) ν α µν ( t − τ ) E ext ν ( τ )d τ. (1)Here and in the following, µ and ν refer to the Cartesiancoordinates of the response and external electric field re-spectively. We emphasize that α ( t ) is a causal responsefunction that is zero in the negative time domain.In frequency space the dipole response is obtained viaFourier transformation and convolution theorem as d µ ( ω ) = (cid:88) ν α µν ( ω ) E ext ν ( ω ) . (2)Hence, the α µν component of the polarizability tensordescribes the strength of the response along µ given afield along ν .In practice, a typical calculation of the polarizabilitytensor within real-time TDDFT (RT-TDDFT) consists ofrecording the induced time-dependent dipole moment fol-lowing a weak electric field impulse. Then, the responsein frequency space is obtained by Fourier transformation d µ ( ω ) = (cid:90) ∞ d µ ( t ) e iωt d t. (3)By choosing the external electric field E ext ν ( t ) = K δ ( t )to be aligned along ν , the µν component of the polariz-ability tensor is obtained as α µν ( ω ) = d µ ( ω ) K . (4)The full polarizability tensor can thus be calculated bycarrying out at most three calculations with externalfields along each Cartesian direction; this number canbe 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 isdominated by the linear response regime.The time-propagation approach allows calculatingmetallic nanoparticles with hundreds of atoms in par-ticular when combined with exchange correlation (XC)functionals that can accurately account for the d -bandposition in noble metal NPs.For completeness, we note that for molecules it is usu-ally more convenient to evaluate the excitation spectrumdirectly in the frequency domain using the Casida ap-proach for linear-response TDDFT (LR-TDDFT). Inthis case the response can be formally obtained by solv-ing a matrix eigenvalue equation, which yields the excita-tion energies ω I and the corresponding transition dipolemoments (cid:104) Ψ I | r µ | Ψ (cid:105) , and the dynamic polarizability canthen be obtained as α µν ( ω ) = (cid:88) I ω I (cid:104) Ψ | r µ | Ψ I (cid:105) (cid:104) Ψ I | r ν | Ψ (cid:105) ω I − ω . (5)To deal with the divergence along the real frequencyaxis 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 widththat is determined by the η parameter.Below, we also consider coupling to a homogeneoussphere described by a local dielectric function ε ( ω ) ac-cording to Mie theory. To this end, we note that in thequasi-static limit the dynamic polarizability tensor of asphere with volume V is diagonal and takes the form α µν ( ω ) = 3 ε V ε ( ω ) − ε ( ω ) + 2 δ µν . (6) B. Dipolar coupling
Given the dynamic polarizability of two or more unitswe can evaluate the response of the coupled system. Thisis a widely known approach for molecular assemblies, seefor example Refs. 39–41, but we present it here in wholefor completeness.Let α ( i )0 ( ω ) be the irreducible polarizability tensor of the individual units. For a system composed of N units, theinduced dipole moments, given the total electric field E ( i )tot ( ω ) at each unit, are obtained as d (1) ( ω ) d (2) ( ω )... d ( N ) ( ω ) (cid:124) (cid:123)(cid:122) (cid:125) d ( ω ) = α (1)0 ( ω ) α (2)0 ( ω ) . . . α ( N )0 ( ω ) (cid:124) (cid:123)(cid:122) (cid:125) α ( ω ) E (1)tot ( ω ) E (2)tot ( ω )... E ( N )tot ( ω ) (cid:124) (cid:123)(cid:122) (cid:125) E tot ( ω ) . (7)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 T iµjν = ∇ µ ∇ ν | R ij | = δ µν | R ij | − R ij,µ R ij,ν | R ij | . (8)2herefore, the total electric field at each polarizable unit is obtained as E (1)tot ( ω ) E (2)tot ( ω )... E ( N )tot ( ω ) (cid:124) (cid:123)(cid:122) (cid:125) E tot ( ω ) = E (1)ext ( ω ) E (2)ext ( ω )... E ( N )ext ( ω ) (cid:124) (cid:123)(cid:122) (cid:125) E ext ( ω ) − T . . . T N T T N ... . . . ... T N T N . . . (cid:124) (cid:123)(cid:122) (cid:125) T d (1) ( ω ) d (2) ( ω )... d ( N ) ( ω ) (cid:124) (cid:123)(cid:122) (cid:125) d ( ω ) . (9)By substituting Eq. (9) into Eq. (7) and solving for theinduced dipole moment we obtain d ( ω ) = [ I + α ( ω ) T ] − α ( ω ) E ext ( ω ) , (10)where the reducible unit-wise polarizability tensor isidentified α ( ω ) = [ I + α ( ω ) T ] − α ( ω ) (11a)= (cid:2) α ( ω ) − + T (cid:3) − . (11b)Each unit contributes to the total dipole moment of thecoupled system. Assuming a uniform external electricfield throughout the coupled system ( E (1)ext = E (2)ext = . . . = E ( N )ext ), the total polarizability tensor for coupledsystem comprising N units is obtained by carrying out adouble summation over all units of the system α µν ( ω ) = N (cid:88) i N (cid:88) j [ α ] iµjν ( ω ) . (12)To present the results, we use the dipole strength func-tion given by the imaginary part of the dynamic polariz-ability S µ ( ω ) = 2 ωπ (cid:61) [ α µµ ( ω )] . (13)The dipole strength function equals the electronic pho-toabsorption spectrum safe for a constant multiplier andsatisfies the Thomas–Reiche–Kuhn sum rule analogouslyto the oscillator strength. C. Computational details
The coupled systems considered in this study compriseone Al NP with 201, 586 or 1289 atoms and one or morebenzene molecules, which have been investigated in wholeusing RT-TDDFT calculations in Ref. 27. The formalismoutlined in the previous section requires material specificinput in the form of the dynamic polarizability of theindividual subsystems, i.e. isolated Al NPs of differentsizes and an isolated benzene molecule, respectively.The data for these systems have been calculated inRefs. 27 and 42 using the PBE XC functional in the adi-abatic limit and RT-TDDFT via the δ -kick technique (Sect. II A) as implemented using linear combination ofatomic orbitals (LCAO) basis sets in the gpaw code. The projector augmented-wave method was employedwith double- ζ polarized (dzp) basis sets as provided in gpaw . The wave functions were propagated up to 30 fsusing a time step of 15 as. Further details of these calcu-lations are given in Refs. 27 and 42.For benzene we also evaluated the first 16 roots ofthe excitation spectrum using LR-TDDFT within theCasida approach as implemented in the NWChem suite. Calculations were carried out using the B3LYPfunctional and the 6-311G ∗ basis set. The excita-tion spectra were subsequently transformed into dynamicpolarizabilities via Eq. (5).All spectra obtained in this work or taken from Ref. 42are broadened using η = 0 . was fitted to the obtained spectra. Thedetails of the fitting scheme are outlined in Supplemen-tary Note 1. III. RESULTS AND DISCUSSIONA. Dynamic polarizability of individual NPsand molecules
In the following we consider Al NPs with 201, 586 and1289 atoms, which exhibit a plasmon resonance at about7.7 eV (Fig. 1a). Below we analyze their coupling to twoor more benzene molecules. The response of the latter isobtained either at the PBE or the B3LYP level yieldingthe first bright excitation at 7.1 eV and 7.3 eV, respec-tively (Fig. 1b).
B. Coupling as a function of distance
Using the dynamic polarizabilities of the isolated unitsin Eq. (12), we first inspect the optical spectrum of a sys-tem comprising an Al
NP and two benzene moleculesas a function of the NP-molecule distance. Referencedata from RT-TDDFT calculations obtained using thePBE XC functional for the full system are available fromRefs. 27 and 42. These full TDDFT calculations take intoaccount not only coupling to all orders but also charge3 (a) Al Al Al Energy (eV) (b)
PBE, GPAWB3LYP, NWC
HEM D i po l e s t r eng t h f un c t i on ( / e V ) FIG. 1. Optical spectra for isolated (a) Al NPs and (b)benzene molecules. All NP data in (a) were obtained from gpaw (RT-TDDFT) using the PBE XC functional. The ben-zene data in (b) were obtained from gpaw (RT-TDDFT) and
NWChem (Casida) calculations using different XC function-als as indicated in the legend. transfer and renormalization of the underlying states andexcitations due to NP-molecule interactions.The spectra obtained in the DC approximation andfrom full TDDFT calculations are in good agreement overthe entire range of distances considered here (Fig. 2a).Since the DC calculations only include dipolar interac-tions, the good agreement suggests that for the systemunder study, higher-order terms, charge transfer, and or-bital hybridization play a rather small role in the part ofthe configuration space considered here.For the shortest distance of 3 ˚A a slight underestima-tion becomes apparent of both position and width ofthe lower polariton, which we attribute to the absenceof orbital overlap and to a lesser extent the neglect ofhigher-order multipole interactions. In the full TDDFTcalculation, the former contribution, which is the moredominant one owing to the fact that higher order multi-poles should be even more blue detuned from the benzenetransition than the dipolar one, results in a broadeningof the benzene transition, which couples with rate g tothe Al NP plasmon. An increase of the decay rate ( γ )of one of the strongly coupled elements causes the Rabisplitting Ω to increase as Ω = (cid:112) g − ( γ pl − γ ex ) (forzero detuning), causing the TDDFT-computed polari-tons to be broader and farther apart than in the cal-culations using the DC approximation. Also, since theAl-NP plasmon is blue-detuned from the benzene tran-sition, the two subsystems do not contribute equally toboth polaritons. The lower polariton has a more benzene-like character and is more influenced by the width of thebenzene transition. Conversely, the upper polariton ismore Al-plasmon-like and is not influenced significantlyby the decay rate of the molecular transition.In spite of the limitations discussed in the previousparagraph, 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 havebeen computed, it is straightforward (and orders of mag-nitude faster than with full TDDFT modeling) to mapout the configuration space. The DC calculations of hy-brid systems can even have certain advantages comparedto full TDDFT calculations using incomplete basis sets,as numerical factors that can affect the accuracy of thelatter are effectively avoided.Following the evolution of the spectrum as a practicallycontinuous function of the NP-molecule distance revealsthe emergence of a clear lower polariton state starting atdistances of about 10 ˚A (Fig. 2b). The coupling strengths g extracted from the DC calculations agree well with thefull TDDFT calculations, both with respect to magni-tude and distance dependence (Fig. 2c). Interestinglyone observes the difference in g between the two types ofcalculations to increase with distance, which might ap-pear unintuitive at first. This increasing difference stemsfrom the different sources of error in the TDDFT and DCcalculations. In the former, the use of localized basis setsthat shift in space between calculations induces numer-ical errors, which exaggerate the blue-shift of the lowerpolariton at large distances (Supplementary Fig. 1). Inthe latter, the missing effect of state hybridization leadsto overly pronounced features, i.e. the polariton peaksare higher and the valley between them deeper. In thecoupled oscillator model, an increased value g both makesthe lower polariton peak more pronounced and the sep-aration between lower and upper polariton larger. C. Coupling as a function of the number ofmolecules and NP size
Having confirmed the ability of the DC approxima-tion to reproduce the distance dependence for an Al NPwith two benzene molecules and the concomitant emer-gence of strong coupling, it is now instructive to analyzethe spectrum as a function of the number N of benzenemolecules, which provides another means for tuning thecoupling. As their number increases the effective mass ofbenzene excitations in the NP-molecule hybrid system in-creases, enhancing the coupling rate proportionally to thesquare root of N . The resulting large coupling strengthsinduce much larger changes of the polariton spectra thanwhen the NP-molecule distance is varied. In particular,one notes that both the lower and upper polaritons ex-hibit significant changes with the Rabi splitting growingmonotonically with N . We consider Al NPs with three sizes containing 201,586 and 1289 atoms coupled to up to eight benzenemolecules. Again the DC approximation works welloverall, although the agreement becomes worse as NPsize and/or the number of benzene molecules increases(Fig. 3).For the Al case with several benzene molecules, oneobserves a gradual shift of spectral weight from the up-4
Energy (eV) D i po l e s t r eng t h f un c t i on ( / e V ) (a) Dipolar couplingFull TDDFT 5 10 15
Distance ( ˚A) E ne r g y ( e V ) (b) Distance ( ˚A) C oup li ng s t r eng t h g ( e V ) (c) Dipolar couplingFull TDDFT160180200220240260280 D i po l e s t r eng t h f un c t i on ( / e V ) FIG. 2. Distance dependence of the optical response from dipolar coupling and full TDDFT calculations. (a) Opticalspectra for a system comprising a Al
NP and two benzene molecules as the NP-molecule distance is varied from RT-TDDFTcalculations for the full system (Full TDDFT) and from Eq. (12) (Dipolar coupling). (b) Spectrum of Al +2 benzene system asa practically continuous function of the NP-molecule distance obtained using the dipolar coupling approximation. (c) Couplingstrength extracted from Bayesian fits of the spectra in (a) and (b) to a coupled oscillator model. The extracted values areshown as round markers for the full system (Full TDDFT) and as a solid line for the dipolar coupling calculations. Al x 0.5 (a) Al + N × Dipolar couplingFull TDDFT05001000 Al x 0.5 (b) Al + N × Dipolar couplingFull TDDFT5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5
Energy (eV) Al x 0.5 (c) Al + N × Dipolar couplingFull TDDFT0.0 0.2 0.4 0.6 0.8 1.00.00.20.40.60.81.0 D i po l e s t r eng t h f un c t i on ( / e V ) FIG. 3. Dependence of the optical response on the number ofmolecules from dipolar coupling and full TDDFT calculations.Optical spectra of a system comprising (a) a Al
NP, (b)a Al or (c) a Al and a variable number of benzenemolecules from RT-TDDFT calculations for the full system(Full TDDFT) and from Eq. (12) (Dipolar coupling). Thesolid gray lines show the spectrum of the isolated NPs forcomparison. per to the lower polariton. For small N it is sensible toattribute the lower (upper) polariton primarily benzene(Al-NP) character. For larger N , however, the distinc-tion between plasmonic (Al-NP) and excitonic character(benzene molecule) becomes invalid as the system forms afully-mixed state (Fig. 3a). This marked change betweenthe balance of the upper and lower polaritons with N isthe result of a gradual shift from a red-detuned moleculartransition of one benzene with respect to the NP-plasmonto a blue-detuned one for N = 8. This change is mostlyascribed to the red-shift of the NP-plasmon with increas-ing N , which may not be captured quantitatively in theDC-approximation. The position of the lower polaritonin this Al sequence is slightly blue-shifted relative tothe full TDDFT data. On the other hand, the DC calcu-lations yield rather sensible results for the width of bothpolaritons.By comparison, for the larger Al and Al NPsone observes a systematic underestimation of both po-sition and width of the lower polariton (Fig. 3b,c). Asnoted above (Sect. III B), the differences between DC andTDDFT calculations at a distance of 3 ˚A are likely relatedto the absence of orbital hybridization in the DC approx-imation. The orbital hybridization in TDDFT may alsobe the cause of other, potentially small, effects that es-cape the DC approximation. In addition to the abovementioned red-shift of the NP-plasmon, the presence ofthe molecule may itself modify the mode volume of thecavity (beyond the explicitly studied plasmon-transitioncoupling) and increase the coupling strength. Dipolar coupling
Full TDDFT (PBE) NP: PBE, GPAWNP: PBE, GPAWNP: Mie (experimental ε ): PBE, GPAW: B3LYP, NWChem: B3LYP, NWChem6.0 6.5 7.0 7.5 8.0 8.5 Energy (eV)
PBE, GPAWMie (experimental ε )x 10 D i po l e s t r eng t h f un c t i on ( / e V ) FIG. 4. Mixing response calculations from different sources.Optical response of a system comprising Al and two ben-zene molecules at a distance of 3 ˚A, where the dynamic polar-izability of the former is alternatively either by using the PBEXC functional and RT-TDDFT within the gpaw code or byusing the Mie approximation in the quasi-static limit with anexperimental dielectric function. The dynamic polarizabilityof the latter is alternatively obtained by using the PBE XCfunctional and RT-TDDFT within the gpaw code or by us-ing the B3LYP XC functional and Casida LR-TDDFT withinthe
NWChem code. The response of the isolated systems isshown for reference at the bottom of the figure using shadedlines.
D. Response functions from multiple sources
While the full TDDFT calculations can capture a vari-ety of contributions that are missed by DC calculations,they are effectively limited by the need to describe thesystem using one common XC functional. Yet, function-als that perform well for metals are often poorly suitedfor molecular systems and vice versa. A distinct advan-tage of the DC approximation is its ability to combineresponse functions from multiple different sources, in-cluding but not limited to different types of electronicstructure calculations as well as classical electrodynam-ics simulations. Here, we specifically exploit this featureto analyze the effect of the XC functional with regard tothe description of the excitation spectrum of benzene. Si-multaneously, we highlight the possibility for inter-codecoupling when using the DC approximation.We consider a system comprised of a Al
NP and twobenzene molecules at a distance of 3 ˚A using dynamicpolarizabilities calculated using different XC function-als (PBE and B3LYP) and different electronic structurecodes ( gpaw and
NWChem ) (Fig. 4). In comparisonto PBE, the B3LYP XC functional gives the first brightexcitation higher in energy, closer to both the experimen-tal value and the plasmon resonance of the Al NP. Thelarger spectral overlap leads to less detuning, which pri-marily manifests itself in a more pronounced transfer ofoscillator strength from the upper to the lower polari- ton. As proof-of-concept, we also demonstrate the DCapproximation using the polarizability of a Mie spherewith the dielectric function of bulk Al by McPeak et al . By matching the diameter of the Mie sphere to the dis-tance between opposing { } facets of the Al NP,16.2 ˚A, the spectra are similar. We note that the de-viation between the Mie spectra based on the dielectricfunctions measured experimentally by McPeak et al . andRakic is considerably larger than the deviation betweenthe present calculation and the spectrum based on thedata from McPeak et al. IV. CONCLUSIONS AND OUTLOOK
In this study we have analyzed the efficacy of the DCapproximation for capturing the emergence of SC in AlNP-benzene hybrid systems. To this end, we comparedoptical spectra obtained from DC calculations with theresults from an earlier study that applied TDDFT to theentire system. We find that overall the DC approxima-tion is able to reproduce the full TDDFT results for thissystem well over a large size range of considered NPs. De-viations become more pronounced at short NP-moleculedistances as orbital overlap and higher-order multipoleinteractions start to play a role.After computation of the dynamic polarizabilities ofthe isolated components, the DC approximation is manyorders of magnitude faster than a full scale TDDFT.This enables not only rapid screening of configurationspace but computations for much larger systems thanthose accessible by fully-fledged TDDFT calculations oncomplete systems (or other first-principles approaches),while still maintaining the underlying accuracy of the in-dividual coupled elements. It thereby becomes possibleto study, e.g., ensembles of multiple NPs, mixtures ofmultiple NPs interacting with multiple molecules or ag-gregates, and – by extension of the formalism – also theproperties of such systems inside of cavities.The approach employed here has a long history pri-marily in the context of molecular aggregates.
Here,we have demonstrated that it is equally applicable toNP and NP-molecule assemblies. More importantly it isshown that by using dynamic polarizabilities from first-principles calculations near-quantitative predictions be-come possible also for strongly coupled systems.We note that calculations via the DC approxima-tion are complementary to electrodynamics simulationsusing either the finite-difference time-domain methodor the discrete dipole approximation, where the latterbears some technical similarities to the present approach.Specifically, DC calculations enable one to address sys-tems with units in the nanometer size range that arecommonly difficult or impossible to access using classicalelectrodynamics.As noted the DC approximation as used here has short-comings that need to be addressed in the future. In par-ticular it will break at short NP-NP or NP-molecule dis-6ances, where higher-order multipoles and charge trans-fer effects become important. One should note, however,that particles in solution are usually protected by a lig-and shell, which prevents very short particle-particle dis-tances. In the future the approach could therefore beextended to include higher-order multipoles and to ac-count for ligand shells around particles. Since systems ofinterest often are composed of particles in solution, thedielectric screening due to the solvation medium couldalso be included in the approach, which can be accom-plished via polarizable continuous medium approachescommonly used in quantum chemistry methods. As afinal note, it is also of interest to include retardation ef-fects, to accurately describe large systems, and to allowfor periodic systems. SUPPLEMENTARY MATERIAL
See supplementary material for details on the fitting ofthe spectra to the coupled oscillator model.
ACKNOWLEDGMENTS
We gratefully acknowledge the Knut and Alice Wal-lenberg Foundation (2019.0140, J.F., P.E.), the SwedishResearch Council (2015-04153, J.F., P.E.), the Academyof Finland (332429, T.P.R; 295602, M.K.), and thePolish National Science Center (2019/34/E/ST3/00359,T.J.A.). The computations were enabled by resourcesprovided by the Swedish National Infrastructure forComputing (SNIC) at NSC, C3SE and PDC partiallyfunded by the Swedish Research Council through grantagreement no. 2018-05973 as well as by the CSC – ITCenter for Science, Finland, by the Aalto Science-ITproject, Aalto University School of Science, and by theInterdisciplinary Center for Mathematical and Computa-tional Modeling, University of Warsaw (Grant
DATA AVAILABILITY
This study used the TDDFT data published in Ref. 42.The new data generated in this study are openly avail-able via Zenodo at http://doi.org/10.5281/zenodo.4095510 , Ref. 56.
SOFTWARE USED
The gpaw package with LCAO basis sets andthe LCAO-RT-TDDFT implementation was used forthe RT-TDDFT calculations. The NWChem suite was used for Casida LR-TDDFT calculations. The ase library was used for constructing and manip-ulating atomic structures. The NumPy, SciPy and Matplotlib Python packages and the VMDsoftware were used for processing and plotting data.The emcee Python package was used to fit spectra tothe coupled oscillator model.
REFERENCES T. W. Ebbesen, Acc. Chem. Res. , 2403 (2016). P. T¨orm¨a and W. L. Barnes, Rep. Prog. Phys. ,013901 (2015). D. G. Baranov, M. Wers¨all, J. Cuadra, T. J. An-tosiewicz, and T. Shegai, ACS Photonics , 24 (2018). A. Frisk Kockum, A. Miranowicz, S. De Liberato,S. Savasta, and F. Nori, Nat. Rev. Phys. , 19 (2019). T. Yoshie, A. Scherer, J. Hendrickson, G. Khitrova,H. M. Gibbs, G. Rupper, C. Ell, O. B. Shchekin, andD. G. Deppe, Nature , 200 (2004). K. M. Birnbaum, A. Boca, R. Miller, A. D. Boozer,T. E. Northup, and H. J. Kimble, Nature , 87(2005). J. Kasprzak, M. Richard, S. Kundermann, A. Baas,P. Jeambrun, J. M. J. Keeling, F. M. Marchetti, M. H.Szymanska, R. Andre, J. L. Staehli, V. Savona, P. B.Littlewood, B. Deveaud, and L. S. Dang, Nature ,409 (2006). J. A. Hutchison, T. Schwartz, C. Genet, E. Devaux,and T. W. Ebbesen, Angew. Chem. Int. Ed. , 1592(2012). F. Herrera and F. C. Spano, Phys. Rev. Lett. ,238301 (2016). J. Galego, F. J. Garcia-Vidal, and J. Feist, Nat. Com-mun. , 13841 (2016). L. A. Mart´ınez-Mart´ınez, R. F. Ribeiro, J. Campos-Gonz´alez-Angulo, and J. Yuen-Zhou, ACS Photonics , 167 (2018). B. Munkhbat, M. Wers¨all, D. G. Baranov, T. J. An-tosiewicz, and T. Shegai, Sci. Adv. , eaas9552 (2018). J. Flick and P. Narang, J. Chem. Phys. , 094116(2020). J. Fregoni, G. Granucci, M. Persico, and S. Corni,Chem , 250 (2020). A. Thomas, J. George, A. Shalabney, M. Dryzhakov,S. J. Varma, J. Moran, T. Chervy, X. Zhong, E. De-vaux, C. Genet, J. A. Hutchison, and T. W. Ebbesen,Angew. Chem. Int. Ed. , 11462 (2016). A. Thomas, L. Lethuillier-Karl, K. Nagarajan, R. M. A.Vergauwe, J. George, T. Chervy, A. Shalabney, E. De-vaux, C. Genet, J. Moran, and T. W. Ebbesen, Science , 615 (2019). K. Hirai, R. Takeda, J. A. Hutchison, and H. Uji-i,Angew. Chem. Int. Ed. , 5332 (2020). J. Feist and F. J. Garcia-Vidal, Phys. Rev. Lett. ,196402 (2015). X. Zhong, T. Chervy, S. Wang, J. George, A. Thomas,J. A. Hutchison, E. Devaux, C. Genet, and T. W.Ebbesen, Angew. Chem. Int. Ed. , 6202 (2016).7 E. T. Jaynes and F. W. Cummings, Proc. IEEE , 89(1963). B. M. Garraway, Phil. Trans. R. Soc. A , 1137(2011). J. Galego, F. J. Garcia-Vidal, and J. Feist, Phys. Rev.X , 041022 (2015). M. Ruggenthaler, J. Flick, C. Pellegrini, H. Appel,I. V. Tokatly, and A. Rubio, Phys. Rev. A , 012508(2014). J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio,Proc. Natl. Acad. Sci. U.S.A. , 3026 (2017). H. Ling Luk, J. Feist, J. J. Toppari, and G. Groenhof,J. Chem. Theory Comput. , 4324 (2017). T. Neuman, R. Esteban, D. Casanova, F. J. Garcia-Vidal, and J. Aizpurua, Nano Lett. , 2358 (2018). T. P. Rossi, T. Shegai, P. Erhart, and T. J. An-tosiewicz, Nat. Commun. , 3336 (2019). J. Flick, D. M. Welakuh, M. Ruggenthaler, H. Appel,and A. Rubio, ACS Photonics , 2757 (2019). E. Orgiu, J. George, J. Hutchison, E. Devaux, J. Dayen,B. Doudin, F. Stellacci, C. Genet, J. Schachenmayer,C. Genes, G. Pupillo, P. Samori, and T. Ebbesen, Nat.Mater. , 1123 (2015). P. Hohenberg and W. Kohn, Phys. Rev. , B864(1964). W. Kohn and L. Sham, Phys. Rev. , A1133 (1965). E. Runge and E. Gross, Phys. Rev. Lett. , 997(1984). K. Yabana and G. F. Bertsch, Phys. Rev. B , 4484(1996). M. Kuisma, A. Sakko, T. Rossi, A. Larsen, J. Enko-vaara, L. Lehtovaara, and T. Rantala, Phys. Rev. B , 115431 (2015). M. E. Casida, in
Recent Advances in Density Func-tional Methods, Part I , edited by D. P. Chong (WorldScientific, Singapore, 1995) p. 155. M. E. Casida, J. Mol. Struct. THEOCHEM , 3(2009). G. Mie, Ann. Phys. , 377 (1908). A. Sihvola, J. Nanomat. , 1 (2007). H. DeVoe, J. Chem. Phys. , 393 (1964). H. Fidder, J. Knoester, and D. A. Wiersma, J. Chem.Phys. , 7880 (1991). B. Augui´e, B. L. Darby, and E. C. L. Ru, Nanoscale , 12177 (2019). T. P. Rossi, T. Shegai, P. Erhart, and T. J. An-tosiewicz, Zenodo (2019), 10.5281/zenodo.3242317. J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996). J. Enkovaara, C. Rostgaard, J. J. Mortensen, J. Chen,M. Du(cid:32)lak, L. Ferrighi, J. Gavnholt, C. Glinsvad,V. Haikola, H. A. Hansen, H. H. Kristoffersen,M. Kuisma, A. H. Larsen, L. Lehtovaara, M. Ljung-berg, O. Lopez-Acevedo, P. G. Moses, J. Ojanen,T. Olsen, V. Petzold, N. A. Romero, J. Stausholm-Møller, M. Strange, G. A. Tritsaris, M. Vanin, M. Wal-ter, B. Hammer, H. H¨akkinen, G. K. H. Madsen, R. M.Nieminen, J. K. Nørskov, M. Puska, T. T. Rantala, J. Schiøtz, K. S. Thygesen, and K. W. Jacobsen, J.Phys.: Condens. Matter , 253202 (2010). P. E. Bl¨ochl, Phys. Rev. B , 17953 (1994). M. Valiev, E. J. Bylaska, N. Govind, K. Kowalski, T. P.Straatsma, H. J. J. Van Dam, D. Wang, J. Nieplocha,E. APhysical Review A, T. L. Windus, and W. A.de Jong, Comput. Phys. Commun. , 1477 (2010). C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B ,785 (1988). A. D. Becke, J. Chem. Phys. , 5648 (1993). A. D. McLean and G. S. Chandler, J. Chem. Phys. ,5639 (1980). R. Krishnan, J. S. Binkley, R. Seeger, and J. A. Pople,J. Chem. Phys. , 650 (1980). X. Wu, S. K. Gray, and M. Pelton, Opt. Express ,23633 (2010). Z.-J. Yang, T. J. Antosiewicz, and T. Shegai, Opt.Express , 20373 (2016). K. M. McPeak, S. V. Jayanti, S. J. P. Kress, S. Meyer,S. Iotti, A. Rossinelli, and D. J. Norris, ACS Photonics , 326 (2015). A. D. Raki´c, Appl. Opt. , 4755 (1995). V. Barone, M. Cossi, and J. Tomasi, J. Chem. Phys. , 3210 (1997). J. Fojt, T. P. Rossi, T. J. Antosiewicz, M. Kuisma, andP. Erhart, Zenodo (2021), 10.5281/zenodo.4095510. J. Mortensen, L. Hansen, and K. Jacobsen, Phys. Rev.B , 035109 (2005). A. Larsen, M. Vanin, J. Mortensen, K. Thygesen, andK. Jacobsen, Phys. Rev. B , 195112 (2009). A. Larsen, J. Mortensen, J. Blomqvist, I. Castelli,R. Christensen, M. Dulak, J. Friis, M. Groves, B. Ham-mer, C. Hargus, E. Hermes, P. Jennings, P. Jensen,J. Kermode, J. Kitchin, E. Kolsbjerg, J. Kubal,K. Kaasbjerg, S. Lysgaard, J. Maronsson, T. Max-son, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard,J. S. O. Sch¨utt, M. Strange, K. Thygesen, T. Vegge,L. Vilhelmsen, M. Walter, Z. Zeng, and K. W. Jacob-sen, J. Phys.: Condens. Matter , 273002 (2017). C. R. Harris, K. J. Millman, S. J. van der Walt,R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser,J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus,S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Hal-dane, J. F. del R´ıo, M. Wiebe, P. Peterson, P. G´erard-Marchant, K. Sheppard, T. Reddy, W. Weckesser,H. Abbasi, C. Gohlke, and T. E. Oliphant, Nature , 357 (2020). P. Virtanen, R. Gommers, T. E. Oliphant, M. Haber-land, T. Reddy, D. Cournapeau, E. Burovski, P. Pe-terson, W. Weckesser, J. Bright, S. J. van der Walt,M. Brett, J. Wilson, K. J. Millman, N. Mayorov,A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J.Carey, ˙I. Polat, Y. Feng, E. W. Moore, J. VanderPlas,D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen,E. A. Quintero, C. R. Harris, A. M. Archibald, A. H.Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0Contributors, Nat. Methods , 261 (2020). J. D. Hunter, Comput. Sci. Eng. , 90 (2007).8 W. Humphrey, A. Dalke, and K. Schulten, J. Mol.Graph. , 33 (1996). J. Stone,
An Efficient Library for Parallel Ray Trac-ing and Animation , Master’s thesis, Computer ScienceDepartment, University of Missouri-Rolla (1998). D. Foreman-Mackey, D. W. Hogg, D. Lang, andJ. Goodman, Publ. Astron. Soc. Pac.125