# Ultrastrong Coupling of a Single Molecule to a Plasmonic Nanocavity: A First-Principles Study

Mikael Kuisma, Benjamin Rousseaux, Krzysztof Czajkowski, Tuomas Rossi, Timur Shegai, Paul Erhart, Tomasz Antosiewicz

► To cite this version:

Mikael Kuisma, Benjamin Rousseaux, Krzysztof Czajkowski, Tuomas Rossi, Timur Shegai, et al. Ultrastrong Coupling of a Single Molecule to a Plasmonic Nanocavity: A First-Principles Study. ACS photonics, 2022, 9, pp.1065 - 1077. 10.1021/acsphotonics.2c00066. http://dx.doi.org/10.1021/acsphotonics.2c00066

HAL Id: hal-03674866

https://u-paris.hal.science/hal-03674866

Submitted on 21 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Ultrastrong Coupling of a Single Molecule to a Plasmonic Nanocavity: A First-Principles Study

Mikael Kuisma, Benjamin Rousseaux, Krzysztof M. Czajkowski, Tuomas P. Rossi, Timur Shegai, Paul Erhart, and Tomasz J. Antosiewicz*

<figure><img src="image_2.png" /></figure>

Cite This: ACS Photonics 2022, 9, 1065–1077

<figure><img src="image_3.png" /></figure>

Read Online

ACCESS

<figure><img src="image_4.png" /></figure>

Metrics & More

<figure><img src="image_5.png" /></figure>

Article Recommendations

<figure><img src="image_6.png" /></figure>

Supporting Information

**ABSTRACT:** Ultrastrong coupling (USC) is a distinct regime of light-matter interaction in which the coupling strength is comparable to the resonance energy of the cavity or emitter. In the USC regime, common approximations to quantum optical Hamiltonians, such as the rotating wave approximation, break down as the ground state of the coupled system gains photonic character due to admixing of vacuum states with higher excited states, leading to ground-state energy changes. USC is usually achieved by collective coherent coupling of many quantum emitters to a single mode cavity, whereas USC with a single molecule remains challenging. Here, we show by time-dependent density functional theory (TDDFT) calculations that a single organic molecule can reach USC with a plasmonic dimer,

consisting of a few hundred atoms. In this context, we discuss the capacity of TDDFT to represent strong coupling and its connection to the quantum optical Hamiltonian. We find that USC leads to appreciable ground-state energy modifications accounting for a non-negligible part of the total interaction energy, comparable to $k_B T$ at room temperature.

**KEYWORDS:** strong coupling, time-dependent density functional theory, plasmonics, nanophotonics, excitons

<figure><img src="image_7.png" /></figure>

Strong light-matter coupling is a regime characteristic of interacting systems which can no longer be treated perturbatively. In such a hybridized state the light and matter subsystems are described as dressed polaritonic states with modified behavior that extends beyond optical properties.¹ These dressed states form when the rate of coupling exceeds the individual damping rates and excitations persist long enough to allow coherent energy exchange between the subsystems. Research in this field² has already led to exciting observations of strong-coupling-induced modifications of exciton transport,³ polaron photoconductivity,⁴ photochemical rates,⁵ and ground-state reactivity⁶ as well as single-photon nonlinearities.⁷ For many of these strong-coupling induced phenomena, it is sufficient that the coupling strength is only slightly larger than the damping rates. Commonly, strong coupling can therefore be reached already at coupling rates that are much smaller than the resonance energies of the constituent components. By contrast, ultrastrong coupling (USC) requires the coupling strength to be comparable to the resonance or transition energies in the case of cavity and matter, respectively.

In the dipole approximation, the magnitude of the coupling strength g of cavity-light-matter interaction is determined by the transition dipole moment of a single element of matter $\mu_1 = |\mu_1|$, the number of these elements N, and the vacuum field of the cavity $E_{vac} = |E_{vac}|$ according to¹

$$g = \sqrt{N}\mu_1 \cdot E_{vac} \quad \text{with} \quad E_{vac} = \sqrt{\frac{\hbar\omega}{2\epsilon_0\epsilon V}} \quad (1)$$

where $E_{vac}$ is evaluated at the quasinormal mode (QNM) frequency $\omega$ corresponding to the (complex) mode volume $V$.⁸⁻¹⁰ The strong coupling formalism expressed in the Jaynes–Cummings (JC) Hamiltonian is simple to utilize for interpreting experimental or numerical results when assuming a single cavity mode. However, when two (or more) QNMs are present and may potentially interfere, such as in systems incorporating lossy plasmonic nanocavities, the JC approach is insufficient.¹¹ Specifically, it becomes necessary to account for dissipative coupling between the QNMs, resulting in a reduction of the cavity-emitter coupling rate,¹² and may result in Fano-like interference.¹¹ In a single QNM system, however, such as in simple plasmonic systems, in which subwavelength-localized longitudinal fields dominate, eq 1 will remain as the characteristic relation of these systems, such that the effective volume V is

Received: January 12, 2022

Published: March 2, 2022

<figure><img src="image_8.png" /></figure>
<figure><img src="image_9.png" /><figcaption>Figure 1. Schematic illustration and photoabsorption spectra of the studied nanoparticle-dimer systems spanning a broad spectral range from visible to UV. From left to right: sodium dimers with 30 Å gap (resonance at 3.0 eV) with nonacene (electronic transition of μ = 15.3 D at 2.9 eV), magnesium dimers with 18 Å gap (4.6 eV) with tetracene (11.1 D, 4.4 eV), and aluminum dimers with 10 Å gap (7.2 eV) with benzene (4.4 D, 7.2 eV).</figcaption></figure>

the volume of a transverse cavity with equal coupling strength,¹³,¹⁴ as long as the imaginary part of V remains negligible compared to its real part.¹⁰ When g is at least on the order of 10% of the molecular transition energy ω<sub>ex</sub> (ζ = g/ω<sub>ex</sub> ≥ 0.1), the system is said to be ultrastrongly coupled. An important motivation for reaching the USC regime is a non-negligible ground state (GS) energy modification

$$ \Delta E = \frac{g^2}{2\omega_{\text{ex}}} \qquad (2) $$

predicted by cavity quantum electrodynamics (cQED) assuming a single transverse optical mode.¹⁵,¹⁶ In the longitudinal case, the GS energy modification due to zero-point energy shifts¹⁷ is attractive and, in the perturbative limit, is known as van der Waals (vdW) attraction. This lowering of the GS energy is in contrast to the repulsive GS energy modification that is obtained in the naive 2-excitation transverse model, as detailed in the Discussion section in which we summarize this dichotomy. Under favorable conditions, ΔE could be on the order of k<sub>B</sub>T at room temperature, a value which would have significant impact on polaritonic chemistry.⁶,¹⁸ By looking at eq 1 and eq 2 one can observe that reaching USC and thus large ΔE becomes easier at lower frequencies. Specifically, ζ ∝ μ₁√N/ω<sub>ex</sub>V, while the expected USC-induced GS energy modification is to the first approximation independent of the transition energy, |Δ<sub>e</sub>| ∝ Nμ<sub>1</sub><sup>2</sup>/V.

In the pursuit of USC with optical cavities at room temperature, recent work has shown that very large coupling strengths in excess of half of the resonance/transition energy can be obtained even with only 4% of a Fabry-Pérot cavity occupied.¹⁹ This is possible owing to the very large oscillator strengths of plasmonic nanorods, and the approach can be engineered further to reach even larger coupling strengths.²⁰ However, molecules have typically much smaller transition dipole moments than the collective oscillations associated with the conduction electrons in metallic nanoparticles. In the case of molecules, one therefore requires much larger numbers (densities) to reach comparable coupling strengths, typically by filling the entirety of a cavity. Coupling a cavity mode with N single-transition molecules leads to N + 1 total states, including two (bright) polaritonic states and N − 1 purely molecular dark states. The dark states play, however, a crucial role in interactions as they entropically undermine strong coupling.

For example, the upper polariton quickly decays into localized dark states.²¹,²² By contrast, in the case of a single molecule (N = 1) the number of dark states is zero. This provides a motivation for limiting the number of molecules per cavity.

Here, we explore an alternative approach: Instead of maximizing the coupling strength by using extremely large values of μ₁ and/or N, we reach USC by increasing the vacuum electric field E<sub>vac</sub>. A small mode volume V on the order of 100 nm³ may be obtained in dimers²³ or particle-on-mirror systems²⁴ despite the large size of the employed optical antennas. Even smaller cavities may be reached by utilizing atomic-scale cavities composed of metal clusters with a few hundred atoms, either alone²⁵ or arranged in nanoscale dimers,²⁶ whose modes can be tuned by the gap size and/or shape.²⁷,²⁸ Building on a recent time-dependent density functional theory (TDDFT) study of strong coupling in a benzene-Al<sub>201</sub> system,²⁵ we utilize a cavity made up of two identical plasmonic nanoparticles, resulting in a highly enhanced electric field in a small gap. We then couple the nanoparticle dimer with several different π-conjugated molecules to quantify the dependence of the coupling strength and ΔE on the material parameters. We emphasize that the TDDFT approach used here takes into account plasmon decay due to Landau damping and the corresponding dissipation (lower Q-values), while electron-electron and electron-phonon scattering processes are taken into account by empirical broadening. This enhances the reliability of the present approach compared to models that rely on two-level systems and/or fully empirical broadening.

