Electrochemically Switchable Multimode Strong Coupling in Plasmonic Nanocavities

The strong-coupling interaction between quantum emitters and cavities provides the archetypical platform for fundamental quantum electrodynamics. Here we show that methylene blue (MB) molecules interact coherently with subwavelength plasmonic nanocavity modes at room temperature. Experimental results show that the strong coupling can be switched on and off reversibly when MB molecules undergo redox reactions which transform them to leuco-methylene blue molecules. In simulations we demonstrate the strong coupling between the second excited plasmonic cavity mode and resonant emitters. However, we also show that other detuned modes simultaneously couple efficiently to the molecular transitions, creating unusual cascades of mode spectral shifts and polariton formation. This is possible due to the relatively large plasmonic particle size resulting in reduced mode splittings. The results open significant potential for device applications utilizing active control of strong coupling.


Sample design:
Our system is the nanoparticle-on-mirror (NPoM) with an ultrasmall gap of ∼ 1 nm and designable properties, filled with a monolayer of spacer cucurbit[n]uril (CB[n]) molecules capable of hosting active quantum emitters in the form of methylene blue (MB) molecules in the nanogap.These macrocyclic CB[n] molecules help to preclude aggregation of the MB molecules and to assemble them with deterministic orientation along the z-direction in the cavity, via supramolecular guest-host chemistry.In particular, CB [7] can accommodate exactly one MB molecule (MB:CB [7] = 1:1) within its hollow, hydrophobic interior, while carbonyl portals at either end of the molecules bind them with their rims flat onto the gold surface.When a monolayer of CB [7] is first deposited on the gold mirror and suitably filled with MB molecules, gold nanoparticles can then bind on top to form the filled nanocavities, with the MB molecules vertically aligned within the nanogap.Previous spectroscopic measurements with empty CB[n] molecules have shown that such a particle-mirror gap is approximately 0.9 nm in extent and presents an effective refractive index of 1.4. 1

Sample preparation:
1 mM CB [7] and 1 mM MB are mixed in aqueous solution, to ensure the ratio of MB:CB [7]   is 1:1. 2 The Au substrate is immersed in this solution, rinsed with water, and dried with nitrogen.The 40 nm citrate capped AuNPs (BBI Solutions) are drop casted onto the Au substrate for 10 s.After this, the substrate is rinsed and dried.

