Ultrastrong coupling between nanoparticle plasmons and cavity photons at ambient conditions
Results .Ultrastrong coupling in plasmon–microcavity systems . The system under study is illustrated in Fig.? 1a . It consists of a sub-diffractive periodic array of gold (Au) nanorods placed at?the antinode of the fundamental Fabry–Pérot (FP) microcavity mode formed by two Au?mirrors and filled by a SiO 2 spacer. The nanorod array couples to the vacuum field of the FP microcavity, thus producing plasmon–cavity polaritons manifested as distinct resonant spectral features emerging in transmission, reflection, and absorption spectra of the coupled system.Fig. 1: Sketch and numerical modeling of the coupled system.a Artistic illustration of the system: an array of plasmonic nanorods positioned in the middle of a Fabry–Pérot cavity formed by two gold mirrors. The cavity interior is filled with SiO 2 . The array couples to the FP cavity mode, exchanging energy at a rate g . b False-color normal-incidence absorption spectra as a function of cavity thickness with an array of 300?nm long plasmonic nanorod (width 50?nm, height 20?nm) positioned in the middle of SiO 2 -filled Fabry–Pérot cavity. The vertical dashed line indicates the nanorod plasmon resonance outside of the cavity. The curved lines indicate resonances of the empty FP cavity, whose even modes are not modified by the coupling. ΩRdenotes plasmon–cavity mode splitting at zero detuning. c The electric field intensity (in the log scale) and the electric field lines in the vertical plane across the middle of the nanorod induced by a normally incident plane wave (polarized in the figure plane) for the coupled system of 400?nm thick cavity and 300?nm long nanorods calculated for the lower and upper polaritons. Full size image To provide initial insight into the behavior of the coupled system, we perform numerical finite-difference time-domain (FDTD) simulations (FDTD Solutions, Lumerical). Figure? 1b shows a map of absorption spectra of coupled FP–nanorod systems at normal incidence with the electric field parallel to the nanowires as a function of the cavity thickness for nanorod lengths L ?=?300?nm and dy ?=?30?nm spacing in the y -direction. For an easy comparison between these coupled system spectra with the uncoupled elements, we plot the bare FP cavity resonances with curved lines. The vertical dashed line marks the bare plasmon nanorod resonance of the array. A comparison clearly shows a rather complicated picture of new eigenmodes’ dispersion in which the even FP modes are practically unperturbed while the odd FP modes are shifted significantly from the bare cavity positions.The 1st order FP mode of an empty cavity intersects the bare nanorod array plasmon resonance around 400?nm cavity thickness resulting in a distinct anticrossing (Fig.? 1b ). The lower polariton (LP) transitions from a plasmon-dominated mode (for a thin cavity) to an FP-dominated mode at large detuning (for a thick cavity). However, the upper polariton (UP) upon acquiring a plasmon-like character at large detuning, crosses the 2nd order FP mode and approaches the spectral position of the 3rd FP cavity mode, which in the coupled system is strongly pushed to the blue due to hybridization with the plasmon. Such qualitative behavior is observed for all the odd FP modes: each odd coupled i th mode is pushed to the blue beyond the subsequent even i ?+?1st mode (which is unperturbed) and approaches the i ?+?2nd odd FP mode. In fact, at no point in the spectral analysis do any of the FP–plasmon–polaritons follow the plasmon dispersion. In contrast, the even modes in the coupled system do not significantly interact with the array because they have a node of the electric field in the center of the cavity, where the rods are positioned. These observations suggest that a multimode character of the FP microcavity is important for a detailed interpretation of our results.Another remarkable feature of the absorption map in Fig.? 1b is the dispersion of the lower polariton in the thin cavity limit: for cavities thinner than about 200?nm, the LP dispersion exhibits a back-bending to extremely low energies, Fig.? 1b , which is likely related to the near-field interaction of discrete plasmonic nanoparticles with the cavity mirrors in this short-range limit. This behavior is not reproduced by the Hamiltonian modeling, and we will not consider it in detail in the following.The spatial distributions of the electric field induced by a normally incident plane wave inside the plasmon–cavity system calculated at the resonant energies for a 400?nm thick cavity (Fig.? 1c ) clearly display the opposite symmetries of the two resonances. While the lower energy mode shows an anti-symmetric combination of cavity and plasmon fields, featuring two saddle points above and below the nanorod, the upper energy mode is a symmetric combination. Such behavior highlights the polaritonic nature of the two resonances of the hybrid system. For a 400?nm thick cavity, corresponding to near-resonant coupling ( ω cav?=? ω pl ??0.8?eV), the Rabi splitting, ΩR, estimated as the energy difference between the two absorption peaks reaches ~1?eV. Thus, assuming that ΩR?=?2 g on resonance, we estimate the normalized coupling strength of g / ω pl?>?0.5, which clearly indicates the ultrastrong coupling regime in the system. In what follows, we perform a more rigorous estimation of the g / ω plvalues in our systems based on a full Hopfield Hamiltonian.The same qualitative behavior is observed for diluted arrays and ones with longer nanorods, as shown in Supplementary Fig.? 1 . Although the coupling strength decreases with a smaller nanorod density, the spectral behavior indicates that all the odd modes mix with the plasmon, yielding a complex polaritonic system. However, for longer rods and/or diluted arrays, in addition to the plasmon, higher-order resonances are observed, which are associated with lattice modes. Additional data, including transmission and reflection spectra and complete data for the incident light polarization perpendicular to the nanorod axis, as well as the uncoupled cavity and array elements, are provided in Supplementary Figs.? 1 – 5 .Samples of coupled plasmon–microcavity systems were fabricated by combination of electron beam evaporation (Au mirrors), plasma-enhanced chemical vapor deposition (dielectric spacers), and electron beam lithography (nanorod arrays) (see “Methods” for details). Figure? 2a shows a bright-field optical microscope image of the fabricated nanorod arrays with lengths ranging from 200 to 400?nm with a step of 50?nm. The nanorods have fixed height of h ?=?20?nm and width of w ?=?50?nm. An exemplary scanning electron microscope (SEM) image of gold nanorods array with length of Lrod ?=?250?nm is shown in Fig.? 2b (see “Methods” for details). Both the figures clearly show high-density plasmonic arrays with an interparticle distance as small as 30?nm, corresponding to the surface filling factor of 60%. More examples are shown in Supplementary Fig.? 6 .Fig. 2: Fabricated samples.a Bright-field optical microscope images of gold nanorod arrays positioned in the middle of a SiO 2 -filled FP cavity (without the top mirror) fabricated by electron beam lithography. Individual nanorods have a fixed height of h ?=?20?nm, width of w ?=?50?nm, and length varying from 200 to 400?nm. The side-to-side distance between the nanorods is 30?nm. The arrays are 250?×?250?μm 2 . b SEM image of the Lrod ?=?250?nm nanorods array. The inset shows a magnified view of the nanorod array. Full size image Next, we proceed to optical measurements of the fabricated plasmon–cavity systems using Fourier transform infrared (FTIR) spectroscopy (see “Methods” for the details of measurements). Figure? 3a, b shows exemplary, measured at normal incidence, reflection, and absorption spectra of an empty 400?nm thick cavity, 300?nm long nanorods array, and those of the coupled system (see Supplementary Figs.? 7 and 8 for measured reflection and absorption spectra of all uncoupled cavities and nanorods). The uncoupled cavity and array resonances overlap spectrally and, when coupled, unambiguously confirm the realization of a giant Rabi splitting in the spectra of the coupled plasmon–cavity systems. We note that the 2nd order Fabry–Pérot mode redshifts in the hybrid system although it cannot interact with the nanorod array located exactly in the middle of the cavity. This behavior is consistently observed for all measured coupled systems, Supplementary Fig.? 9 . It could be explained by the hybrid cavities having slightly larger thickness due to the presence of the plasmonic array that could redshift all the uncoupled Fabry–Pérot modes.Fig. 3: Measurements of the fabricated samples.a , b Measured reflection ( a ) and absorption ( b ) spectra of an empty Lcav ?=?400?nm cavity, bare Lrod ?=?300?nm long plasmonic nanorods, and those of the coupled system with the electric field polarization parallel to the major rod axis. c Measured dispersion of the reflection spectra of the coupled plasmon–cavity system with Lrod ?=?300?nm plasmonic nanorods as a function of the cavity thickness revealing an anti-crossing between the two polaritonic modes. Dashed lines show the positions of the bare plasmon array mode and the bare Fabry–Pérot modes of the 1st and 2nd order. ΩRdenotes the minimal observed splitting between the two polaritonic modes. Full size image Dispersion of measured normal-incidence reflection spectra from coupled systems with 300?nm long nanorods and varying cavity thickness displays a clear anticrossing between the 1st order Fabry–Pérot mode and the plasmon mode of the array, Fig.? 3c (see Supplementary Fig.? 9 for dispersions of reflection and absorption spectra vs cavity thickness for all nanorod lengths). The spectra also reveal the 2nd order Fabry–Pérot mode (third dip from the left), which does not interact with the nanorods due to the electric field node in the center of the cavity. As revealed by the Hamiltonian analysis of the spectra in the next section, the nanorods array mode additionally redshifts from ~0.8 to ~0.7?eV due to the presence of a dielectric medium around the rods. For the thinnest 100?nm thick plasmon-loaded cavities, neither reflection nor absorption spectra show the exact position of the upper polariton, which was beyond the detection range of the FTIR microscope used for this set of measurements. For this reason, we performed additional reflectivity measurements in the visible range for the 100?nm thick samples using a normal optical microscope to capture the spectral feature of the upper polariton (see Supplementary Fig.? 10 ). Additional reflection spectra at normal incidence in the visible range were collected using a 20× objective (Nikon, NA?=?0.45), directed to a fiber-coupled spectrometer and normalized with reflection from a standard dielectric-coated silver mirror. Based on these spectra, the vacuum Rabi splitting taken as the energy difference between the two reflection dips at zero detuning ( ω cav?=? ω pl, 500?nm thick cavity), reaches ~0.8?eV at the resonant energy of ~0.7?eV, Fig.? 3c . Thus, the Rabi splitting in our samples exceeds both the bare cavity and bare plasmon resonance frequencies, indicating that the hybrid plasmon–cavity system is deep into the USC regime.Analysis of the ultrastrong coupling using Hopfield Hamiltonian . We now turn to a more thorough analysis of the experimental data. Since a rough estimation already reveals that the Rabi splitting in our system is comparable to the transition energy of uncoupled oscillators, the usual Jaynes–Cummings or Rabi-type coupled Hamiltonians are invalid, and a more general Hamiltonian must be used. Therefore, to analyze our system we employ the full Hopfield Hamiltonian including both the fast-rotating and the quadratic A2 terms, which capture the essential physical characteristics of an ultrastrongly coupled system9. We will focus on the two lowest modes of the plasmon–cavity structure, hence we will consider only coupling of two oscillators: the 1st order normal incidence FP mode of the cavity?with energy \(\hbar \omega_{cav}\), and the collective long-axis plasmon mode of the array?with energy \(\hbar \omega_{pl}\). Here, the cavity mode plays the role of the light component of the system, whereas the plasmonic nanorod array mode plays the role of the matter component. The total Hamiltonian thus reads: $$\hat H = \hbar \omega _{cav}\left( {\frac{1}{2} + \hat a^\dagger \hat a} \right) + \hbar \omega _{pl}\left( {\frac{1}{2} + \hat b^\dagger \hat b} \right) + \hat H_{int},$$ (1)where \(\hat a\)and \(\hat b\)are the microcavity and collective plasmon annihilation operators, respectively, and \(\hat H_{int}\)is the interaction Hamiltonian. If we were to consider individual nanoparticle plasmons interacting with each other instead of the collective array mode, the Hamiltonian would also yield additional eigenstates weakly interacting with light16. As long as we work away from the Rayleigh modes of the array30, which is ensured by sub-diffraction periodicity, all the plasmon–plasmon interaction effects can be absorbed into the collective plasmon frequency ω pl.The interaction part can be written differently depending on the gauge in which the electromagnetic field is treated. The two options that are often used are the Coulomb gauge and its dipole representation. The latter can be obtained from the Coulomb gauge by performing the Power–Zienau–Woolley transformation31. When a cavity couples to a two-level system, the two representations are not gauge-invariant because of the two-level approximation32 , 33. However, since we are considering coupling of two harmonic oscillators, the two pictures provide identical spectra14 , 16. We?will therefore use the Coulomb gauge, in which the single-mode interaction Hamiltonian can be written as34 , 35: $$\hat H_{int} = \hbar g_C\left( {\hat a^\dagger + \hat a} \right)\left( {\hat b^\dagger + \hat b} \right) + \frac{{\hbar g_C^2}}{{\omega _{pl}}}\left( {\hat a^\dagger + \hat a} \right)^2,$$ (2)where \(\hbar g_C = \mu _{pl}\sqrt {a^2\rho } {\cal{E}}_{vac}\frac{{\omega _{pl}}}{{\omega _{cav}}}\)is the coupling strength with μ plbeing the transition dipole moment of the plasmonic nanorod, ρ the plasmonic nanoparticles density per unit area a2 ( \(\sqrt {a^2\rho }\)thus has a familiar \(\sqrt N\)scaling), and \({\cal{E}}_{vac} = \sqrt {\frac{{\hbar \omega _{cav}}}{{2\varepsilon \varepsilon _0a^2L_{\mathrm{eff}}}}}\)the vacuum electric field of the cavity with Leff being the effective cavity mode transverse thickness16. The first term in Eq. ( 2 ) is the usual Rabi-type interaction including both slow and fast-rotating terms. The second term is the so-called A2 term, which arises from the expansion of the minimal coupling Hamiltonian \(\left( {{\mathbf{p}} - \frac{e}{c}{\mathbf{A}}} \right)^2\)and “protects” the coupled system from the superradiant phase transition12 , 13, as well as stabilizes the spectrum against the square-root singularity36. Supplementary Figure? 12 shows the spectrum of Hamiltonian ( 1 ) for ω pl?=? ω cav?=?1?eV as a function of the coupling constant g C; it also demonstrates that neglecting the A2 term as well as fast-rotating terms leads to incorrect and unphysical spectra.We want to emphasize that although our system is essentially classical, we choose to use the quantum Hamiltonian because it provides a convenient description in terms of the modes amplitudes via creation operators from start. Moreover, the use of the quantum Hamiltonian allows us to obtain the characteristics of the ground state of the system in a straightforward way, which we analyze in the following. However, the linear response and the energy spectrum of the system can be equally obtained from a classical description not involving any operator algebra, which is demonstrated above by the FDTD simulations.In a classical optical experiment, the outcome of which is some response function of the system, such as elastic scattering, reflection, or absorption, one cannot access directly the ground-state energy. However, spectral positions of the resonant features in reflection or absorption spectra reflect approximately the transition energies between the ground and first excited states of the system \(\hbar \omega _ \pm = E_{ \pm 1} - E_0\). Therefore, to model the system with the Hopfield Hamiltonian framework, we fit the measured dispersions of reflection dips with calculated transition energies \(\hbar \omega _ \pm\)of the Hopfield Hamiltonian34.The resulting Hamiltonian fit of a coupled system’s resonant transitions as a function of the bare cavity energy is presented in Fig.? 4a for Lrod ?=?300?nm nanorod arrays. For each cavity thickness, the bare cavity energy was determined from the spectral position of its reflection dip (Supplementary Fig.? 7 ). By assuming that the effective cavity thickness scales as \(L_{\mathrm{eff}} = \frac{{\lambda _{cav}}}{{4n}}\)with n being refractive index of the cavity medium, we arrive at the coupling strength in the Coulomb gauge \(g_C = \omega _{pl}\mu _{pl}\sqrt {\frac{{\hbar \rho }}{{\pi \varepsilon _0nc}}}\), which is independent of the cavity's thickness and energy. Hence, we fit the polaritonic dispersion by freely varying plasmon frequency ω pland the coupling strength g C. For the Lrod ?=?300?nm nanorod arrays, the fitting yields the plasmon frequency of 640?meV and the coupling strength of 300?meV, resulting in Rabi splitting of exactly 2 g C?=?600?meV at resonance, ω pl?=? ω cav(see Supplementary Fig.? 13 for Hamiltonian fits of other coupled systems, and Supplementary Table? I for extracted plasmon energies and coupling strengths). For all five nanorod lengths, we consistently obtain normalized coupling strength values g C/ ω plin the range from 0.4 to 0.56, Fig.? 4b , which unambiguously indicate the USC regime of interaction between the nanorods and the cavity modes9. Furthermore, we notice that the normalized coupling strength \(g_C/\omega _{pl} = \mu _{pl}\sqrt {\frac{{\hbar \rho }}{{\pi \varepsilon _0nc}}}\)is a function of the plasmon transition?dipole moment and the particle density only. Therefore, if the product \(\mu _{pl}\sqrt \rho\)grows with increasing nanorod length, we may expect even higher values of g C/ ω plfor longer rods resonating at lower energies.Fig. 4: Analysis of the experimental data.a Fitting of the measured polaritonic dispersion of the coupled plasmon–cavity system ( Lrod ?=?300?nm) with Hopfield Hamiltonian transition energies. Dots show resonant energies of the coupled system extracted as experimental reflection dips, lines are Hopfield polaritons dispersion, gray dashed lines are the bare cavity and bare plasmon energies. ΩRdenotes Rabi splitting between the two polaritonic modes at the zero-detuning point. b Normalized coupling strength g C/ ω plat zero detuning versus nanorod length. Solid line is the linear interpolation. Full size image We also compare the resulting fits with those obtained by applying the multimode Hopfield Hamiltonian accounting for all the normal-incidence modes of a Fabry–Pérot cavity, which can be solved analytically35(Supplementary Note? 1 ). The most prominent difference is that in the multimode picture the upper polariton crosses the bare plasmon energy exactly at the point where it also crosses the even cavity mode (see Supplementary Fig.? 14 ), in perfect agreement with FDTD simulations (Fig.? 1b ). Nevertheless, the resulting fits and coupling strengths are very close to the results of the single-mode Hamiltonian analysis (see Supplementary Table? II ), which justifies its validity in our case. Of course, the single-mode Hamiltonian only yields correct spectra as long as all higher-order (odd) cavity modes are detuned from the array mode, which is the case in our experimental configuration. The multimode Hamiltonian picture, however, may become important in other cases, as was also mentioned in the discussion of numerical FDTD results in Fig.? 1b .It is further instructive to compare the obtained values with an estimation for g Cthat can be deduced directly from the geometry of the system. The vacuum electric field can be calculated as \({\cal{E}}_{vac} = \sqrt {\frac{{\hbar \omega _{cav}}}{{2\varepsilon \varepsilon _0a^2L_{\mathrm{eff}}}}}\)where for the effective cavity mode thickness one can use \(L_{\mathrm{eff}} = \frac{{\lambda _{cav}}}{{4n}}\). The plasmon transition dipole moment can be estimated from the scattering cross-section of a single nanorod in free space (see Supplementary Note? 2 ) and by applying the classical Larmor formula for the decay rate of a point dipole37. This yields the value of around 3.4?×?10 4 Debye for the 400?nm long Au nanorod (corresponding to the radiative decay rate of ~67?meV, see Supplementary Table? III ). Combining these values with the nanorod density ρ ?=?(430?nm ?80?nm) ?1 , one obtains the resonant ( ω pl?=? ω cav) coupling strength of around g C?≈?0.3?eV, which agrees perfectly with the results of fitting.We further illustrate the importance of keeping the quadratic term by analyzing the data with simpler, albeit a priori incorrect, Hamiltonians. An attempt to fit the experimental with eigenvalues of Hopfield Hamiltonian without the A2 term does not yield any reasonable result with a region of the energy spectrum becoming imaginary, Supplementary Fig.? 17 , and slightly overestimated coupling strengths. This imaginary spectrum is a fundamental property of the coupled oscillators Hamiltonian without any kind of quadratic stabilizing term36 , 38. Fitting the data with no A2 Hopfield Hamiltonian under RWA (i.e., also without fast-rotating terms), although seems to give a better fit, yields regions with negative LP energy spectrum and largely overestimated coupling strength, Supplementary Fig.? 17 .Besides dispersions of polaritonic energies, a remarkable feature of all measured reflection spectra is that the lower (and even upper) polaritons are much narrower than the bare plasmon mode (see, e.g., Fig.? 3a and Supplementary Fig.? 9 ). This behavior contrasts the non-Hermitian Jaynes–Cummings and Rabi Hamiltonians, where the decay rates of the photonic mode and electronic transition are “shared” equally at zero detuning between lower and upper polaritons: Im ω± ?=??( γ cav?+? γ x)/4, where γ xis the linewidth of the electronic transition3 , 4. A similar suppression of the polariton linewidth has been observed with a cyclotron resonance coupled to a Fabry–Pérot cavity39. While we cannot directly introduce a non-Hermitian part into Hamiltonian ( 1 ), since it will render the ground-state energy complex-valued, we can qualitatively describe the decay processes in the coupled system. As we showed above, the linewidth of a single plasmonic nanorod outside the cavity is largely determined by its radiative loss. For nanorod arrays outside the cavity the radiative decay is even more dominant, reaching ~95% of the total linewidth (see Supplementary Note? 2 ). However, when the dipolar oscillator, such as our plasmonic array, is placed between the mirrors, its radiation does not instantaneously leave the system—instead, it first bounces between the mirrors, and leaves the system only at the cavity’s leakage rate γ cav. Therefore, fast radiative decay of the plasmon array mode becomes irrelevant, and the polariton linewidths are mostly determined by the total cavity’s and non-radiative array’s decay rates. This qualitative argument explains the apparent narrowness of the observed polaritonic bands. For a more rigorous description of the polaritons linewidth in the ultrastrongly coupled system (which is outside the scope of this work), the master equation approach might be needed40.Ground-state energy and photonic occupancy . Having performed the fitting of the experimental data, we can analyze how the ground state of the system \(\left. {G} \right\rangle\)is modified by the ultrastrong coupling. Again, we will restrict ourselves to the single-mode model of the system. In the uncoupled case, the global ground state is a direct product of the zero-photon and zero-plasmon states \({ { {G} \rangle = 0_{cav}} \rangle \otimes 0_{pl}} \rangle\), and the energy of this state is \(E_G = \langle {0H_{cav} + H_{pl}0} \rangle = \frac{{\hbar} }{2}(\omega _{cav} + \omega _{pl})\), correspondingly. The USC modifies the global ground state \(\left. {\tilde G} \right\rangle\)by admixing the states with different number of excitations, i.e., the global ground state with the higher excited states15, thus modifying the ground-state energy. Since after diagonalization, the coupled system comprises two new harmonic oscillators, its ground-state energy is \(\tilde E_G = \frac{\hbar }{2}(\omega _ + + \omega _ - )\).The A2 term modifies not only the energies of the polaritonic excited states41, but also the ground-state energy of the system. Since the ground-state energy of a harmonic oscillator (or a set thereof) is half the transition energy (sum of those), its modification can be calculated as \(\delta E_G = \tilde E_G - E_G = \frac{\hbar }{2} (\omega _ + + \omega _ - - \omega _{cav} - \omega _{pl})\). By expanding the solution of Eq. ( 6 ) (see “Methods”) into a Taylor series near g C?=?0, the ground-state energy modification can be approximated by: $$\delta E_G = \, \frac{{\omega _{cav}}}{{\sqrt 2 \omega _{pl}}}\frac{{\sqrt {A+B} - \sqrt {A-B} }}{{B}} g_C^2 + O\left( {g_C^4} \right),$$ (3)where we notated? \(A=\omega_{cav}^2+\omega_{pl}^2\)and \(B=\omega_{cav}^2-\omega_{pl}^2\).? At?zero detuning ( ω pl?=? ω cav) this expression?yields \(\delta E_G = \frac{{g_C^2}}{{2\omega _{cav}}} + O\left( {g_C^4} \right)\).The ground-state energy change at zero cavity–plasmon detuning can be estimated as \(\delta E_G \approx \frac{{g_C^2}}{{2\omega _{cav}}}\), which for g C/ ω pl?≈?0.5 becomes roughly δE G?≈? g C/4?≈?75?meV accounting for about 12% of the unperturbed ground-state energy E G. This value of the relative energy modification is smaller than what could be obtained in the system studied in ref.25with \(\frac{{g_C}}{{\omega _0}} \approx 0.73\). However, our absolute value is much greater because our system exhibits interacting resonances in the near-IR to visible range with resonant energies around 1?eV, while the characteristic energies of the system in ref.25lie 1 order of magnitude lower at around 100?meV (additionally the plasmonic nanoparticle array in our case fills only about 4% of the cavity interior, as opposed to ref.25where the active material fully saturates the mode volume). Thus, the absolute ground-state energy change in our system is several times greater than k B T at room temperature. This, in turn, implies that such ground-state energy modification might be important in practice and may show up in realistic USC-related effects even at room temperature.The normalized ground-state energy variation \(\frac{{\delta E_G}}{{E_G}} = \frac{{\tilde E_G - E_G}}{{E_G}} = \frac{{\omega _ + + \omega _ - }}{{\omega _{cav} + \omega _{pl}}} - 1\)calculated using the obtained coupling strengths and analytical expressions for polariton energies ω± , Fig.? 5a , predicts up to ?10% modification of the ground-state energy for normal incidence Fabry–Pérot mode upon coupling with the plasmonic array (see Supplementary Fig. 18 for the ground-state energy modification for other nanorod lengths). We want to emphasize that this vacuum energy is not an arbitrary reference level for all the higher energy states of the system. If the coupling constant can freely vary, for example, by allowing a single particle to move across a landscape with varying coupling strength, it will come to a state with the lowest vacuum energy even if the system is not coherently or thermally excited.Fig. 5: Modification of the vacuum state by the ultrastrong coupling.a The normalized vacuum energy variation in the coupled plasmon–cavity system as a function of the bare cavity energy for normal incidence eigenmodes calculated with the coupling strength obtained from fitting of the Lrod ?=?300?nm system, as well as the vacuum energy variation calculated directly from the measured polariton energies (error bars show 95% confidence intervals for indirect measurements). b The photonic occupancy of the modified ground state in the coupled system as a function of the bare cavity energy calculated with the coupling strength obtained from fitting of the Lrod ?=?300?nm system. Full size image These theoretical values follow the experimentally obtained trend (circles), which was obtained using the measured cavity and polariton energies with only the bare plasmon frequency ω pladopted from the fitting. The theory predicts a relatively slow dependence of the normalized ground-state energy change on the detuning, whereas the experiment is more sensitive to that. This can be explained by the error in the determination of polaritons energies: in particular, the UP energy has been extracted with the error of up to ±0.1?eV, which already constitutes a few percent of the total ground state energy. Calculating the difference of two close values \(\tilde E_G - E_G\)makes the relative error even worse. An additional possible source of disagreement is that the true polariton energies do not exactly correspond to the extrema of a response function, such as transmission or reflection, but lie rather close to them due to the multi-mode nature of the system and the Fano resonance mechanism6. Despite the non-ideal agreement, however, we stress that both theoretical predictions and experimental reflectivity data signal the ground-state energy modification of the order of 10% in our plasmon–microcavity systems. Such a modification is a clear hallmark of ultrastrong coupling, since in the conventional strong coupling picture, where g? ω , the additive coupled and uncoupled energies are exactly the same, i.e., ω+ ?+? ω? ?=? ω cav?+? ω pl, as can be seen from the Jaynes–Cummings model.Lastly, we study the photonic occupancy \(\tilde n_{\mathrm{phot}} = \left\langle {\tilde G\hat a^\dagger \hat a\tilde G} \right\rangle\)of the modified ground state \(\left. {\tilde G} \right\rangle\)(the plasmon occupancy of the ground state \(\left\langle {\tilde G\hat b^\dagger \hat b\tilde G} \right\rangle\)equals the photonic one15). In the USC regime, the ground state of the system acquires a non-zero photonic component due to the aforementioned admixing of states with different excitation numbers15. The photonic occupancy calculated with the use of the extracted coupling strength for the Lrod ?=?300?nm coupled systems, shown in Fig.? 5b , suggests that the ground state of the ultrastrongly coupled system may contain up to 0.06 bare cavity photons for cavities resonant with the nanorod array ( ω pl?=? ω cav ??0.5?eV). This is smaller than 0.37 photons estimated for Landau polaritons in the THz range21, but it is still a feasible number for converting to real photons by fast modulation of the coupling strength. The photonic occupancies calculated for other plasmonic nanorods predict almost identical values, Supplementary Fig.? 19 . .