Since such an ab initio approach allows us to quantitatively investigate longitudinal plasmonic cavities, it is insightful to compare and highlight their differences to transverse cavities. The necessity of the quadratic diamagnetic terms in the transverse Hamiltonian has been discussed in the literature extensively.²⁹⁻³¹ For example, they protect nongapped systems from infrared divergences³² and, in the case of a cavity, offer protection from superradiant phase transitions. Recently Schäfer et al. suggested that the same effects should be present in plasmonic cavities.³¹ Here, we therefore discuss the manifestation of these effects based on quantitative ab initio calculations as well as theoretical analysis. We find that they, however, do not emerge upon introduction of coupling per se, but rather result from the intraparticle effect of plasmonic self-polarization, which also protects from superradiant transitions in interparticle coupling.
<figure><img src="image_10.png" /><figcaption>**Figure 2.** (a) Photoabsorption spectra of Al<sub>147</sub> dimers with a 10 Å gap coupled strongly to benzene (solid line) and reference spectra (dotted). (b) Induced charge density of benzene, Al<sub>147</sub> dimer, and the lower/upper polaritons showing the in-phase and out-of-phase combinations. The 3D maps are plotted at the relevant resonant frequencies which are also labeled in the subsequent panels. (c–f) Transition contribution maps (TCMs) of the individual and coupled systems, where red and blue indicate positive and negative contributions, respectively. The TCMs of both polaritons consist of that of benzene and the dimer. At the LP the benzene transition (circled) is in phase with the Al plasmon, while at the UP it is out-of-phase, consistent with the in-phase and out-of-phase nature of LP and UP, respectively.</figcaption></figure>

## RESULTS

We model cavities made of sodium (Na), magnesium (Mg), and aluminum (Al) nanoparticle dimers, which exhibit plasmon resonances at approximately 3, 5, and 7 eV, respectively, covering a broad spectral range. These cavities are paired with aromatic hydrocarbon molecules that exhibit highest occupied molecular orbital (HOMO)-lowest unoccupied molecular orbital (LUMO) electronic excitations close to the nanoparticle-dimer plasmon resonances (Figure 1). Their resonance with the respective plasmons ensures adequate spectral overlap even if molecules fill the gap and cause a red shift of the plasmon,³³ which occurs in parallel to USC. The molecular excitation is tuned to each system by extending the chain from a single ring in benzene to ten rings in decacene. This elongation results in an increase of the transition dipole moment, but simultaneously forces a wider dimer gap to accommodate a longer molecule. Hence, the resulting coupled spectra are determined by a number of factors, some of which are beyond explicit control, as they are dictated by tuning the cavity resonance frequency. All the modeled combinations of dimers, molecules, and spatial arrangements are tabulated in *Supporting Information* (SI) Table S1, while representative dimer-molecule geometries are displayed in Figure 1 and Figure S1.

Most of the coupled systems considered here feature a single molecule, whose size sets the minimum gap size of the nanoparticle dimer. For example, for a sodium dimer coupled to nonacene the minimum gap is 28 Å. As we will show in the following, these single-molecule-nanoparticle-dimer systems are capable of reaching USC. To gain additional insight into the

intricacies of nanoscale polaritons, we also investigate the effects of smaller gap sizes by utilizing closely spaced parallel dimers of shorter molecules. These artificial molecular dimers are designed to have their joint molecular transition at a similar energy as that of a longer molecule. For example, in sodium dimers this condition is fulfilled by two pentacene molecules placed 2 Å apart owing to their mutual interaction. These structures allow us to investigate the limits of USC in nanoscale systems. Furthermore, for selected dimers such as Al<sub>147</sub> and Mg<sub>309</sub>, we also analyze the impact of nanoparticle orientation on the coupling strength by considering facet-to-facet, edge-to-edge, and corner-to-corner alignments, owing to the fact that such nanoscale details play a crucial role in determining the modal structure of atomic-scale structures.<sup>27,28,34</sup>

**Aluminum.** Figure 2a shows typical TDDFT-calculated photoabsorption spectra, specifically of Al<sub>147</sub> dimers coupled to a single benzene molecule. The Rabi splitting already for a single benzene is approximately 1 eV and both polaritons become more pronounced with increasing number of molecules. For dimers composed of the larger Al<sub>309</sub> cluster (Figure S2), despite the same gap size of 10 Å, the splitting of the two peaks is weaker. In fact, for the latter case the spectrum barely exhibits a dip between upper polariton (UP) and lower polariton (LP) with only one molecule involved in the coupling, although for increasing *N* much clearer polaritons develop.

To show that the two observed peaks are indeed the strongly coupled bonding and antibonding of, respectively, the LP and UP, we plot the induced electronic charge densities at the resonances/transitions of benzene, the nanoparticle dimer, and
<figure><img src="image_11.png" /><figcaption>**Figure 3.** Photoabsorption spectra of strongly coupled molecule–nanoparticle-dimer systems (solid lines, gray dotted lines mark reference spectra): (a,b) Na nanodimers coupled to nonacene (reference Na dimer gap 30 Å), (c,d) Mg nanodimers to tetracene (reference Mg dimer gaps 18 Å). (e,f) Induced charge densities at the lower and upper polaritons for Na and Mg dimers coupled to molecules showing the bonding and antibonding character of a strongly coupled system. (g) Fitted coupling strength $g_{\text{fit}}$ for all systems versus the fitted molecular transition energy $\omega_{\text{ex}}$. The results are arranged in three groups: Na dimers (squares), Mg dimers (crosses), and Al dimers (circles). The insets in the legend mark the different cases; for a single shape/color the varying parameter is either the gap size (continuous variation of points) or the number of molecules (scattered points). The dotted line marks $0.1\omega_{\text{ex}}$ the commonly adopted lower limit of USC. The coupling strengths vary significantly with gap size. For all nanoparticle dimers the USC regime can be reached through a suitable combination with molecules.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 4.** (a–c) The electric field enhancement of empty Mg₃₀₉ dimers with 18 Å gap differs qualitatively with the relative orientation of the nanoparticles, including significant quantitative differences in their maximum amplitudes. In contrast, coupling to tetracene, leads to induced fields at (d–f) lower polaritons and (g–i) upper polaritons that are qualitatively similar, which illustrates a significant modification of the cavity by the presence of the molecule that occurs in addition to the USC of the molecular electronic transition to the nanoscale-dimer plasmons.</figcaption></figure>

both polaritons (Figure 2b). The dipolar character of the transition in benzene is evident, as is the bright (dipolar) mode of the dimer.²⁷ At both polaritons the polarization of the dimer is the same. At the LP and UP benzene is polarized in-phase and

out-of-phase, respectively, demonstrating the two orthogonal mixed states of benzene and the dimer.

Further, by analyzing the transition contribution maps (TCMs) (Figure 2c–f), which visualize the Kohn–Sham
(KS) electron-hole transition contributions to photoabsorption.³⁵ The benzene transition at 7.19 eV is doubly degenerate with an energy difference between the occupied and unoccupied states of approximately 6.5 eV, while the Al₁₄₇-dimer plasmon consists of a large number of coherent transitions from just below to just above the Fermi energy with a predominant energy difference of approximately 2 eV. The TCMs of the two polaritons, which are formed by mixing of the molecular exciton and plasmon, show the same characteristic distribution of the KS transitions.²⁵ However, at the LP the benzene transitions are in-phase with those of the plasmon, while at the UP they are out-of-phase, hence screening the plasmon. These observations are qualitatively identical for other studied systems.

**Sodium.** Figure 3a,b presents additional photoabsorption spectra of selected coupled systems (also see Figure S2), clearly showing LP and UP. For the lowest energy case, sodium Na₁₄₇ dimer in Figure 3a, the Rabi splitting is on the order of 0.5 eV for the smallest dimer gap with LP and UP being of equal amplitude, indicating only minor detuning between the molecular and Na₁₄₇ resonances. In Figure 3e we show the corresponding induced charge densities. For the coupled Na₃₀₉ dimer (Figure 3b), the UP is more pronounced than the LP due to larger detuning, although the plasmon peak positions are similar in Na₁₄₇ and Na₃₀₉. The similar Rabi splitting in these systems indicates that the coupling strength is sensitive to the gap size but not to the size of the nanoparticles.