Spectro-electrochemical cells:
The cells are made by 3D printed chambers with adhering glass coverslips using double sided tapes.The supporting electrolyte used is phosphate buffer solution (PBS) composed of 0.1 M monobasic and dibasic potassium phosphates.The NPoM samples are immersed in the electrolyte, and placed under the coverslips.The Au substrate is used as the working electrode.A Ag/AgCl reference electrode (3 M KCl, eDAQ ET072, Green Leaf Scientific) and a Pt mesh counter electrode (Alfa Aeser) are used.A potentiostat (Autolab PGSTAT204,

Optical measurements:
Photoluminescence spectroscopy is recorded on a modified Olympus BX51 with a 633 nm laser at the power of 100 µW.The spectra are recorded with an Andor camera coupled to a Triax 320 spectrometer.
Dark-field scattering spectroscopy is performed on a modified Olympus BX51 with a halogen lamp.The spectra (lower case 's') are recorded using a fibre-coupled OceanOptics QE65000.
Automated tracking microscopy was performed using a Python particle-tracking code. 1 A standard diffuser (Labsphere) is used as a reference to normalise white light scattering.
For both spectroscopy methods, excitation and collection happen through a 0.8-NA Olympus LWD BF/DF objective.
The 50 recorded scattering spectra in Fig. S1 without MB (a, uncoupled) and with the blueemitting MB molecules in (b) evidence a clear difference with and without strong coupling between the MB molecules and NPoM modes.In addition, slight variations from nanoparticle to nanoparticle show that all cavities are slightly different due the particle sizes and facets, which makes an exact reproduction in simulation impossible.

Cyclic voltammetry (CV):
The reduction potential of MB hosted in CB within the NPoM geometry is measured by a typical CV measurement.As MB reduction at a Au electrode completes at −0.5 V, 3 the potential window is chosen from −0.5 to +0.5 V.The recorded CV in Fig. S2 shows reduction of MB completes at −0.5 V, and this value is used for chronoamperometry in Fig. 1 in the main text.

NPoM density:
In Fig. S3 we present an exemplary image of the collection area (a) and NPoMs (b).We find less dense areas where a single NPoM response can be collected.

FDTD method:
Our full-dimensionality finite-difference time-domain (FDTD) simulations 4 of the NPoM nanocavity are performed using Ansys Lumerical software (2022 R1.1 release).The gold (Au) nanoparticle (depicted in Fig. 1a of the main text) is modelled as a truncated sphere of diameter D = 80 nm, with a facet width of w = 20 nm.The CB [7] layer is treated as a uniform dielectric medium with refractive index 1.4 1 and thickness d = 1 nm, separating the Au nanoparticle from an infinite Au substrate.For the permittivity of Au, Au , we employ the data of Johnson and Christy, 5 but to account for additional dissipation mechanisms stemming from inhomogeneities or imperfections of the NPoM constructs in the experiment, we enhance the damping through multiplication of the imaginary part of Au by 1.2 (as in Fig. 2 and Fig. 4).The ambient dielectric environment of the NPoM is assumed to be water with a refractive index of n = 1.33.Our FDTD calculations employ a 1 nm mesh size in a region of dimensions (100 × 100 × 100) nm 3 containing the NPoM structure, which is progressively refined to a mesh size of 0.25 nm in the nanogap region, ensuring an accurate description of the spatial variations in the plasmonic near-field therein.The injection of an exciting plane wave into the simulation domain (with angle of incidence θ i = 90 • ) and the calculation of the scattering cross-section in a spectral range from 300 nm to 1000 nm are performed using the total-field scattered-field (TFSF) technique; 4 in particular, the totalfield region is defined by a cuboidal volume of dimensions (200 × 200 × 180) nm 3 containing the NPoM and the scattering cross-section is proportional to the total power transmission across the TFSF boundary.The entire simulation domain is (1 × 1 × 1) µm 3 , and open boundaries are described by means of perfectly matched layers 4 using standard profiles.

Two-level Maxwell-Bloch model:
In order to study the interaction of the NPoM nanocavity modes with an active layer of MB emitters, we employ a self-consistent, two-level Maxwell-Bloch description 6,7 with Hamiltonian Here, Ĥe = hω e σ † σ is the field-free Hamiltonian of the two-level molecular emitter, in which ω e denotes the transition frequency, while σ † and σ are the raising and lowering operators, respectively.The operator V (t) = − μ • E(t) describes the dipolar interaction between an emitter and the electromagnetic field, in which μ = µ(σ † + σ) is the dipole moment operator and µ is a dipole matrix element.For our simulations, we assume that the dipole moment is oriented in the z-direction, so that µ = µe z .
We regard each two-level emitter as an open quantum system which, in addition to its dipole interaction with an external field, experiences relaxation and dephasing.The temporal evolution of the emitter density operator ρ under these conditions can be described by means of the following quantum master equation in Lindblad form, where σz = σ † σ − σσ † , while γ r and γ d are incoherent relaxation and pure dephasing rates, respectively.Forming matrix elements of the left-and right-hand sides of Eq. ( 2) in a basis S7 comprising the emitter ground and excited states, we find the following coupled system in which we have introduced the total dephasing rate Γ = γ d + γ r /2.
For an ensemble of emitters distributed with density N d , we define the ground-and excited-state population densities as N 1 = N d ρ 11 and N 2 = N d ρ 22 respectively, and calculate the macroscopic polarization via By separating Eq. ( 4) according to the real and imaginary parts of the off-diagonal density matrix elements, and using Eq. ( 5), we can rewrite the coupled system of equations in the form In this work, we solve numerically MaxwellâĂŹs equations with the optical Bloch equations [Eqs.( 6) and ( 7)] in a completely self-consistent fashion, where the local electric field (computed using the FDTD method) drives a polarization response from the quantum emitters (determined by solving the optical Bloch equations), which in turn couples back to the field.The entire CB [7]:MB layer is regarded as the active medium, with a nominal, single-

Metal-insulator-metal waveguide model:
To develop a minimal model that captures the evolution of the polariton pairs for different cavity modes, we invoke the metal-insulator-metal (MIM) waveguide analysis from Ref. 8,9 Specifically, we treat the NPoM nanogap as a one-dimensional waveguide structure comprising an insulator of thickness d and dielectric permittivity g = n 2 g , confined between two metal plates with common permittivity m .The dispersion of the MIM modes in such a structure is given by 10 where k 0 = 2π/λ is the vacuum wavenumber and γ = {(2 g )/(k 0 d m (ω)]} 2 .We use this effective refractive index to match the Fabry-Pérot boundary condition for transverse magnetic resonances in a finite-size cylindrical gap of height d, which leads to the problem of self-consistently solving the equation 11 in which α mn is the nth root of the mth order Bessel function J m and β mn is a modedependent phase shift.The latter can be used to account for differences in the frequencies of the confined modes stemming from the more complex geometry of our NPoM system. 12,13 mimic the resonant behaviour of the active material in the nanogap (i.e., the spectrally sharp emitters), we adopt the following model permittivity, where B = 1.4 2 is the background dielectric constant and f is an overall coupling factor, while ω 0 and Γ are the energy and linewidth of the two-level medium, respectively.
This method is of course inadequate to match the exact spectral positions of the polari- 500 600 700 800 900 1000  9) and the dashed horizontal lines the constants α mn − β mn .(b) Polariton wavelengths for multimode strong coupling, according to the MIM model, as the coupling constant f increases from 0 to 6.The black dashed line indicates the emitter transition wavelength, while the colored lines show how the polariton positions evolve with increasing strength of interaction between the emitters and the plasmonic modes l 1 (purple), l 2 (yellow) and l 3 (green).

