Strong light–matter interactions facilitate not only emerging applications in quantum and non-linear optics but also modifications of properties of materials. 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 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.

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–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)32 calculations for studying polariton 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 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 Eext(t),

dμ(t)=ναμν(tτ)Eνext(τ)dτ.
(1)

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

dμ(ω)=ναμν(ω)Eνext(ω).
(2)

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,

dμ(ω)=0dμ(t)eiωtdt.
(3)

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

αμν(ω)=dμ(ω)K0.
(4)

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, K0, 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 atoms34 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

αμν(ω)=I2ωIΨ0|rμ|ΨIΨI|rν|Ψ0ωI2ω2.
(5)

To deal with the divergence along the real frequency axis and to obtain (artificial) broadening of the spectrum, it is common to set ωω + 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.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 form38 

αμν(ω)=3ε0Vε(ω)1ε(ω)+2δμν.
(6)

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

d(1)(ω)d(2)(ω)d(N)(ω)d̲(ω)=α0(1)(ω)α0(2)(ω)α0(N)(ω)α̲0(ω)Etot(1)(ω)Etot(2)(ω)Etot(N)(ω)E̲tot(ω).
(7)

The Coulomb interaction between units i and j, located at Ri and Rj, is given in atomic units by 1/|Rij|, where Rij=RiRj. The dipole–dipole interaction between point dipoles at Ri and Rj is given by a tensor

Tiμjν=μν1|Rij|=δμν|Rij|33Rij,μRij,ν|Rij|5.
(8)

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

Etot(1)(ω)Etot(2)(ω)Etot(N)(ω)E̲tot(ω)=Eext(1)(ω)Eext(2)(ω)Eext(N)(ω)E̲ext(ω)0T12T1NT210T2NTN1TN20T̲d(1)(ω)d(2)(ω)d(N)(ω)d̲(ω).
(9)

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

d̲(ω)=I̲+α̲0(ω)T̲1α̲0(ω)E̲ext(ω),
(10)

where the reducible unit-wise polarizability tensor is identified,

α̲(ω)=I̲+α̲0(ω)T̲1α̲0(ω)
(11a)
=α̲0(ω)1+T̲1.
(11b)

Each unit contributes to the total dipole moment of the coupled system. Assuming a uniform external electric field throughout the coupled system [Eext(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,

αμν(ω)=iNjN[α̲]iμjν(ω).
(12)

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

Sμ(ω)=2ωπI[αμμ(ω)].
(13)

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 functional43 in the adiabatic limit and RT-TDDFT via the δ-kick technique33 (Sec. II A), as implemented using linear combination of atomic orbital (LCAO) basis sets34 in the GPAW code.44 The projector augmented-wave45 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 functional47,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 model51 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)].

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 benzene data in (b) were obtained from GPAW (RT-TDDFT) and NWChem (Casida) calculations using different XC functionals, as indicated in the legend.

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 benzene data in (b) were obtained from GPAW (RT-TDDFT) and NWChem (Casida) calculations using different XC functionals, as indicated in the legend.

Close modal

Using the dynamic polarizabilities of the isolated units in Eq. (12), we first inspect the optical spectrum of a system comprising an Al201 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.

FIG. 2.

Distance dependence of the optical response from dipolar coupling and full TDDFT calculations. (a) Optical spectra for a system comprising a Al201 NP and two benzene molecules as the NP–molecule distance is varied from RT-TDDFT calculations for the full system (full TDDFT) and from Eq. (12) (dipolar coupling). (b) Spectrum of Al201 +2 benzene system as a practically continuous function of the NP–molecule distance obtained using the dipolar coupling approximation. (c) Coupling strength extracted from Bayesian fits of the spectra in (a) and (b) to a coupled oscillator model. The extracted values are shown as round markers for the full system (full TDDFT) and as a solid line for the dipolar coupling calculations.

FIG. 2.