**Magnesium.** In Figure 3c, Mg₂₀₁ dimers are coupled to tetracene with smaller gaps than in the Na dimers, while a plot of an exemplary induced charge density at the LP/UP is shown in Figure 3f. The coupling strength for the smallest gap is 0.6 eV, a significant fraction of the 4.4 eV transition energy. For the larger Mg₃₀₉ dimer with tetracene (Figure S2) the coupling strengths are up to 15% smaller than for Mg₂₀₁ for the same gap size, which is in line with the larger mode volume of Mg₃₀₉. Additionally, the coupling strength should be affected by the relative orientation of the dimer components. For example, a dimer in corner-on-corner orientation localizes the electric field more strongly than the facet-on-facet configuration (Figure 4a,c),²⁷ which suggests that a single-atom protrusion may be considered a picocavity.³⁴,³⁶ Also, at the level of the empty dimers one observes slightly different plasmon energies for larger gaps as well as the appearance of a charge transfer plasmon in addition to the bright dipolar mode.³⁴ Combined, these differences between the differently oriented dimers result in different localization of the modes with variations by up to a factor of 4²⁷ with the induced electric fields shown in Figure 4a-c. Such differences should theoretically result in coupling strengths varying by as much as 50%. This is, however, not observed. Indeed, the Rabi splittings in all three dimer orientations are nearly identical when coupled to tetracene (Figure 3d). The only differences are small variations of the relative amplitudes of the LP and UP, which are caused by the slight differences between the cavity resonance energies and different detuning with respect to the electronic transition of the molecule. The origin of the equal coupling strength in all three cases is made clear via the induced electric fields of LPs and UPs, especially in the gaps. Despite the different orientation of the dimer elements, both polaritons are quantitatively very similar (Figure 4d-i), demonstrating that in such coupled systems the molecule changes the character of the cavity, overriding the unique individual modal distributions when forming hybridized modes. Thus, the picoscale electric field localization offered by single atoms³⁴,³⁶ may disappear with increasing coupling, resulting in a mode volume that is

determined by the molecule(s) rather than the atomic features of the metal antenna(s).²³,³⁷

**Single-Molecule USC.** The photoabsorption spectra are fitted with the velocity-coupled harmonic oscillator model, Methods section and eq 7, to obtain the coupling strengths g<sub>fit</sub> (fitting parameter), which are plotted in Figure 3g as a function of the fitted molecular excitation energy ω<sub>ex</sub> (see Figure S3 for per-molecule coupling strength). Grouped by cavity material from lowest to highest resonance energy, the results show that very strong interaction is possible even with a single molecule occupying the gap. When a single molecule is replaced by a pair of molecules (an artificial molecular dimer; Figure S1d,f) to enable the use of smaller gaps, g becomes even larger. Already the coupling strengths of the single-molecule cases exceed 0.1ω<sub>ex</sub> (dotted line), reaching similar values (ζ = g/ω<sub>ex</sub>) for all three systems in the range 0.11 to 0.13. For two molecules in the gap ζ is even larger, reaching ζ = 0.22 for Na, ζ = 0.18 for Al, and ζ = 0.16 for Mg (see Supplementary Table S2).

## DISCUSSION

**Effective Vacuum Fields.** In the strong coupling regime, and especially in the USC limit, discussing the system in terms of undressed molecules and bare dimers is not appropriate, since the underlying components are thoroughly mixed into the new polaritonic states. This is demonstrated by the qualitative similarity of the field distribution at LP and UP (Figure 4) in different configurations, as local atomic-scale variations in the individual components become much less pronounced in the coupled systems. To gain intuition, it is instructive to refer to the basic characteristics of a strongly coupled system as expressed via eq 1, by calculating the effective vacuum fields and comparing them to the theoretical values (see Methods) of the corresponding empty dimers. The latter requires an evaluation of the volume of the corresponding QNMs.⁸

For our discussion we define the effective vacuum field E<sub>vac</sub><sup>eff</sup> based on the coupling strengths g<sub>fit</sub> obtained from fitting the spectra with the coupled harmonic oscillator model in eq 7, as the magnitude of the electric field necessary to obtain g<sub>fit</sub> when acting onto the transition dipole moment of N molecules coupled to the dimer. It is calculated as

$$E_{\text{vac}}^{\text{eff}} = \frac{g_{\text{fit}}}{\sqrt{N} \mu_1} \qquad (3)$$

based on eq 1. Thus, E<sub>vac</sub><sup>eff</sup> accounts not only for the inhomogeneous field of the plasmonic cavity, but also for other coupling-induced modifications of the interacting dimer-molecule system that are not captured by simple models. The additional effects include contributions from higher energy transitions of the molecule that reshape the cavity mode, the interaction between KS states of the interacting entities or shift of the cavity resonance.

As TDDFT deals directly with electronic properties of matter rather than derived optical quantities such as permittivities, direct use of the existing computational QNM³⁸ formalism is not feasible. Hence, to compute the QNMs as outlined in ref 39 (see Methods for additional details), we approximate the atomic clusters by Drude nanospheres whose permittivity is tailored to match the peak position, width, and amplitude obtained with TDDFT and the radius is determined by the physical size of the cluster (see Table 1 for parameter values). This allows us to subsume all nonlocal and quantum size effects⁴⁰ into an effective permittivity, which is a good approximation for the auxiliary role
<table><caption>**Table 1. Parameters Used in the Classical Electromagnetic Calculation of the Quasinormal Modes of Nanodimers:**

**Nanosphere Radius, Drude Permittivity:** $\epsilon(\omega) = 1 - \omega_p^2/(\omega^2 + i\gamma\omega)$</caption><thead><tr><th>structure</th><th>radius (Å)</th><th>ħω<sub>p</sub> (eV)</th><th>ħγ (eV)</th></tr></thead><tbody><tr><td>Al<sub>147</sub></td><td>7.4</td><td>13.62</td><td>0.75</td></tr><tr><td>Mg<sub>201</sub></td><td>9.8</td><td>5.30</td><td>0.47</td></tr><tr><td>Na<sub>147</sub></td><td>10.1</td><td>8.03</td><td>0.47</td></tr></tbody></table>

of employing the QNMs to verify the effective vacuum fields obtained from TDDFT calculations and $g_{fit}$. The mode volumes, being a position-dependent quantity, are calculated assuming the probing dipole is placed parallel to the dimer axis in the geometrical center of an empty dimer gap. To account for a finite size of the molecules, the dipole is displaced from the center to 1/fourth of the gap size or at the position of the end of a molecule. Next, we use the QNM volumes (Figure S7) to calculate the corresponding QNM vacuum field $E_{vac}^{QNM}$ that acts on a point-like molecule. The imaginary parts $Im\{V\}$ of the QNM volumes are also shown to be at least ten times smaller than the real parts $Re\{V\}$, Figure S8).

The computed effective vacuum field values are plotted in Figure 5a as a function of the fitted electronic transition energy of the molecule. $E_{vac}^{eff}$ approaches $10 \text{ V} \cdot \text{nm}^{-1}$ for Al$_{147}$, $5 \text{ V} \cdot \text{nm}^{-1}$ for Mg$_{201}$, and $2 \text{ V} \cdot \text{nm}^{-1}$ for Na$_{147}$. The ratio of the maximum and minimum effective vacuum fields for the three metals is approximately 20 and is markedly larger than the corresponding ratio for the coupling strengths, which approaches 6. This can be rationalized by the competition between a decreasing mode confinement when going from Al to Mg to Na in order to accommodate increasingly longer molecules with an increasing transition dipole moment. This decrease in mode confinement with increasing gap is clearly visible Figure 5b-d, where we compare the effective and theoretical vacuum fields.

In general, the vacuum fields corresponding to the bright dipolar bonding QNMs, $E_{vac}^{QNM}$, increase for decreasing gaps, as expected with the probing dipole moving closer to the

nanoparticles. They also increase for a fixed gap, when the dipole is brought closer to one of the dimer components. For all three Na, Mg, and Al dimers the vacuum fields are similar for any given gap size and span a range from ca. 0.2 to 7 eV, although the latter's $E_{vac}$ is the largest of these three. This is, however, mainly caused by an increase of the frequency of the QNM, since for a given gap size the QNM volumes of all dimers are similar (Figure S7). Also, due to the very small size of the particles (ca. 1 nm), any higher order modes are much weaker and their mode volumes are a few orders of magnitude larger yielding negligible interaction with the molecules/dipoles.

The effective vacuum fields, which are derived from the fitted coupling strength, are larger than the ones based on the QNMs regardless of the probing position, be it in the middle of the gap or at the position of the end of the molecule. This demonstrates the complex interactions that clearly go beyond a simple two-by-two coupling model, highlighting the fact that the molecule itself reshapes the cavity mode. This further is demonstrated in Figure S4, which shows the TDDFT-calculated field profiles for Mg$_{201}$ dimers with gap sizes of 15 and 30 Å, respectively (Supplementary Note S1). The average energy density of the induced electric field in a volume corresponding to a centrally placed molecule for the 30 Å gap is on the order of 10% of the maximum value (observed near one of the Mg atomic clusters). However, the ratio of $E_{vac}^{eff}/E_{vac}^{QNM}$ is approximately 2 for both one tetracene and two naphthalene molecules. For the smaller gap of 15 Å, the average normalized field intensity is ≈0.4, while in the coupled system the ratio $E_{vac}^{eff}/E_{vac}^{QNM}$ is 2 for tetracene and 3 for two naphthalenes. Clearly, the molecules couple to the cavity mode with greater efficiency than predicted by the vacuum field of the bare cavity and position of the molecule.