S10
tons, and the coupling strengths between the various plasmonic modes and the emitters are, in general, different.As such, using only a single coupling constant f for all plasmonic modes is insufficient, but it nevertheless suffices to explore qualitative trends in multimode strong coupling, as demonstrated in Fig. S5.The quantities α 0n − β 0n for the plasmonic modes l n , where n = 1, 2, and 3, are shown by the dashed lines in Fig. S5a, and the black/grey curves show the dispersion πdn eff (λnm) λmn for 0 ≤ f ≤ 3.As mentioned before, the constant β 0n is introduced in the literature to account for mode-dependent boundary behaviour and is used here to adjust the mode frequencies according to our FDTD simulation results.The three points of intersection between the dispersion function with f = 0 and the dashed lines [i.e., the solutions of Eq. ( 9)] determine the spectral positions of the three plasmonic modes l n in our NPoM system, while they correspond to the polariton peak locations for non-zero f .
As shown in Fig. S5b, with increasing coupling constant f the points of intersection to the left of λ 0 = 2πc/ω 0 are blue-shifting, while those to the right are red-shifting.If the emitter transition is assumed spectrally sharp without any damping, we obtain a unique solution λ 0 = 650 nm, which is the nominal transition wavelength.However, for non-vanishing broadening as considered here, we take the outermost points of intersection to mark the polariton peak locations, designated n + and n − .
The polariton spectral positions are extracted for coupling constants in the range 0 ≤ f ≤ 6, as shown in Fig. S5b.The colored solid lines indicate how the polariton wavelengths change with increasing strength of interaction between the emitters and each cavity mode l 1 (purple), l 2 (yellow) and l 3 (green).The splitting for each pair of polaritons becomes approximately linear and is qualitatively similar to the trends discussed in connection with Fig. 3 of the main text.
Classical coupled-oscillators model: Simplifying our analysis one step further and neglecting the multimode character of the coupling, we can employ a classical, coupled-oscillators model for each pair of emitter transition and NPoM cavity mode.In this elementary model, we consider a 2 × 2 Hamiltonian matrix to describe the coupling between a plasmonic mode l n and emitters (ω 0 ) with potentially different dissipation rates κ n , Γ and frequencies ω Since the coupling strength is proportional to the dipole moment, we can write g n = g µ is the coupling strength per unit dipole moment and assumed to be a constant for each plasmonic mode and a fixed geometry.By diagonalisation, we obtain two eigenvectors of this Hamiltonian corresponding to a pair of polaritons, For the numerical simulations, we set ω 0 = 1.91 eV and Γ = 6.6 meV.In order to obtain the spectral position ω a and κ n for plasmonic mode l n , we employ a classical dipole source positioned at the center of the NPoM nanogap (as in Fig. 3), attributing the radiative decay to the dominantly radiative l 1 and l 2 cavity modes.By fitting with Lorentzian line-shapes, our peak positions for l 1 and l 2 are found to be ω We fit the polariton peak positions from our full numerical simulations in Fig. 3a of the main text using the functions defined by Eq. ( 12), where the coupling strength per unit dipole moment g µ between the emitters and each plasmonic mode l n is the quantity to be estimated.Note that higher energy modes usually have a broadening effect and overlap with each other, which renders it difficult to extract the peak positions accurately.Therefore, we focus on the coupling between the emitters and the l 1 , l 2 modes for our coupled-oscillators treatment.Higher-order hybrid states: For smaller magnitudes µ of the dipole moment, one can assume that the quantum emitters couple to multiple plasmonic modes simultaneously and independently with different effective coupling strengths, as we discuss in the main text.As the polariton resonances shift with different slopes when increasing µ, a situation can arise in which different polaritons approach the same energy as µ becomes very large.We study this situation in Fig. S7, where we plot the scattering cross-section on a logarithmic scale for better visibility.We clearly see that the polaritons 2 + and 1 + approach 3 + and form anti-crossings with it around µ = 15 D and 20 D, respectively, as marked by the dashed lines.Our numerical results thus point to the formation of higher-order hybrid states from these polaritons.