Distance dependence of the optical response from dipolar coupling and full TDDFT calculations. (a) Optical spectra for a system comprising a Al201 NP and two benzene molecules as the NP–molecule distance is varied from RT-TDDFT calculations for the full system (full TDDFT) and from Eq. (12) (dipolar coupling). (b) Spectrum of Al201 +2 benzene system as a practically continuous function of the NP–molecule distance obtained using the dipolar coupling approximation. (c) Coupling strength extracted from Bayesian fits of the spectra in (a) and (b) to a coupled oscillator model. The extracted values are shown as round markers for the full system (full TDDFT) and as a solid line for the dipolar coupling calculations.

Close modal

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),3 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.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).

FIG. 3.

Dependence of the optical response on the number of molecules from dipolar coupling and full TDDFT calculations. Optical spectra of a system comprising (a) a Al201 NP, (b) a Al586 NP, or (c) a Al1289 NP and a variable number of benzene molecules from RT-TDDFT calculations for the full system (full TDDFT) and from Eq. (12) (dipolar coupling). The solid gray lines show the spectrum of the isolated NPs for comparison.

FIG. 3.

Dependence of the optical response on the number of molecules from dipolar coupling and full TDDFT calculations. Optical spectra of a system comprising (a) a Al201 NP, (b) a Al586 NP, or (c) a Al1289 NP and a variable number of benzene molecules from RT-TDDFT calculations for the full system (full TDDFT) and from Eq. (12) (dipolar coupling). The solid gray lines show the spectrum of the isolated NPs for comparison.

Close modal

For the Al201 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.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 Al201 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 Al586 and Al1289 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.52 

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 Al201 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.53 By matching the diameter of the Mie sphere to the distance between opposing {100} facets of the Al201 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ć54 is considerably larger than the deviation between the present calculation and the spectrum based on the data from McPeak et al.

FIG. 4.

Mixing response calculations from different sources. 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 using shaded lines.

FIG. 4.

Mixing response calculations from different sources. 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 using shaded lines.

Close modal

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 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 medium55 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 package44,57 with LCAO basis sets58 and the LCAO-RT-TDDFT implementation34 was used for the RT-TDDFT calculations. The NWChem suite46 was used for Casida LR-TDDFT calculations. The ASE library59 was used for constructing and manipulating atomic structures. The NumPy,60 SciPy,61 and Matplotlib62 Python packages and the VMD software63,64 were used for processing and plotting data. The emcee65 Python package was used to fit spectra to the coupled oscillator model.