The interaction in such nanoscale polaritonic systems extends beyond coupling the cavity mode to the electronic transition, as it also involves a non-negligible modification of the cavity and facilitates reaching USC (Figure 4). Indeed, based on our findings we ascertain that the cavity is modified by electronic states, which are not involved in the molecular HOMO-LUMO transition that couples directly to the plasmon. Molecules not

<figure><img src="image_13.png" /><figcaption>**Figure 5.** (a) Effective vacuum field $E_{vac}^{eff}$ of the dimer cavities interacting with the molecules for Na (squares), Mg (crosses), and Al (circles) for various molecules and gap sizes. A larger $E_{vac}^{eff}$ for the same colored symbol typically corresponds to a smaller gap. For Na the effective vacuum field is small and does not vary by a lot, being a consequence of large gap sizes and small overlap between mode and molecule. This applies also for the systems involving Mg, except for a Mg dimer coupled to a naphthalene dimer, for which $E_{vac}^{eff}$ changes significantly. For Al the vacuum field spans a broad range of values. (b,c,d) Comparison of effective vacuum fields (symbols) of dimer-molecule systems with the QNM vacuum field (lines) of empty dimers versus gap size with the probing dipole placed in the center of the gap (1/2, black) or displaced toward one of the dimer components to 1/4 (light gray) of the gap. The colored lines mark vacuum fields of empty dimers for dipoles placed at an end of the corresponding molecule (whose color in the plot matches the color of the line).</figcaption></figure>
only couple as dipolar resonators but form a kind of dielectric bridge focusing the field of the cavity mode and increasing the effective coupling in the system. This modification of the cavity beyond a two-level model by the *body* of a molecule is reminiscent of the effect of background permittivity of a quantum dot interacting with a plasmonic bow-tie antenna.²³

**USC Contribution to Ground State Energy Shifts.** We now estimate the ground state energy shifts in the USC regime. As in recent quantum electrodynamics DFT (QEDFT) calculations,⁴¹ which require the quadratic diamagnetic or self-polarization terms to be accounted for consistency,³¹ in our case also, analogous terms appear and ensure consistent results as discussed in the previous section. Hence, the computed spectra and underlying coupling strengths are sound, as is the case also in classical calculations,¹⁹ and the results can be interpreted in the framework of quantum mechanical Hamiltonians, namely the Hopfield model.¹⁵ One of the results of the Hopfield Hamiltonian is a predicted modification of a system's GS due to admixing states with higher numbers of excitations, which in the USC regime may constitute a significant fraction of $k_B T$ at room temperature.¹⁹,²⁰ However, due to the small cavity sizes in our study, we use the longitudinal Hamiltonian as derived in **Supplementary Note S2** as a simple model to interpret the results to an adequate approximation. This is justified since the imaginary parts of the mode volumes (responsible for non-Hermiticity and dissipative couplings) were shown to be at least ten times smaller than the respective real parts (Figure S7). The equivalence of this model to TDDFT is discussed in the following section with technical details presented in **Supplementary Note S3**. At zero cavity-exciton detuning the resulting shift of the ground state energy is $\Delta E_L = -g^2/2\omega_{ex}$ (Supplementary Note S2). We emphasize that our calculated longitudinal zero-point energy shifts are concrete random phase approximation (RPA)-correlation energies describing vdW forces in a low coupling limit as discussed in **Supplementary Note S2**.

Using the fitted coupling strengths (Figure 3e), the longitudinal zero-detuning cQED result is used to calculate the expected GS modifications (Figure 6). The calculated $\Delta E_L$

<figure><img src="image_14.png" /><figcaption>**Figure 6.** Predicted ground state energy modification based on cQED estimate in the single molecule USC regime is on the order of $k_B T$ at room temperature for all considered systems. Al dimers yield the largest GS of up to 60 meV, but this is observed only for gaps ≲10 Å, in which the large coupling strength is obtained by a combination of small gap size and molecule-induced modification of the mode volume. In Mg and Na dimers the predicted GS modification is smaller than for Al dimers.</figcaption></figure>

are on the order of a few tens of meV, values which constitute a significant fraction of $k_B T$ at room temperature. The maximum $\Delta E_L$ of up to -60 meV is predicted for Al owing to the very large coupling strengths, especially for gaps on the order of ≲10 Å. For Mg and Na the expected GS modifications are smaller by a factor of 2 due to a smaller coupling strength which is not fully offset by their respective lower transition energies. While these numbers are small in comparison to the unperturbed ground state energies, they are obtained for a single molecule coupled to an optical cavity. Although our cavities are model systems, they are still based on ab initio TDDFT and representative for what occurs in particle-on-mirror geometries²⁴ or picocavities³⁶ and these results can provide estimates for the order of magnitude of USC modifications to the ground state energy landscape. One of the relevant impacts of these changes is modification of ground state chemical reactivity via USC⁴² and correspondence to vdW forces.⁴³ We thus look into the energy scales of these interactions, keeping in mind that USC modifications calculated here refer to a single optical mode correction, while vdW accounts for all possible modes and polarizations.

We compare USC modifications and vdW for Mg₂₀₁ dimers with a single tetracene molecule. In addition to the USC GS modification, we calculate the total energy contributions as a function of three parameters: gap size and tetracene rotation angle (see Figure S5 for g). The largest contribution is calculated with conventional (static) density functional theory (DFT) using the *semilocal* PBE⁴⁴ exchange-correlation functional. As PBE does not include vdW interactions, which are associated with *nonlocal* correlation,⁴⁵ we evaluate their contributions to the energy by calculating dispersion corrections using the DFT-D3 method.⁴⁶ We use DFT-D3 as opposed to a nonlocal functional such as vdW-DF⁴⁵ as it provides a simple means for quantifying the magnitude of the vdW contribution, enabling direct comparison with the magnitude of the USC contribution. The energy contributions from DFT-PBE and vdW are calculated as $\Delta E = E_{2Mg_{201}+\text{tetr}} - E_{2Mg_{201}} - E_{\text{tetr}}$ (Figure 7).

First, we evaluate the energy as a function of the gap size *d* and separate the individual contributions from DFT-PBE and vdW, as well as the USC estimation (Figure 7a). The sum of DFT-PBE ($\Delta E_{\text{DFT}}$) and vdW ($\Delta E_{\text{vdW}}$) shows that the Mg₂₀₁-tetracene-Mg₂₀₁ system is bound at *d* = 17.2 Å with a potential well of about 350 meV. The USC correction is small amounting to approximately 30 meV, which is only 1 order of magnitude smaller than the other two terms. Thus, with respect to distance, USC has a small impact on the total energy at equilibrium and if separated out, it would be responsible for a small shift of the potential well on the order of single picometers. However, the USC correction decays less rapidly than $\Delta E_{\text{DFT}}$ or $\Delta E_{\text{vdW}}$ with the gap size, reaching 20% of their joint value (Figure S6k). Thus, a single strongly coupled mode can account for a significant portion of the full vdW interaction.

Next, we consider the rotations of the tetracene molecule in the gap. The center of the molecule is fixed to the middle of the gap with *d* = 17.2 Å (Figure 7b,c). For rotations about both *y* and *z* (see inset of Figure 7b,c), $\Delta E_{\text{vdW}}$ is smallest when the molecule is aligned along the dimer axis. $\Delta E_{\text{DFT}}$, on the other hand, has different minima for the different rotations. For the *z*-axis rotation, the energy also shows a local maximum caused by corner hydrogen atoms moving closer to Mg during rotation. The $\Delta E_{\text{USC}}$, however, is an order of magnitude smaller at -30 meV at 0° rotation angle and can only slightly change the angular positions of the minima.
<figure><img src="image_15.png" /><figcaption>**Figure 7.** Comparison of predicted ground state energy modification in the single molecule USC regime based on the longitudinal quantum optical Hamiltonian (solid orange line) with other contributions to the energy landscape, namely DFT-PBE (dashed-dotted blue line) and vdW (dashed dark orange line), for a Mg<sub>201</sub> dimer coupled to tetracene as a function of (a) dimer gap size, (b) y-axis rotation angle of tetracene, and (c) z-axis rotation angle of tetracene. The gap size in (b,c) is 17.2 Å, the minimum of the energy landscape in (a). In general, the USC energy modification plays a small part in modifying the energy landscape, but due to the longitudinal coupling its contribution results in a more stable system.</figcaption></figure>

We also calculated the rotation-dependence for other gap sizes (Figure S6). Overall, the contribution of the single mode USC to the GS is small, typically between 5–20% of the total energy. The largest relative values occur for larger gaps when $\Delta E_{DFT}$ and $\Delta E_{vdW}$ are small, or for small gaps when $\Delta E_{DFT}$ and $\Delta E_{vdW}$ cancel each other out. However, overall the largest potential impact of USC on the complex energy landscape occurs for large g which take place for small gaps.