Figure 1 FigureFigure S2 :
FigureS4shows an absorption spectrum of an ensemble of MB molecules with a significant

FigureFigure S4 :
Figure S3: (a) Collection area identified by illumination along the collection path, circled in red for clarity.(b) To ensure only one NPoM is in the collection area (circled white), regions of sparsely distributed NPoMs are selected.NP diameter = 80 nm.100x, 0.8 NA objective.
emitter transition wavelength of λ e = 2πc/ω e = 650 nm (where c is the speed of light in vacuum), a dipole moment of µ = 4 D, a total dephasing rate Γ = 5 × 10 13 rad/s and an emitter density of N d = 1.5 × 10 26 /m 3 , corresponding to N ≈ 50 MB molecules under the nanoparticle facet.

Figure
Figure S5: (a) Examples of how the polariton wavelengths are determined.The black/grey curves represent the dispersion in Eq. (9) and the dashed horizontal lines the constants α mn − β mn .(b) Polariton wavelengths for multimode strong coupling, according to the MIM model, as the coupling constant f increases from 0 to 6.The black dashed line indicates the emitter transition wavelength, while the colored lines show how the polariton positions evolve with increasing strength of interaction between the emitters and the plasmonic modes l 1 (purple), l 2 (yellow) and l 3 (green).
Fig. S6 presents the fits, where black dots indicate the energy of the uncoupled emitters, and the colored dots show the energies of the polaritons extracted from the FDTD simulations, with purple for l 1 and yellow for l 2 .The colored solid lines represent fits of moment g (1) µ = 16.1 meV/D and g (2) µ = 16.3 meV/D.By comparing with the linewidth, we find that the strong coupling regime for l 1 and l 2 , where g n > (κ n + Γ)/2, is achieved for µ > 3.4 D and µ > 3 D, respectively.

Figure S6 :
Figure S6: Fitting of the numerical simulation results (colored dots) according to a classical oscillator model (colored solid lines), for coupling between the emitters and (a) the l 1 and (b) the l 2 modes.The black dots indicate the uncoupled emitter energy.

Figure S7 :
Figure S7: Simulated scattering spectra with varying magnitude of the dipole moment from µ = 0 D to 42 D. The dashed white lines act as a guide to the eye for the polariton resonances.