1.
T. W.
Ebbesen
, “
Hybrid light–matter states in a molecular and material science perspective
,”
Acc. Chem. Res.
49
,
2403
2412
(
2016
).
2.
P.
Törmä
and
W. L.
Barnes
, “
Strong coupling between surface plasmon polaritons and emitters: A review
,”
Rep. Prog. Phys.
78
,
013901
(
2015
).
3.
D. G.
Baranov
,
M.
Wersäll
,
J.
Cuadra
,
T. J.
Antosiewicz
, and
T.
Shegai
, “
Novel nanostructures and materials for strong light–matter interactions
,”
ACS Photonics
5
,
24
42
(
2018
).
4.
A.
Frisk Kockum
,
A.
Miranowicz
,
S.
De Liberato
,
S.
Savasta
, and
F.
Nori
, “
Ultrastrong coupling between light and matter
,”
Nat. Rev. Phys.
1
,
19
40
(
2019
).
5.
T.
Yoshie
,
A.
Scherer
,
J.
Hendrickson
,
G.
Khitrova
,
H. M.
Gibbs
,
G.
Rupper
,
C.
Ell
,
O. B.
Shchekin
, and
D. G.
Deppe
, “
Vacuum Rabi splitting with a single quantum dot in a photonic crystal nanocavity
,”
Nature
432
,
200
203
(
2004
).
6.
K. M.
Birnbaum
,
A.
Boca
,
R.
Miller
,
A. D.
Boozer
,
T. E.
Northup
, and
H. J.
Kimble
, “
Photon blockade in an optical cavity with one trapped atom
,”
Nature
436
,
87
90
(
2005
).
7.
J.
Kasprzak
,
M.
Richard
,
S.
Kundermann
,
A.
Baas
,
P.
Jeambrun
,
J. M. J.
Keeling
,
F. M.
Marchetti
,
M. H.
Szymańska
,
R.
André
,
J. L.
Staehli
,
V.
Savona
,
P. B.
Littlewood
,
B.
Deveaud
, and
L. S.
Dang
, “
Bose–Einstein condensation of exciton polaritons
,”
Nature
443
,
409
414
(
2006
).
8.
J. A.
Hutchison
,
T.
Schwartz
,
C.
Genet
,
E.
Devaux
, and
T. W.
Ebbesen
, “
Modifying chemical landscapes by coupling to vacuum fields
,”
Angew. Chem., Int. Ed.
51
,
1592
1596
(
2012
).
9.
F.
Herrera
and
F. C.
Spano
, “
Cavity-controlled chemistry in molecular ensembles
,”
Phys. Rev. Lett.
116
,
238301
(
2016
).
10.
J.
Galego
,
F. J.
Garcia-Vidal
, and
J.
Feist
, “
Suppressing photochemical reactions with quantized light fields
,”
Nat. Commun.
7
,
13841
(
2016
).
11.
L. A.
Martínez-Martínez
,
R. F.
Ribeiro
,
J.
Campos-González-Angulo
, and
J.
Yuen-Zhou
, “
Can ultrastrong coupling change ground-state chemical reactions?
,”
ACS Photonics
5
,
167
176
(
2018
).
12.
B.
Munkhbat
,
M.
Wersäll
,
D. G.
Baranov
,
T. J.
Antosiewicz
, and
T.
Shegai
, “
Suppression of photo-oxidation of organic chromophores by strong coupling to plasmonic nanoantennas
,”
Sci. Adv.
4
,
eaas9552
(
2018
).
13.
J.
Flick
and
P.
Narang
, “
Ab initio polaritonic potential-energy surfaces for excited-state nanophotonics and polaritonic chemistry
,”
J. Chem. Phys.
153
,
094116
(
2020
).
14.
J.
Fregoni
,
G.
Granucci
,
M.
Persico
, and
S.
Corni
, “
Strong coupling with light enhances the photoisomerization quantum yield of azobenzene
,”
Chem
6
,
250
265
(
2020
).
15.
A.
Thomas
,
J.
George
,
A.
Shalabney
,
M.
Dryzhakov
,
S. J.
Varma
,
J.
Moran
,
T.
Chervy
,
X.
Zhong
,
E.
Devaux
,
C.
Genet
,
J. A.
Hutchison
, and
T. W.
Ebbesen
, “
Ground-state chemical reactivity under vibrational coupling to the vacuum electromagnetic field
,”
Angew. Chem., Int. Ed.
55
,
11462
11466
(
2016
).
16.
A.
Thomas
,
L.
Lethuillier-Karl
,
K.
Nagarajan
,
R. M. A.
Vergauwe
,
J.
George
,
T.
Chervy
,
A.
Shalabney
,
E.
Devaux
,
C.
Genet
,
J.
Moran
, and
T. W.
Ebbesen
, “
Tilting a ground-state reactivity landscape by vibrational strong coupling
,”
Science
363
,
615
619
(
2019
).
17.
K.
Hirai
,
R.
Takeda
,
J. A.
Hutchison
, and
H.
Uji-i
, “
Modulation of Prins cyclization by vibrational strong coupling
,”
Angew. Chem., Int. Ed.
59
,
5332
5335
(
2020
).
18.
J.
Feist
and
F. J.
Garcia-Vidal
, “
Extraordinary exciton conductance induced by strong coupling
,”
Phys. Rev. Lett.
114
,
196402
(
2015
).
19.
X.
Zhong
,
T.
Chervy
,
S.
Wang
,
J.
George
,
A.
Thomas
,
J. A.
Hutchison
,
E.
Devaux
,
C.
Genet
, and
T. W.
Ebbesen
, “
Non-radiative energy transfer mediated by hybrid light-matter states
,”
Angew. Chem., Int. Ed.
55
,
6202
6206
(
2016
).
20.
E. T.
Jaynes
and
F. W.
Cummings
, “
Comparison of quantum and semiclassical radiation theories with application to the beam maser
,”
Proc. IEEE
51
,
89
109
(
1963
).
21.
B. M.
Garraway
, “
The Dicke model in quantum optics: Dicke model revisited
,”
Philos. Trans. R. Soc., A
369
,
1137
1155
(
2011
).
22.
J.
Galego
,
F. J.
Garcia-Vidal
, and
J.
Feist
, “
Cavity-induced modifications of molecular structure in the strong-coupling regime
,”
Phys. Rev. X
5
,
041022
(
2015
).
23.
M.
Ruggenthaler
,
J.
Flick
,
C.
Pellegrini
,
H.
Appel
,
I. V.
Tokatly
, and
A.
Rubio
, “
Quantum-electrodynamical density-functional theory: Bridging quantum optics and electronic-structure theory
,”
Phys. Rev. A
90
,
012508
(
2014
).
24.
J.
Flick
,
M.
Ruggenthaler
,
H.
Appel
, and
A.
Rubio
, “
Atoms and molecules in cavities, from weak to strong coupling in quantum-electrodynamics (QED) chemistry
,”
Proc. Natl. Acad. Sci. U. S. A.
114
,
3026
3034
(
2017
).
25.
H. L.
Luk
,
J.
Feist
,
J. J.
Toppari
, and
G.
Groenhof
, “
Multiscale molecular dynamics simulations of polaritonic chemistry
,”
J. Chem. Theory Comput.
13
,
4324
4335
(
2017
).
26.
T.
Neuman
,
R.
Esteban
,
D.
Casanova
,
F. J.
García-Vidal
, and
J.
Aizpurua
, “
Coupling of molecular emitters and plasmonic cavities beyond the point-dipole approximation
,”
Nano Lett.
18
,
2358
2364
(
2018
).
27.
T. P.
Rossi
,
T.
Shegai
,
P.
Erhart
, and
T. J.
Antosiewicz
, “
Strong plasmon-molecule coupling at the nanoscale revealed by first-principles modeling
,”
Nat. Commun.
10
,
3336
(
2019
).
28.
J.
Flick
,
D. M.
Welakuh
,
M.
Ruggenthaler
,
H.
Appel
, and
A.
Rubio
, “
Light–matter response in nonrelativistic quantum electrodynamics
,”
ACS Photonics
6
,
2757
2778
(
2019
).
29.
E.
Orgiu
,
J.
George
,
J. A.
Hutchison
,
E.
Devaux
,
J. F.
Dayen
,
B.
Doudin
,
F.
Stellacci
,
C.
Genet
,
J.
Schachenmayer
,
C.
Genes
,
G.
Pupillo
,
P.
Samorì
, and
T. W.
Ebbesen
, “
Conductivity in organic semiconductors hybridized with the vacuum field
,”
Nat. Mater.
14
,
1123
1129
(
2015
).
30.
P.
Hohenberg
and
W.
Kohn
, “
Inhomogeneous electron gas
,”
Phys. Rev.
136
,
B864
B871
(
1964
).
31.
W.
Kohn
and
L. J.
Sham
, “
Self-consistent equations including exchange and correlation effects
,”
Phys. Rev.
140
,
A1133
A1138
(
1965
).
32.
E.
Runge
and
E. K. U.
Gross
, “
Density-functional theory for time-dependent systems
,”
Phys. Rev. Lett.
52
,
997
1000
(
1984
).
33.
K.
Yabana
and
G. F.
Bertsch
, “
Time-dependent local-density approximation in real time
,”
Phys. Rev. B
54
,
4484
4487
(
1996
).
34.
M.
Kuisma
,
A.
Sakko
,
T. P.
Rossi
,
A. H.
Larsen
,
J.
Enkovaara
,
L.
Lehtovaara
, and
T. T.
Rantala
, “
Localized surface plasmon resonance in silver nanoparticles: Atomistic first-principles time-dependent density-functional theory calculations
,”
Phys. Rev. B
91
,
115431
(
2015
).
35.
M. E.
Casida
, “
Time-dependent density functional response theory for molecules
,” in
Recent Advances in Density Functional Methods, Part I
, edited by
D. P.
Chong
(
World Scientific
,
Singapore
,
1995
), p.
155
.
36.
M. E.
Casida
, “
Time-dependent density-functional theory for molecules and molecular solids
,”
J. Mol. Struct.: THEOCHEM
914
,
3
18
(
2009
).
37.
G.
Mie
, “
Beiträge zur Optik trüber Medien, speziell Kolloidaler Metallösungen
,”
Ann. Phys.
330
,
377
445
(
1908
).
38.
A.
Sihvola
, “
Dielectric polarization and particle shape effects
,”
J. Nanomat.
2007
,
1
9
.
39.
H.
DeVoe
, “
Optical properties of molecular aggregates. I. Classical model of electronic absorption and refraction
,”
J. Chem. Phys.
41
,
393
400
(
1964
).
40.
H.
Fidder
,
J.
Knoester
, and
D. A.
Wiersma
, “
Optical properties of disordered molecular aggregates: A numerical study
,”
J. Chem. Phys.
95
,
7880
7890
(
1991
).
41.
B.
Auguié
,
B. L.
Darby
, and
E. C.
Le Ru
, “
Electromagnetic interactions of dye molecules surrounding a nanosphere
,”
Nanoscale
11
,
12177
12187
(
2019
).
42.
T. P.
Rossi
,
T.
Shegai
,
P.
Erhart
, and
T. J.
Antosiewicz
(
2019
). “
Data and code for strong plasmon-molecule coupling at the nanoscale revealed by first-principles modeling
,” Zenodo.
43.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
, “
Generalized gradient approximation made simple
,”
Phys. Rev. Lett.
77
,
3865
3868
(
1996
).
44.
J.
Enkovaara
,
C.
Rostgaard
,
J. J.
Mortensen
,
J.
Chen
,
M.
Dułak
,
L.
Ferrighi
,
J.
Gavnholt
,
C.
Glinsvad
,
V.
Haikola
,
H. A.
Hansen
,
H. H.
Kristoffersen
,
M.
Kuisma
,
A. H.
Larsen
,
L.
Lehtovaara
,
M.
Ljungberg
,
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.
Walter
,
B.
Hammer
,
H.
Häkkinen
,
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
, “
Electronic structure calculations with GPAW: A real-space implementation of the projector augmented-wave method
,”
J. Phys.: Condens. Matter
22
,
253202
(
2010
).
45.
P. E.
Blöchl
, “
Projector augmented-wave method
,”
Phys. Rev. B
50
,
17953
17979
(
1994
).
46.
M.
Valiev
,
E. J.
Bylaska
,
N.
Govind
,
K.
Kowalski
,
T. P.
Straatsma
,
H. J. J.
Van Dam
,
D.
Wang
,
J.
Nieplocha
,
E.
Apra
,
T. L.
Windus
, and
W. A.
de Jong
, “
NWChem: A comprehensive and scalable open-source solution for large scale molecular simulations
,”
Comput. Phys. Commun.
181
,
1477
1489
(
2010
).
47.
C.
Lee
,
W.
Yang
, and
R. G.
Parr
, “
Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density
,”
Phys. Rev. B
37
,
785
789
(
1988
).
48.
A. D.
Becke
, “
Density-functional thermochemistry. III. The role of exact exchange
,”
J. Chem. Phys.
98
,
5648
5652
(
1993
).
49.
A. D.
McLean
and
G. S.
Chandler
, “
Contracted Gaussian basis sets for molecular calculations. I. Second row atoms, z = 11−18
,”
J. Chem. Phys.
72
,
5639
5648
(
1980
).
50.
R.
Krishnan
,
J. S.
Binkley
,
R.
Seeger
, and
J. A.
Pople
, “
Self-consistent molecular orbital methods. XX. A basis set for correlated wave functions
,”
J. Chem. Phys.
72
,
650
654
(
1980
).
51.
X.
Wu
,
S. K.
Gray
, and
M.
Pelton
, “
Quantum-dot-induced transparency in a nanoscale plasmonic resonator
,”
Opt. Express
18
,
23633
23645
(
2010
).
52.
Z.-J.
Yang
,
T. J.
Antosiewicz
, and
T.
Shegai
, “
Role of material loss and mode volume of plasmonic nanocavities for strong plasmon-exciton interactions
,”
Opt. Express
24
,
20373
20381
(
2016
).
53.
K. M.
McPeak
,
S. V.
Jayanti
,
S. J. P.
Kress
,
S.
Meyer
,
S.
Iotti
,
A.
Rossinelli
, and
D. J.
Norris
, “
Plasmonic films can easily be better: Rules and recipes
,”
ACS Photonics
2
,
326
333
(
2015
).
54.
A. D.
Rakić
, “
Algorithm for the determination of intrinsic optical constants of metal films: Application to aluminum
,”
Appl. Opt.
34
,
4755
4767
(
1995
).
55.
V.
Barone
,
M.
Cossi
, and
J.
Tomasi
, “
A new definition of cavities for the computation of solvation free energies by the polarizable continuum model
,”
J. Chem. Phys.
107
,
3210
3221
(
1997
).
56.
J.
Fojt
,
T. P.
Rossi
,
T. J.
Antosiewicz
,
M.
Kuisma
, and
P.
Erhart
(
2021
). “
Data and code for “dipolar coupling of nanoparticle-molecule assemblies: An efficient approach for studying strong coupling
,” Zenodo.
57.
J. J.
Mortensen
,
L. B.
Hansen
, and
K. W.
Jacobsen
, “
Real-space grid implementation of the projector augmented wave method
,”
Phys. Rev. B
71
,
035109
(
2005
).
58.
A. H.
Larsen
,
M.
Vanin
,
J. J.
Mortensen
,
K. S.
Thygesen
, and
K. W.
Jacobsen
, “
Localized atomic basis set in the projector augmented wave method
,”
Phys. Rev. B
80
,
195112
(
2009
).
59.
A. H.
Larsen
,
J. J.
Mortensen
,
J.
Blomqvist
,
I. E.
Castelli
,
R.
Christensen
,
M.
Dulak
,
J.
Friis
,
M. N.
Groves
,
B.
Hammer
,
C.
Hargus
,
E. D.
Hermes
,
P. C.
Jennings
,
P. B.
Jensen
,
J.
Kermode
,
J. R.
Kitchin
,
E. L.
Kolsbjerg
,
J.
Kubal
,
K.
Kaasbjerg
,
S.
Lysgaard
,
J. B.
Maronsson
,
T.
Maxson
,
T.
Olsen
,
L.
Pastewka
,
A.
Peterson
,
C.
Rostgaard
,
J. S. O.
Schütt
,
M.
Strange
,
K. S.
Thygesen
,
T.
Vegge
,
L.
Vilhelmsen
,
M.
Walter
,
Z.
Zeng
, and
K. W.
Jacobsen
, “
The atomic simulation environment–A Python library for working with atoms
,”
J. Phys.: Condens. Matter
29
,
273002
(
2017
).
60.
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.
Haldane
,
J. F.
del Río
,
M.
Wiebe
,
P.
Peterson
,
P.
Gérard-Marchant
,
K.
Sheppard
,
T.
Reddy
,
W.
Weckesser
,
H.
Abbasi
,
C.
Gohlke
, and
T. E.
Oliphant
, “
Array programming with NumPy
,”
Nature
585
,
357
362
(
2020
).
61.
P.
Virtanen
,
R.
Gommers
,
R.
Gommers
,
T. E.
Oliphant
,
M.
Haberland
,
T.
Reddy
,
D.
Cournapeau
,
E.
Burovski
,
P.
Peterson
,
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
,
İ.
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.0 Contributors
, “
SciPy 1.0: Fundamental algorithms for scientific computing in Python
,”
Nat. Methods
17
,
261
272
(
2020
).
62.
J. D.
Hunter
, “
Matplotlib: A 2D graphics environment
,”
Comput. Sci. Eng.
9
,
90
95
(
2007
).
63.
W.
Humphrey
,
A.
Dalke
, and
K.
Schulten
, “
VMD – Visual molecular dynamics
,”
J. Mol. Graph.
14
,
33
(
1996
).
64.
J.
Stone
, “
An efficient library for parallel ray tracing and animation
,” M.S. thesis,
Computer Science Department, University of Missouri-Rolla
,
1998
.
65.
D.
Foreman-Mackey
,
D. W.
Hogg
,
D.
Lang
, and
J.
Goodman
, “
emcee: The MCMC hammer
,”
Publ. Astron. Soc. Pac.
125
,
306
(
2013
).

Supplementary Material