In the studied cases the USC modifications to the ground state can reach significant values. Hence, one may ask if this fraction could be larger in other nanoscale systems. Overall, it is a question of balancing repulsive and attractive forces vs distance, rotation angle, and potentially molecule number and confronting these with the achievable ultrastrong coupling. A simple answer is not straightforward as the number of parameters to investigate/optimize is very large and, to a large extent, all $\Delta E$ contributions are susceptible to any of the above parameter changes. A key point to be addressed further, is the possibility of decoupling the induced electric field responsible for USC from the structural influences on PBE and the rest of vdW energy changes. Such decoupling would allow one to modify USC semi-independently of $\Delta E_{PBE}$ and $\Delta E_{vdW}$ although presently it is unclear to what extent it would be achievable. However, based on these initial studies, enhancing USC to significantly modify the energy landscape of nanoscale systems in the single-molecule regime at the scale of 100 meV appears challenging.

**Relation between TDDFT and Longitudinal Hamiltonian.** In the previous section, we employed the coupling coefficients g obtained by TDDFT to estimate the shift in the ground state energy using the correction that is obtained from a simple Hopfield Hamiltonian. This assumes that the g values obtained by TDDFT carry the same physical meaning as the coupling coefficient that enters the Hopfield Hamiltonian. It is thus warranted to discuss the validity of this approximation.

TDDFT is an exact theory that captures the dynamics of any electronic Hamiltonian, where the electrons are coupled only longitudinally via the Coulomb potential.⁴⁷ To the extent that we can neglect transverse photonic modes in the nanoscale dimers considered here, our calculations therefore provide an exact model up to the exchange-correlation approximations used. This also implies that the cQED extension of DFT/TDDFT to transverse photons is not required.³¹,⁴¹

It is not straightforward to compare model optical Hamiltonians rooted in second quantization with TDDFT, which is a classical field theory derivable from a classical Lagrangian operating on auxiliary densities. There is, however, plenty of common ground, as the TDDFT equations of motion may be written using classical mechanics with Poisson brackets.⁴⁸ Also, it is customary to bosonize the Fermionic polarization, which allows one to approximate a polaritonic system as a quadratic bosonic Hamiltonian. Bohm and Pines⁴⁹ applied this approach to the collective degree of freedoms (DOFs) in metals, establishing the foundation of the modern RPA, while Hopfield⁵⁰ used it for the treatment of localized dielectric response. These two examples represent the two extremes: free electrons associated with a plasmonic current and confined localized polarization.

Thus, the common denominator is the symplectic structure and similar linear equations of motion derivable either with commutators or Poisson brackets. There is thus a relation via the RPA, which corresponds to the TDDFT-exchange correlation (XC) kernel being zero $f_{xc} = 0$ and allows one to write the RPA-TDDFT equations of motion in second quantization.

To recognize the relation between model optical Hamiltonians and TDDFT it is useful to recast the latter in a similar 2 × 2 form as the Hopfield Hamiltonian. A suitable starting point is the Casida (frequency space) formulation of TDDFT.⁵¹ In this approach the response of the system is obtained by solving the following equation (see Supplementary Note S3)

$$U\Omega^2 U^T = \Delta^2 + 2\sqrt{\Delta} K \sqrt{\Delta} \quad (4)$$

where U is an orthonormal matrix diagonalizing the Casida matrix, and $\Omega$ is a diagonal matrix containing Casida eigenvalues, $\Delta$ are the electron–hole excitation energies, and K represents the coupling between the excitations. To obtain an estimate of the magnitude of K, which corresponds to the self-polarization energy, neglecting dissipation (including hot carrier generation⁵²), we can transform the system such that each dipolar plasmonic mode is represented by a 1 × 1-block in the Casida form $\Omega_p^2 = \Delta_p^2 + 2\sqrt{\Delta_p} K_p \sqrt{\Delta_p}$. Since $\Delta_p$ is a weighted sum of very low energy electron–hole transitions contributing to the (collective) plasmon, $\Delta_p^2 \approx 0$, and one is left with

$$\Omega_p^2 \approx 2\sqrt{\Delta_p} K_p \sqrt{\Delta_p} \quad (5)$$
For nanoparticles in the size between approximately one hundred and 1000 atoms, the average of the Kohn–Sham excitations $\Delta_P$ is on the order of 1 eV and decreases with increasing size.²⁵,³⁵,⁵² The plasmon energy on the other hand is $\Omega_P = 3.7$ eV for $Ag^{35,52}$ and $\Omega_P = 7.8$ eV for Al.²⁵ This gives self-polarization energies $K_P$ of 7 and 30 eV for Ag and Al, respectively. Although the magnitude of $\Delta_P$ decreases for large nanoparticles, the plasmon energy remains large and finite. We can thus represent the coherent low energy KS excitations as a single plasmonic mode, which we can identify as being analogous to the cavity element of the Hopfield Hamiltonian. These modes can then be connected to lowest order via dipolar coupling, yielding an equivalent form to the Hopfield Hamiltonian. This provides the basis for the dipolar coupling of subsystems, specifically nanoparticles and molecules, where the individual components are already diagonalized TDDFT systems.⁵³

We note that in the bulk limit ($q \to 0$) $\Delta_P$ goes to zero linearly in $q$ while $K_P$ diverges as $1/q^2$. $\Omega_P$ does not, however, diverge due to the presence of noninteger occupation numbers (which were omitted in eq 5 above) and approaches the bulk plasmon. This is equivalent to the protection against the infrared divergence familiar from Hopfield Hamiltonians.³¹

In Supplementary Note S3, we start from the RPA second quantization Hamiltonian (excluding exchange) of two subsystems (for example two plasmonic nanoparticles or an ultrastrongly coupled molecule and a nanocavity), and diagonalize the subsystems. At this stage, the intrasystem correlation energy of the subsystems may be obtained. This demonstrates a key difference between the self-polarization in the transverse system, where it is introduced due to coupling, and the plasmonic self-polarization here, which occurs regardless of intersystem coupling and is not caused by it. With this setup, we heuristically demonstrate a mechanism for a similar “no-go” theorem for the superradiant phase transition in purely longitudinal plasmonic systems, as Bernardis et al. presented for a two level system for a particular gauge.⁵⁴ Here, we discuss this in terms of coupling of a single plasmon to a molecule. Most of the oscillator strength of the plasmon, which is due to coherent collective excitations of intraband electrons (~1 eV KS energy) in the subsystems, ends up in the plasmon of energy $\Omega_P$.

Strengthening of the interaction between a cavity and a molecule requires more oscillator strength [see eq 1]. As a result of the Thomas–Reiche–Kuhn-sum rule,⁵⁵,⁵⁶ more oscillator strength requires a larger number of occupied levels. This increases the electron density in a confined volume and pushes the plasmon energy upward due to self-polarization, thus protecting the systems from a superradiant phase transition upon introduction of intersystem coupling. This is demonstrated in the SI by making a canonical transformation to the subsystems, and observing their dipolar coupling at the new renormalized energies. This is analogous to the transverse cavity case, where the $g_c^2/\omega_m(\hat{a} + \hat{a}^\dagger)^2$ term protects from a superradiant phase transition,¹⁹ and equivalently, this positive shift (renormalization) may be incorporated in the cavity eigenvalues

$$ \omega'_{cav} \rightarrow \sqrt{\omega_{cav}^2 + \frac{4g_c^2\omega_{cav}}{\omega_{ex}}} \quad (6) $$

by means of a Bogoliubov rotation.¹⁵ Here, the canonical transformation to the TDDFT subsystems is in analogy to the Bogoliubov transformation (which is also a symplectic trans-

formation). After we have diagonalized the subsystems, i.e., obtained a strong dipolar plasmon, we diagonalize the full system and obtain the zero-point energy shifts corresponding to the vdW interaction. Thus, we note that the zero-point energy shifts due to plasmonic self-polarization should not be mistaken for intersystem effects. This is in contrast to transverse coupling, where the $A^2$ term is introduced upon the coupling itself (care must be taken here as well in a single cavity mode approximation, since in a cavity, compared to free vacuum, there is not only an extra cavity mode, but absence of free vacuum modes below the plasmonic cutoff of the cavity).

## CONCLUSIONS AND OUTLOOK

In this work, we have explored the possibility of reaching the ultrastrong coupling regime with single molecules by coupling them to idealistic optical plasmonic nanocavities. Our calculations indicate that single-molecule-based USC is viable across a broad spectral range, but requires balancing the size of the optical mode and its vacuum electric field with the size and transition dipole moment of the molecule. Herein, we set the plasmonic cavity sizes and their resulting QNM volumes to the range of 1 to 10 nm³. This ensured good matching with the size of the molecule, allowing for very efficient use of the enhanced vacuum field. Specifically, for $Mg_{201}$ and $Al_{147}$ dimers with the smallest possible gaps, the mode volume of coupled systems approaches 1 nm³, which is close to the size of the molecules. In fact, squeezing the mode volumes even further could become counterproductive for reaching single-molecule USC, since single-atom picocavities³⁴,³⁶ could couple to individual atoms in nearby molecules, rather than the entire molecule, and induce, e.g., local modifications of the subgroups of a molecule.

Simultaneously, sub-10 nm³ mode volumes are associated with very large vacuum fields, here calculated to be in excess of 10 V·nm⁻¹. It is only in such large fields that the USC regime can be reached with a single molecule. Furthermore, the ultrastrongly coupled dressed light-matter systems studied here clearly manifest effects that extend beyond two-level models. This pertains specifically to the modification of the cavity by higher-energy excitations of the molecule, which contribute to the dielectric response of the molecule. They focus the mode in the vicinity of the molecule and increase the coupling strength, inducing shifts of the resonant energies of the components of the system. For comparison of magnitudes, real electric fields of comparable amplitude are obtained with attosecond lasers with peak intensities on the order of or exceeding 1TW·cm⁻². Naturally, the vacuum fluctuations in ultrasmall cavities should not be mistaken for real fields, but it might be interesting to study their implications in the future with TDDFT beyond the linear response regime, including cQED in the computational framework.⁴¹

In summary, the good matching of the spatial dimensions of the cavity and coupled molecules as well as the cooperation of the vacuum field and transition dipole moment of the molecules can enable single-molecule USC across the entire studied spectrum. For Al nanoparticles the dominant contribution responsible for their USC to benzene comes from the very strong vacuum field, which can exceed 10 V·nm⁻¹. Mg and Na nanoparticles exhibit, respectively, gradually weaker vacuum fields and rely more on the increasing transition dipole moment of the molecule to reach USC. However, even for Na dimers the vacuum fields can easily exceed 2 V·nm⁻¹. These enable single-molecule coupling strengths on the order of 13% of the molecular excitation energy. In turn, the predicted USC ground
state energy modification estimated by using cQED reaches 30–60 meV for a single molecule, which is comparable to $k_B T$ at room temperature. Such ground state energy modifications could have significant implications for strong-coupling assisted chemistry and other material properties.

## METHODS

**DFT and TDDFT Calculations.** The DFT and TDDFT calculations were carried out using the PBE⁴⁴ exchange-correlation functional in the adiabatic limit. The photo-absorption spectra are calculated using the δ-kick technique⁵⁷ in the linear-response regime and employing the dipole approximation for light-matter interaction. The spectrum is presented as the dipole strength function that is equivalent to photoabsorption cross section apart from a constant multiplier. The default projector augmented-wave (PAW)⁵⁸ data sets and double-ζ polarized (dzp) basis sets provided in GPAW were used for Al, C, and H. The dzp basis set of Al includes diffuse 3p functions, which are important for describing plasmon resonances.⁵⁹ For Na and Mg the corresponding p-valence basis sets were used with only one and two, respectively, electrons in the 3s orbital considered explicitly, while the lower electrons were treated as a frozen core within PAW. This simplification had only a minute impact on the photoabsorption spectrum as verified against the larger PAW setups. In general, while the used basis sets might not be adequate for yielding numerical values at the complete-basis-set limit, they are expected to be sufficient for the purposes of the present work. For real-time TDDFT a grid spacing parameter of 0.3 Å was chosen to represent densities and potentials, and the molecules/particles were surrounded by a vacuum region of at least 6 Å. The Hartree potential was evaluated on a larger grid with at least 100 Å vacuum around the system and a coarser grid spacing of 1.2 Å, and subsequently refined to the original grid. For the time propagation, we used a time step of $\Delta t = 15$ as and total propagation time of at least $T = 30$ fs. The spectra were broadened using Lorentzian spectral broadening with 0.1 eV corresponding to a full width at half-maximum of 0.2 eV. The ground state total energy calculations were calculated using the finite difference mode with wave functions expanded on a real space grid with a mesh spacing of 0.2 Å and a vacuum region around the molecules of 8 Å. The contributions of vdW interactions were evaluated by calculating dispersion corrections as an add-on to standard DFT.⁴⁶ All individual metal nanoparticles and molecules were relaxed independently to the point where all atomic forces are below 0.05 eV/Å, the grid spacing was 0.2 Å. Once relaxed, the dimers and molecule-dimer systems were assembled, but not relaxed further.

**Fitting of Absorption Spectra.** All photoabsorption spectra were fitted with the velocity-coupled harmonic oscillator model

$$ \sigma_{\text{abs}}(\omega) = \text{Im} \left[ \frac{a\omega[\omega_{\text{ex}}^2 - (\omega + \frac{\gamma_{\text{ex}}}{2})^2]}{[\omega_{\text{pl}}^2 - (\omega + \frac{\gamma_{\text{pl}}}{2})^2)(\omega_{\text{ex}}^2 - (\omega + \frac{\gamma_{\text{ex}}}{2})^2) - 4g^2\omega^2} \right] \quad (7) $$

to obtain the coupling strength $g$ as well as the resonance positions and widths of plasmon and molecular excitation, $\omega_{\text{pl}}$, $\gamma_{\text{pl}}$, $\omega_{\text{ex}}$ and $\gamma_{\text{ex}}$ respectively. The model assumes that the entire oscillator strength is given by the uncoupled plasmon with amplitude $a$. This is a good approximation, as the molecular photoabsorption spectra are about 10 to 100 times less intense than those of the nanoparticle dimers. This assumption allows us to obtain coupling strengths directly from a single calculation without overfitting.

**Quasinormal Modes.** The QNMs of nanosphere dimers are calculated using an approach stemming from the QNMEig solver described in ref 39, which is here based on the Wave Optics Module of COMSOL, a commercial software implementing finite-element method for electromagnetic modeling. The solver allows for the finding of QNM frequencies and normalized QNM fields for dispersive materials with a Drude–Lorentz permittivity. Here, we use a fitted Drude model to match the absorption spectrum of single clusters calculated with TDDFT, effectively constructing an effective, size-dependent permittivity for each metal cluster/particle. The nanoparticle radii $r$ correspond to half of the maximum distance between centers of opposite sides of icosahedral (Na, Al) or truncated octahedra (Mg) clusters modeled with TDDFT. The parameters of the Drude permittivity models and nanosphere radii are given in Table 1. In the calculations we use an extra fine mesh setting for the simulation domain and override the settings for spherical domains representing nanoclusters so that the maximum mesh element size is $r/3$, while the minimal mesh element size is $r/6$. We calculate the first 8 eigenfrequencies around the bright dipole mode frequency, which is obtained by modeling nanosphere dimers using SMUTHI, a T-matrix method code.⁶⁰ The normalized fields are then used to calculate the mode volume⁸ for dipole positions in the center of the gap (center of molecule), shifted to 1/4th of the gap and at the position of the end of a molecule. Only the bright dipolar mode contributes, since higher order (e.g., quadrupole) modes are negligibly small for the considered sphere sizes, while the on-axis component of the mode field vanishes at the dimer axis for other dipole modes.

**Software Used.** DFT calculations were carried out using the GPAW package⁶¹,⁶² with localized basis sets (LCAO mode)⁶³ and with uniform real-space grids with the finite difference approximation. TDDFT calculations were conducted using the LCAO-RT-TDDFT implementation in GPAW.⁶⁴ The ASE library⁶⁵ was used for constructing the atomic structures. The NumPy,⁶⁶ SciPy,⁶⁷ and Matplotlib⁶⁸ Python packages, and Inkscape were used for processing data and generating figures.

## ASSOCIATED CONTENT

### Supporting Information

The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acsphotonics.2c00066.

Figures S1–S4: Additional modeling details, structures, calculation and fitting results, induced electric fields; Figures S5–S6: Additional energy dependencies on structure parameters; Figures S7–S8: QNM volume calculation results; Tables S1 and S2: Modeled systems and fitted coupling strengths; Supplementary Note S1: Additional discussion of Figure S4 on mode profiles and induced fields; Supplementary Note S2: Derivation of the longitudinal QED Hamiltonian; Supplementary Note S3: Canonical transformation of subsystems (PDF)

## AUTHOR INFORMATION

### Corresponding Author

Tomasz J. Antosiewicz – Faculty of Physics, University of Warsaw, PL-02-093 Warsaw, Poland; orcid.org/0000-0003-2535-4174; Email: tomasz.antisiewicz@fuw.edu.pl
## Authors

Mikael Kuisma - Department of Chemistry, University of Jyväskylä, FI-40014 Jyväskylä, Finland; orcid.org/0000-0001-8323-3405

Benjamin Rousseaux - Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris, F-75005 Paris, France

Krzysztof M. Czajkowski - Faculty of Physics, University of Warsaw, PL-02-093 Warsaw, Poland; orcid.org/0000-0001-9106-2837

Tuomas P. Rossi - Department of Applied Physics, Aalto University, FI-00076 Aalto, Finland

Timur Shegai - Department of Physics, Chalmers University of Technology, SE-412 96 Gothenburg, Sweden; orcid.org/0000-0002-4266-3721

Paul Erhart - Department of Physics, Chalmers University of Technology, SE-412 96 Gothenburg, Sweden; orcid.org/0000-0002-2516-6061

Complete contact information is available at: https://pubs.acs.org/10.1021/acsphotonics.2c00066

## Author Contributions

M.K. and B.R. performed the analytical analysis. K.M.C. performed the QNM calculations. T.P.R. and T.J.A. performed the TDDFT calculations and data analysis. T.S., P.E., and T.J.A. conceived the general approach. T.J.A. coordinated the project. All authors wrote the paper.

## Funding

We acknowledge financial support from the Swedish Research Council (VR Miljö, Grant No: 2016-06059), the Knut and Alice Wallenberg Foundation (Grant No: 2019.0140), the Polish National Science Center (projects 2019/34/E/ST3/00359 and 2019/35/B/ST5/02477). T.P.R. acknowledges support from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Grant Agreement No. 838996 and support from the Academy of Finland under the Grant No. 332429. T.J.A. acknowledges support from the Project HPC-EUROPA3 (INFRAIA-2016-1-730897), with the support of the EC Research Innovation Action under the H2020 Programme.

## Notes

The authors declare no competing financial interest.

The atomic structures, the photoabsorption spectra, and basic codes are available at 10.5281/zenodo.3985706. Other data that support the findings of this study are available from the authors upon reasonable request.

## ACKNOWLEDGMENTS

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 and by the Interdisciplinary Centre for Mathematical and Computational Modelling at University of Warsaw (Grant #G55-6).

## REFERENCES

(1) Baranov, D. G.; Wersäll, M.; Cuadra, J.; Antosiewicz, T. J.; Shegai, T. Novel Nanostructures and Materials for Strong Light-Matter Interactions. ACS Photonics 2018, 5, 24–42.

(2) Garcia-Vidal, F. J.; Ciuti, C.; Ebbesen, T. W. Manipulating matter by strong coupling to vacuum fields. Science 2021, 373, eabd0336.

(3) Schachenmayer, J.; Genes, C.; Tignone, E.; Pupillo, G. Cavity-Enhanced Transport of Excitons. Phys. Rev. Lett. 2015, 114, 196403.

(4) Krainova, N.; Grede, A. J.; Tsokkou, D.; Banerji, N.; Giebink, N. C. Polaron Photoconductivity in the Weak and Strong Light-Matter Coupling Regime. Phys. Rev. Lett. 2020, 124, 177401.

(5) Munkhbat, B.; Wersäll, M.; Baranov, D. G.; Antosiewicz, T. J.; Shegai, T. Suppression of photo-oxidation of organic chromophores by strong coupling to plasmonic nanoantennas. Sci. Adv. 2018, 4, eaas9552.

(6) Thomas, A.; Lethuillier-Karl, L.; Nagarajan, K.; Vergauwe, R. M. A.; George, J.; Chervy, T.; Shalabney, A.; Devaux, E.; Genet, C.; Moran, J.; Ebbesen, T. W. Tilting a ground-state reactivity landscape by vibrational strong coupling. Science 2019, 363, 615–619.

(7) Englund, D.; Faraon, A.; Fushman, I.; Stoltz, N.; Petroff, P.; Vucković, J. Controlling cavity reflectivity with a single quantum dot. Nature 2007, 450, 857–861.

(8) Sauvan, C.; Hugonin, J. P.; Maksymov, I. S.; Lalanne, P. Theory of the Spontaneous Optical Emission of Nanosize Photonic and Plasmon Resonators. Phys. Rev. Lett. 2013, 110, 237401.

(9) Lalanne, P.; Yan, W.; Vynck, K.; Sauvan, C.; Hugonin, J.-P. Light interaction with photonic and plasmonic resonances. Laser Photonics Rev. 2018, 12, 1700113.

(10) Wu, T.; Gurioli, M.; Lalanne, P. Nanoscale Light Confinement: the Q's and V's. ACS Photonics 2021, 8, 1522–1538.

(11) Franke, S.; Hughes, S.; Dezfouli, M. K.; Kristensen, P. T.; Busch, K.; Knorr, A.; Richter, M. Quantization of Quasinormal Modes for Open Cavities and Plasmonic Cavity Quantum Electrodynamics. Phys. Rev. Lett. 2019, 122, 213910.

(12) Franke, S.; Richter, M.; Ren, J.; Knorr, A.; Hughes, S. Quantized quasinormal-mode description of nonlinear cavity-QED effects from coupled resonators with a Fano-like resonance. Phys. Rev. Res. 2020, 2, 033456.

(13) Feist, J.; Fernández-Domínguez, A. I.; García-Vidal, F. J. Macroscopic QED for quantum nanophotonics: emitter-centered modes as a minimal basis for multiemitter problems. Nanophotonics 2020, 10, 477–489.

(14) Buhmann, S. Y. Casimir-Polder forces on atoms in the presence of magnetoelectric bodies. Ph.D. thesis, Friedrich-Schiller-Universität: Jena, 2007.

(15) Frisk Kockum, A.; Miranowicz, A.; De Liberato, S.; Savasta, S.; Nori, F. Ultrastrong coupling between light and matter. Nat. Rev. Phys. 2019, 1, 19–40.

(16) Ciuti, C.; Bastard, G.; Carusotto, I. Quantum vacuum properties of the intersubband cavity polariton field. Phys. Rev. B 2005, 72, 115303.

(17) Mahan, G. D. Van der Waals Forces in Solids. J. Chem. Phys. 1965, 43, 1569–1574.

(18) Li, X.; Mandal, A.; Huo, P. Cavity frequency-dependent theory for vibrational polariton chemistry. Nat. Commun. 2021, 12, 1315.

(19) Baranov, D. G.; Munkhbat, B.; Zhukova, E.; Bisht, A.; Canales, A.; Rousseaux, B.; Johansson, G.; Antosiewicz, T. J.; Shegai, T. Ultrastrong coupling between nanoparticle plasmons and cavity photons at ambient conditions. Nat. Commun. 2020, 11, 2715.

(20) Mueller, N. S.; Okamura, Y.; Vieira, B. G. M.; Juergensen, S.; Lange, H.; Barros, E. B.; Schulz, F.; Reich, S. Deep strong light-matter coupling in plasmonic nanoparticle crystals. Nature 2020, 583, 780–784.

(21) Wersäll, M.; Munkhbat, B.; Baranov, D. G.; Herrera, F.; Cao, J.; Antosiewicz, T. J.; Shegai, T. Correlative dark-field and photoluminescence spectroscopy of individual plasmon-molecule hybrid nanostructures in a strong coupling regime. ACS Photonics 2019, 6, 2570–2576.

(22) Tichauer, R. H.; Feist, J.; Groenhof, G. Multi-scale dynamics simulations of molecular polaritons: The effect of multiple cavity modes on polariton relaxation. J. Chem. Phys. 2021, 154, 104112.

(23) Yang, Z.-J.; Antosiewicz, T. J.; Shegai, T. Role of material loss and mode volume of plasmonic nanocavities for strong plasmon-exciton interactions. Opt. Express 2016, 24, 20373–20381.

(24) Chikkaraddy, R.; de Nijs, B.; Benz, F.; Barrow, S. J.; Scherman, O. A.; Rosta, E.; Demetriadou, A.; Fox, P.; Hess, O.; Baumberg, J. J. Single-
molecule strong coupling at room temperature in plasmonic nano-cavities. *Nature* **2016**, *535*, 127–130.

(25) Rossi, T. P.; Shegai, T.; Erhart, P.; Antosiewicz, T. J. Strong plasmon-molecule coupling at the nanoscale revealed by first-principles modeling. *Nat. Commun.* **2019**, *10*, 3336.

(26) Zhang, P.; Feist, J.; Rubio, A.; García-González, P.; García-Vidal, F. Ab initio nanoplasmonics: The impact of atomis structure. *Phys. Rev. B* **2014**, *90*, 161407.

(27) Barbry, M.; Koval, P.; Marchesin, F.; Esteban, R.; Borisov, A.; Aizpurua, J.; Sánchez-Portal, D. Atomistic Near-Field Nanoplasmonics: Reaching Atomic-Scale Resolution in Nanooptics. *Nano Lett.* **2015**, *15*, 3410–3419.

(28) Marchesin, F.; Koval, P.; Barbry, M.; Aizpurua, J.; Sánchez-Portal, D. Plasmonic Response of Metallic Nanojunctions Driven by Single Atom Motion: Quantum Transport Revealed in Optics. *ACS Photonics* **2016**, *3*, 269–277.

(29) García-Ripoll, J. J.; Peropadre, B.; De Liberato, S. Light-matter decoupling and A(2) term detection in superconducting circuits. *Sci. Rep.* **2015**, *5*, 16055–16055.

(30) Rokaj, V.; Welakuh, D. M.; Ruggenthaler, M.; Rubio, A. Light-matter interaction in the long-wavelength limit: no ground-state without dipole self-energy. *J. Phys. B: At., Mol. Opt. Phys.* **2018**, *51*, 034005.

(31) Schäfer, C.; Ruggenthaler, M.; Rokaj, V.; Rubio, A. Relevance of the Quadratic Diamagnetic and Self-Polarization Terms in Cavity Quantum Electrodynamics. *ACS Photonics* **2020**, *7*, 975–990.

(32) Rokaj, V.; Penz, M.; Sentef, M. A.; Ruggenthaler, M.; Rubio, A. Quantum Electrodynamical Bloch Theory with Homogeneous Magnetic Fields. *Phys. Rev. Lett.* **2019**, *123*, 047202.

(33) Antosiewicz, T. J.; Apell, S. P.; Claudio, V.; Käll, M. A simple model for the resonance shift of localized plasmons due to dielectric particle adhesion. *Opt. Express* **2012**, *20*, 524–533.

(34) Urbieta, M.; Barbry, M.; Zhang, Y.; Koval, P.; Sánchez-Portal, D.; Zabala, N.; Aizpurua, J. Atomic-Scale Lightning Rod Effect in Plasmonic Picocavities: A Classical View to a Quantum Effect. *ACS Nano* **2018**, *12*, 585–595.

(35) Rossi, T. P.; Kuisma, M.; Puska, M. J.; Nieminen, R. M.; Erhart, P. Kohn-Sham decomposition in real-time time-dependent density-functional theory: An efficient tool for analyzing plasmonic excitations. *J. Chem. Theory Comput.* **2017**, *13*, 4779–4790.

(36) Carnegie, C.; Griffiths, J.; de Nijs, B.; Readman, C.; Chikkaraddy, R.; Deacon, W. M.; Zhang, Y.; Szabó, I.; Rosta, E.; Aizpurua, J.; Baumberg, J. J. Room-Temperature Optical Picocavities below 1 nm3 Accessing Single-Atom Geometries. *J. Phys. Chem. Lett.* **2018**, *9*, 7146–7151.

(37) Zhang, C.; Chen, B.-Q.; Li, Z.-Y. Optical origin of subnanometer resolution in tip-enhanced Raman mapping. *J. Phys. Chem. C* **2015**, *119*, 11858–11871.

(38) Bai, Q.; Perrin, M.; Sauvan, C.; Hugonin, J. P.; Lalanne, P. Efficient and intuitive method for the analysis of light scattering by a resonant nanostructure. *Opt. Express* **2013**, *21*, 27371–27382.

(39) Yan, W.; Faggiani, R.; Lalanne, P. Rigorous modal analysis of plasmonic nanostructures. *Phys. Rev. B* **2018**, *97*, 205422.

(40) Dezfouli, M. K.; Tserkezis, C.; Mortensen, N. A.; Hughes, S. Nonlocal quasinormal modes for arbitrarily shaped three-dimensional plasmonic resonators. *Optica* **2017**, *4*, 1503–1509.

(41) Flick, J.; Welakuh, D. M.; Ruggenthaler, M.; Appel, H.; Rubio, A. Light-Matter Response in Nonrelativistic Quantum Electrodynamics. *ACS Photonics* **2019**, *6*, 2757–2778.

(42) Martínez-Martínez, L. A.; Ribeiro, R. F.; Campos-González-Angulo, J.; Yuen-Zhou, J. Can Ultrastrong Coupling Change Ground-State Chemical Reactions? *ACS Photonics* **2018**, *5*, 167–176.

(43) Haugland, T. S.; Schäfer, C.; Ronca, E.; Rubio, A.; Koch, H. Intermolecular interactions in optical cavities: An ab initio QED study. *J. Chem. Phys.* **2021**, *154*, 094113.

(44) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. *Phys. Rev. Lett.* **1996**, *77*, 3865–3868.

(45) Hyldgaard, P.; Jiao, Y.; Shukla, V. Screening nature of the van der Waals density functional method: a review and analysis of the many-body physics foundation. *J. Phys.: Condens. Matter* **2020**, *32*, 393001.

(46) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. *J. Chem. Phys.* **2010**, *132*, 154104.

(47) Runge, E.; Gross, E. K. U. Density-Functional Theory for Time-Dependent Systems. *Phys. Rev. Lett.* **1984**, *52*, 997–1000.

(48) Chernyak, V.; Mukamel, S. Density-matrix representation of nonadiabatic couplings in time-dependent density functional (TDDFT) theories. *J. Chem. Phys.* **2000**, *112*, 3572–3579.

(49) Bohm, D.; Pines, D. A Collective Description of Electron Interactions: III. Coulomb Interactions in a Degenerate Electron Gas. *Phys. Rev.* **1953**, *92*, 609–625.

(50) Hopfield, J. Theory of the contribution of excitons to the complex dielectric constant of crystals. *Phys. Rev.* **1958**, *112*, 1555.

(51) Casida, M. E. Time-dependent density-functional theory for molecules and molecular solids. *J. Mol. Struct.: THEOCHEM* **2009**, *914*, 3–18.

(52) Rossi, T. P.; Erhart, P.; Kuisma, M. Hot-Carrier Generation in Plasmonic Nanoparticles: The Importance of Atomic Structure. *ACS Nano* **2020**, *14*, 9963–9971.

(53) Fojt, J.; Rossi, T. P.; Antosiewicz, T. J.; Kuisma, M.; Erhart, P. Dipolar coupling of nanoparticle-molecule assemblies: An efficient approach for studying strong coupling. *J. Chem. Phys.* **2021**, *154*, 094109.

(54) De Bernardis, D.; Pilar, P.; Jaako, T.; De Liberato, S.; Rabl, P. Breakdown of gauge invariance in ultrastrong-coupling cavity QED. *Phys. Rev. A* **2018**, *98*, 053819.

(55) Kuhn, W. Über die Gesamtstärke der von einem Zustande ausgehenden Absorptionslinien. *Z. Phys.* **1925**, *33*, 408–412.

(56) Reiche, F.; Thomas, W. Über die Zahl der Dispersionselektronen, die einem stationären Zustand zugeordnet sind. *Z. Phys.* **1925**, *34*, 510–525.

(57) Yabana, K.; Bertsch, G. F. Time-dependent local-density approximation in real time. *Phys. Rev. B* **1996**, *54*, 4484–4487.

(58) Blöchl, P. E. Projector augmented-wave method. *Phys. Rev. B* **1994**, *50*, 17953–17979.

(59) Rossi, T. P.; Lehtola, S.; Sakko, A.; Puska, M. J.; Nieminen, R. M. Nanoplasmonics simulations at the basis set limit through completeness-optimized, local numerical basis sets. *J. Chem. Phys.* **2015**, *142*, 094114.

(60) Egel, A.; Czajkowski, K. M.; Theobald, D.; Ladutenko, K.; Kuznetsov, A. S.; Pattelli, L. SMUTHI: A python package for the simulation of light scattering by multiple particles near or between planar interfaces. *J. Quant. Spectrosc. Radiat. Transfer* **2021**, *273*, 107846.

(61) Mortensen, J.; Hansen, L.; Jacobsen, K. Real-space grid implementation of the projector augmented wave method. *Phys. Rev. B* **2005**, *71*, 035109.

(62) Enkovaara, J.; et al. Electronic structure calculations with GPAW: a real-space implementation of the projector augmented-wave method. *J. Phys.: Condens. Matter* **2010**, *22*, 253202.

(63) Larsen, A.; Vanin, M.; Mortensen, J.; Thygesen, K.; Jacobsen, K. Localized atomic basis set in the projector augmented wave method. *Phys. Rev. B* **2009**, *80*, 195112.

(64) Kuisma, M.; Sakko, A.; Rossi, T.; Larsen, A.; Enkovaara, J.; Lehtovaara, L.; Rantala, T. Localized surface plasmon resonance in silver nanoparticles: Atomistic first-principles time-dependent density-functional theory calculations. *Phys. Rev. B* **2015**, *91*, 115431.

(65) Hjorth Larsen, A.; Jørgen Mortensen, J.; Blomqvist, J.; Castelli, I. E.; Christensen, R.; Dulak, M.; Friis, J.; Groves, M. N.; Hammer, B.Ør.; Hargus, C.; Hermes, E. D.; Jennings, P. C.; Bjerre Jensen, P.; Kermode, J.; Kitchin, J. R; Leonhard Kolsbjerg, E.; Kubal, J.; Kaasbjerg, K.; Lysgaard, S.; Bergmann Maronsson, J.; Maxson, T.; Olsen, T.; Pastewka, L.; Peterson, A.; Rostgaard, C.; Schøtz, J.; Schutt, O.; Strange, M.; Thygesen, K. S.; Vegge, T.; Vilhelmsen, L.; Walter, M.; Zeng, Z.; Jacobsen, K. W The Atomic Simulation Environment – A
Python library for working with atoms. *J. Phys.: Condens. Matter* **2017**,
29, 273002.

(66) Harris, C. R.; et al. Array programming with NumPy. *Nature* **2020**, *585*, 357–362.

(67) Virtanen, P.; et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. *Nat. Methods* **2020**, *17*, 261–272.

(68) Hunter, J. D. Matplotlib: A 2D graphics environment. *Comput. Sci. Eng.* **2007**, *9*, 90–95.