Plasmon mediated coherent population oscillations in molecular aggregates

The strong coherent coupling of quantum emitters to vacuum fluctuations of the light field offers opportunities for manipulating the optical and transport properties of nanomaterials, with potential applications ranging from ultrasensitive all-optical switching to creating polariton condensates. Often, ubiquitous decoherence processes at ambient conditions limit these couplings to such short time scales that the quantum dynamics of the interacting system remains elusive. Prominent examples are strongly coupled exciton-plasmon systems, which, so far, have mostly been investigated by linear optical spectroscopy. Here, we use ultrafast two-dimensional electronic spectroscopy to probe the quantum dynamics of J-aggregate excitons collectively coupled to the spatially structured plasmonic fields of a gold nanoslit array. We observe rich coherent Rabi oscillation dynamics reflecting a plasmon-driven coherent exciton population transfer over mesoscopic distances at room temperature. This opens up new opportunities to manipulate the coherent transport of matter excitations by coupling to vacuum fields.


Introduction
The strong couplings of quantum emitters to light 1,2 emerges as a critical instrument for directing the optical 3,4 and electronic transport [5][6][7][8] properties of nanomaterials by all-optical means 3 .These phenomena can be exploited to modify the outcome of photochemical reactions in the electronic ground 9 and excited states 7,[10][11][12] , for creating new and unusual states of condensed matter systems such as polariton condensates 13 at room temperature 14 or for designing entirely new optical device concepts such as polaritonic lasers 15 operating at room temperature 16 .In particular the strong coupling of quantum emitters to surface plasmon excitations in metallic nanostructures [17][18][19][20] has caught much attention since the strong nanometer-scale spatial confinement of the plasmonic mode [21][22][23][24] offers a direct path for locally enhancing vacuum field fluctuations 25 and -thus -the coupling to quantum matter.Remarkably, it has led to the demonstration of strong coupling to a single molecule 26 or quantum dot 19,27 , fundamental hallmarks of quantum plasmonics 28 .Optical excitations in metallic nanostructures are inherently short lived due to substantial Ohmic and radiative losses 22 .Consequently, quantum coherence in metal-based hybrid systems is usually lost after a few tens to hundreds of femtoseconds 21 .This has, so far, limited experimental work to mostly linear optical studies, either on ensembles or single quantum structures.Advanced twodimensional coherent electronic spectroscopies (2DES) 29 , have emerged as powerful tools for probing quantum-coherent couplings in the time domain [30][31][32] even in strongly dephasing media 33,34 and particularly for accessing many-body excitations [35][36][37][38] in quantum systems.However, these techniques are yet to be applied to hybrid plasmonic systems in the strong coupling regime.In 2DES, the excitation with a pair of phase-locked, short optical pulses allows for the selective excitation of different resonances in the system.A third, time-delayed probe pulse generates twodimensional energy-energy maps of the optical response that correlate the optically excited and detected resonances.Distinct cross peaks in these maps, oscillating as a function of the time delay between excitation pair and optical probe are unambiguous signatures of strong, quantum-coherent couplings in the hybrid systems.Here, we demonstrate such temporally oscillating cross peaks in the 2DES spectra of a model molecular aggregate system revealing collective strong coupling between excitons and surface plasmon polaritons.Our results show that these oscillations are dominated by a robust coherent population transfer between spatially separated strongly and weakly coupled excitons in different regions of the nanostructure, mediated by the plasmonic field.This suggests light-driven coherent transport of matter excitations in nanosystems over mesoscopic distances at ambient conditions.

Results and Discussion
To probe strong couplings between excitons and plasmons, we performed angle-resolved 2DES on a hybrid plasmonic cavity, a gold nanoslit array covered with a J-aggregated thin film (Fig. 1a).This molecular aggregate is based on squaraine monomers (ProSQ-C16, inset in Fig. 1b).Their electronic properties are reasonably well described within phenomenological essential state models 39 , showing that only the lowest-lying electronically excited state is relevant for the present work.When deposited on gold, they form well-ordered J-aggregated thin films 40,41 .Dipolar coupling among neighboring molecules results in a delocalization of the optical excitation across ~ 20 -30 monomers at room temperature and in the formation of strongly red-shifted and spectrally narrow superradiant exciton |⟩ resonances at around 1.59 eV 41 (Fig. 1b).In 2DES, this results in a spectrally sharp and well isolated exciton peak with a dispersive line shape along the detection axis (Fig. 1c).Such 2DES maps are recorded by exciting the sample with a pair of collinear pump pulses with time delay  and by probing the pump-induced change in sample reflectivity with a broadband probe pulse that is delayed by a waiting time  with respect to the second pump.Fourier transform of the resulting differential reflectivity spectra recorded for various , at fixed , results in the 2DES map in Fig. 1c (see Methods).In the following, we label the peaks in these maps as (Ex,Det), where Ex and Det denote the excited and detected resonances, respectively.The dispersive line shape arises from a superposition of three contributions: ground state bleaching (GSB) and stimulated emission (SE) of the one-quantum |⟩ resonance overlap with an excited state absorption (ESA) from |⟩ to two-exciton, |⟩ , states.In J-aggregates, the two-quantum |⟩ resonances are blue-shifted (Δ) since Pauli-blocking in each monomer dictates that the lowest-lying delocalized exciton state in every aggregate can be populated only once 42 .Slight line broadenings result from higher-lying aggregated exciton states.Exciton relaxation within the disordered aggregates leads to a partial decay of the |⟩ peak on a 100-fs time scale (Fig. 1d).We deposit such J-aggregated thin films, with a thickness of 10 nm, onto a plasmonic nanoslit array, milled into a 200-nm thick gold film.The nanoslits form a cavity that locally confines optical near fields and strongly enhances their coupling to the in-plane component of the exciton dipole moment (Fig. S11).Width, height (45 nm) and period (530 nm) of the array are chosen to create sharp surface plasmon polariton (SPP) resonances of the grating with an energy that can be tuned across the exciton resonance by varying the incidence angle  (Fig. S1a).Angle-resolved linear reflectivity spectra (Fig. 1e) show that the collective dipolar coupling between excitons and SPPs results in the formation of mixed upper (UP) and lower (LP) polariton branches.From the avoided crossing, we deduce a normal mode splitting of ~60 meV, twice the Rabi energy ℏΩ  .The polariton branches are superimposed by an angle-independent exciton peak.It is commonly thought to arise from "uncoupled" excitons which are only weakly interacting with the SPP field, e.g., because they lie in regions outside the slits with much reduced local field enhancement 21,43,44 .This interpretation of the linear spectra is well supported by finite difference time domain (FDTD) simulations of Maxwell's equations (Fig. 1f and S1b).Most of these resonances also appear in angle-dependent 2DES spectra recorded at  = 0 fs (Fig. 2a-c and Fig. S3).Along their diagonal, we observe strong (LP,LP) and (X,X) peaks with dispersive line shapes along   .In contrast, the (UP,UP) peak is much weaker and appears only at angles below the crossing at   = 23°.In addition, we find pronounced cross peaks between LP and X, both below and above the diagonal.Their dispersive line shapes are best seen for  = 27° in Fig. 2c.Cross peaks between UP and both, LP and X, are much weaker in amplitude and are only resolved for  <   .While (X,UP) and (UP,X) have dispersive line shapes, the other weak UP peaks appear absorptive in shape.Resonance energies and spectral line shapes deduced from 2DES are supported by angle-resolved pump-probe measurements, shown in Fig. S5.The quite pronounced cross peaks between "uncoupled" excitons and polaritons are unexpected since the uncoupled excitons are thought to be spatially well separated and thus essentially uncorrelated with those excitons that are hybridized with the SPP field.Hence, it is not obvious that their excitation should not result in a polariton nonlinearity.The observation of dispersive line shapes for both diagonal and cross peaks now allows us to correlate one-quantum (1Q) resonances, characterized by a positive GSB and SE peak, and twoquantum (2Q) resonances with a negative ESA signal 41,45 .We deduce 1Q energies from peak maxima along   , while 2Q energies are taken as the zero crossing of a dispersive peak along   .The resulting energies are plotted in Fig. 2d as open circles, together with the 1Q dispersion (solid lines) deduced from angle-resolved reflectivity.In addition, the 2Q dispersions are estimated by adding the energies of the contributing 1Q states without further corrections.Obviously, the 1Q energies obtained from 2DES match those deduced from linear spectroscopy, while the 2Q dispersions show several new features.Since the experiment probes the collective coupling of many excitons to a single plasmonic mode, we expect, from a commonly employed Tavis Cummings (TC) model, to observe three distinct 2Q states 35,45 .The model predicts doubly excited 2LP and 2UP polaritons and a mixed UP/LP state, while all other states remain optically inactive ("dark") 45 .Indeed, these resonances have been seen in the 2Q dispersions measured for semiconductor microcavities 35 and a TC model has recently also been used to discuss organic microcavity polaritons 45 .As a result of the fermionic nonlinearity induced by the exciton part of the wavefunction, the energies of the 2Q states are slightly blue-shifted with respect to twice the 1Q transition.This blue shift is proportional to the two-exciton fraction of the 2Q wavefunction 35,45 .Since the doubly excited X state is uncoupled from the plasmon branch, an angle-independent XX contribution is expected(Fig.2).In the experiments, the 2LP resonances are clearly resolved, while 2UP and UP/LP are apparently lacking.The mixed resonances (grey and black circles) follow the UP/X and LP/X dispersions with a distinct avoided crossing at   .The appearance of those resonances goes beyond the TC model and requires further discussion.Before that, however, we inspect the waiting time dynamics of the 2DES spectra.This is exemplarily done in Fig. 3a for  = 27°.Indeed, we observe pronounced coherent oscillations of the amplitude on all diagonal and cross peaks in the 2DES map, except for the (X,X) peak.Their oscillation period   = 2/(  −   ) matches the splitting between the X and LP resonances.Coherent oscillations with the same period are also observed in pump-probe spectra (Fig. 3e, Fig. S5).The traditional understanding of strong X-SPP coupling would, instead, predict Rabi oscillations with a period   = 2/(  −   ) given by the normal mode splitting of the polaritons 1,21,22,46 .Angle-dependent pump-probe transients, detected at the LP resonance (Fig. 3e and Fig. S5), again reveal   oscillations with a period that decreases monotonically with increasing angle.For  ≥ 21°, also the 2DES maps show persistent amplitude oscillations, except for (X,X) (Fig. S4).These oscillations appear predominantly with a period given by the X-LP splitting, as can be seen by comparing the measured oscillation periods (blue squares in Fig. 3f) to those predicted by the linear dispersion (blue line).Only at one selected angle, close to the crossing, we find an oscillation at the anticipated UP-LP splitting (red circle in Fig. 3f and Fig. S18).To explain these observations, we first introduce a phenomenological extension of the TC model that takes the spatial characteristics of our sample into account.For this, we consider two classes of spatially separated J-aggregated excitons.Excitons in the slit region,   , collectively interact with the plasmon field with a coupling strength   .In contrast, those excitons,   , that lie between the slits, on the flat gold film, interact with   .Both plasmons and excitons are treated as bosonic oscillators 35 .A nonlinearity of the system arises by introducing a finite blue-shift Δ of both twoexciton states (Fig. S14).Using this model, 2DES spectra based are simulated by solving the Lindblad master equation for the density matrix of the coupled system.As can be seen in Figs.3b,d, the model quantitatively accounts for our experimental observations.Specifically, the simulations show dispersive peaks with pronounced amplitude oscillations at the period T X , given by the X-LP splitting, while Rabi oscillations at   are much weaker in amplitude.As in the experiment, the oscillations are basically absent at (X,X).Reasonable agreement between experiment and simulation is achieved when choosing   ≃ 3  , with ℏΩ  = √  2 +   2 , and   ≃   ≃ 3  (see Section 9 of the Supplementary Information).
To rationalize the microscopic dynamical processes that give rise to these transient 2DES spectra we further invoke an elementary Frenkel exciton model 41,42 .We consider a disordered chain of squaraine molecules, each treated as a fermionic two-level system.Neighboring molecules are dipole-coupled via their optical near fields and interact with the plasmonic mode that is delocalized along the chain.In agreement with FDTD simulations of our sample (Figs.S6 -S12), we consider a spatially inhomogeneous plasmon field with a local field in the slit region that is ten times larger than in the region between the slits (Fig. 4a).In the absence of the plasmon field (Fig. S13), the nearest-neighbor coupling results in the formation of few superradiant, moderately localized Jaggregated excitons, strongly red-shifted in energy and extending over ~25 molecules, together with a large number of dark excitons.Between the slits, the wavefunctions of these localized excitons () remain basically unchanged in the presence of the coupling to the plasmon modeexcept for a minor admixture of excitons inside the slits.The plasmon contribution to their wavefunction is small.In contrast, the superradiant excitons inside the slits couple strongly to the plasmonic mode, resulting in a LP mode carrying substantial contributions from   and  and much weaker contributions from all   .For LP, all wavefunctions interfere constructively while for the  states, the contributions from   and   interfere destructively.The resulting linear optical absorption (Fig. 4b) shows strong contributions from the energetically isolated LP state, while the X peak is inhomogeneously broadened.The UP absorption is much weaker since, in our sample, the dipole moment of  (µP) and of the sum of all excitons (µW and µS) are of similar magnitude.Hence, their emission interferes destructively for the UP peak.We now discuss the dynamics of the coupled X-SPP system.For this, we impulsively excite all optical resonances with a spatially homogeneous laser field and follow the spatiotemporal evolution of the excited state populations within the chain of squaraine monomers and in the plasmon mode.The plasmon mode shows the expected Rabi oscillations with   .Out-of-phase oscillations at   are most pronounced for the slit excitons   .They are superimposed, however, with slower oscillations of the   population with a period   .While these slower oscillations are completely absent in the plasmon dynamics, they reappear, phase-shifted by , for those excitons,   , that are localized between the slits.These two distinct types of population oscillations are most clearly seen when spatially integrating over the localized exciton populations   and   (Fig. 4c).Now, the Rabi oscillations on  are perfectly matched by out-of-phase oscillations of the total population of all excitons,   +   (red line in Fig. 4c).The plasmon-mediated coherent population oscillations (CPO) between   and   are only seen in the individual exciton subsystems, while they are absent in the net exciton population.These model calculations suggest that the dipolar coupling to the plasmon induces spatial oscillations in the exciton density, from outside the slits into the slit region and back.These oscillations appear at the period   , given by the energy splitting between X and LP.They clearly dominate the exciton dynamics in the region between the slits.Here, the effect of the "traditional" Rabi oscillations with period   is weak due to the small local plasmon field amplitude.These conclusions are largely corroborated by Fourier transforms of the population dynamics (Fig. 4d).They emphasize the presence of fast Rabi oscillations with   and absence of   oscillations in the dynamics of the plasmons and of the total exciton population,   +   , (blue and red line), respectively.In contrast, the slower oscillations with   between the two distinct classes of excitons become apparent when examining the individual exciton dynamics.These Frenkel exciton simulations form a convincing microscopic basis for the phenomenological extension of the TC model introduced above.Essentially, we can explain the suppression of "traditional" excitonplasmon Rabi oscillations (  ) and the emergence of CPOs with a longer period   in the waiting time dynamics of the 2DES maps by considering two spatially distinct classes of localized Jaggregated excitons that are mutually coupled to a spatially structured plasmonic mode.This model convincingly accounts for the rich spectral and dynamic features in all angle-dependent 2DES maps.Most importantly, it explains how a spatially delocalized plasmon mode induces a coherent real-space energy transport between spatially separated exciton sites that persists during the coherence time of the strongly coupled system.This maps the complex real space dynamics of a nanostructured system of molecular excitons and plasmons onto an effective three-level system, in which two of the levels,   , the excitons inside the slits, and   , the excitons between the slits, are both coupled to a third state, the plasmon.The collective coupling strength of the excitons inside the slits is roughly three times larger than that for those outside the slits.This spatially structured coupling gives rise to CPOs between two of these states without affecting the dynamics of the third.Previously, such CPOs have been discussed in atomic and molecular three-level systems in the context of slow light generation 47 and light storage 48 .Here, we report CPOs in a prototypical allsolid 3-level-system at room temperature and demonstrate how they enable an efficient coherent transport of excitons over mesoscopic distances, from regions outside to inside the slits and back.Atomic and molecular 3-level and 4-level systems offer exciting resources for quantum state manipulation and information processing.The reduction in the speed of light by electromagnetically induced transparency 49 , the coherent trapping of population in optically dark states by stimulated Raman adiabatic passage 50 or lasing without inversion 51 are among the manifestations of the control of optical information that can be achieved.We therefore anticipate that the demonstration of related coherent phenomena in all-solid-state systems will open up new avenues towards optical information processing in strongly coupled exciton plasmon systems.Such systems offer new prospects for manipulating coherent quantum transport by light extending to mesoscopic length scales.As such they emerge as an important class of hybrid materials for exploring and exploiting these opportunities.

Sample preparation
A polycrystalline 200-nm gold film was deposited on a fused silica substrate and subsequently annealed 52,53 .Nanoslit arrays (80 x 150 µm²) with a grating period of 530 nm and a slit depth and width of 45 nm were fabricated using focused Ga ion beam milling (Helios NanoLab 600i, FEI).The nanostructured sample was coated with a 10-nm thick J-aggregate thinfilm by spin-coating a solution of squaraine molecules, dissolved in chloroform, following the procedure reported previously 41 .The squaraine molecules, abbreviated as ProSQ-C16, are (S,S)-enantiomers of 2,4-Bis[4-((S))-2-(hexadecyloxymethyl)-pyrrolidone-2,6-dihydroxyphenyl] and have been synthesized as described by Schulz et al. 40,54 .Angle-dependent linear optical characterization of the fabricated samples have been performed using a supercontinuum white light source (SC400-4, Fianium) that was focused onto the structured area under an angle  with a polarization perpendicular to the slit orientation (as depicted in Fig. 1a).The reflected light was measured with a fiber spectrometer (FLAME, OceanOptics) and normalized to the reflection from bare gold.

Experimental 2DES setup
Angle-dependent pump-probe and 2DES data were recorded with a home-built high-repetition rate setup which allows for rapid data acquisition, resulting in a high signal-to-noise ratio of the measured data within short measurement times.To this aim we employ a high-repetition rate laser system (Tangerine V2, Amplitude Systèms), delivering 260-fs pulses (full width at half maximum of the pulse intensity), centered around 1030 nm, at a repetition rate of 175 kHz.A fraction of the 60 µJ pulse energy is used to pump a home-built non-collinear optical parametric amplifier (NOPA), previously described in ref 41 and based on a design published in ref 55 .The NOPA outputs 650-900 nm pulses with a measured pulse duration of ~12 fs, as characterized using second harmonic frequency resolved optical gating (Fig. S2).The NOPA pulses are used in a home-built 2DES setup, also previously reported in ref 41 .A phase-stable and collinear pump-pulse pair with variable delay  (coherence time) is generated by an interferometer based on birefringent wedges (Translating Wedge-based Identical pulse eNcoding System, TWINS 56 ).The pump pulses are periodically switched on and off at 43.75 kHz by a mechanical chopper system (MC2000B, Thorlabs) equipped with a custom-made blade (500 slots).The vertically aligned and p-polarized pump and probe beams are focused onto the sample to a spot size of ~60x60 µm² under the same angle of incidence  and with a small angle mismatch in the orthogonal direction (see Fig. 1a).The sample is mounted on a rotation stage and the axis of rotation is tuned to coincide with the location of the nanostructured sample to allow for conveniently and accurately tuning  for both pump and probe in the angle-dependent 2DES and pump-probe experiments.The reflected probe beam is then collected and sent to a monochromator (Acton SP2150i, Princeton Instruments) with an attached fast line camera (Aviiva EM4, e2v) allowing to rapidly record probe spectra S at an acquisition rate of 87.5 kHz (i.e., at half the laser repetition rate) as a function of the detection energy E det .We thus record a differential reflectivity spectrum Δ/ from a set of 4 laser shots by taking the difference between spectra with (S on ) and without (S off ) pump Here, the waiting time T denotes the delay between the second pump and the probe pulse.This delay is controlled using a motorized linear translation stage (M126.DG1, Physik Instrumente).For the 2DES measurements, at each waiting time T, the coherence time is scanned and differential spectra are recorded on the fly, whereas for the pump-probe experiments, we fix the coherence time to τ = 0 fs and scan only the waiting time T.
To calculate absorptive 2DES maps 57 a Fourier transform of the differential reflectivity is performed along the coherence time to obtain the 2DES spectra as a function of the excitation energy axis E ex (ℏ Planck's reduced constant, Θ(τ) Heaviside step function).

Frenkel exciton simulations
We model the coupling of the molecular J-aggregate with the plasmon mode by performing qualitative microscopic simulations based on a disordered Frenkel exciton model.We choose N = 300 squaraine monomer states at E SQ = 1.923 eV with Gaussian disorder in site energy of σ = 15.4 meV 41 .We set the monomer transition dipole moment to µ SQ = 1/√N.
The plasmon mode is introduced as a single state |P⟩ with energy E P and transition dipole moment µ  .The Frenkel exciton Hamiltonian reads where V XP,n = − ⃗  ⋅  ⃗⃗ , denotes a local exciton plasmon coupling between the monomer state |n⟩ at site n and the plasmon field  ⃗⃗ , at the same site.We choose a dipolar coupling between monomer states of J = −166 meV, and limit this coupling to nearest neighbors by setting  , =  ,±1 .Periodic boundary conditions are applied.For V XP,n = 0, Eq. ( 3) thus describes a superradiant J-aggregate exciton state with E X = 1.59 eV determined by the nearest-neighbor coupling between the monomers (Fig. S13b) 41,42 .
We choose the local X-P couplings V XP,n to be proportional to the amplitude of the plasmon field at each monomer site n.We assume that the local coupling is governed by the in-plane components of the local SPP field.To account for the spatial profile of the SPP field the X-P coupling is divided into two regions: V S,n for the strong coupling inside the slits and V W,n for a weaker coupling in the region in-between two slits.V S,n and V W,n are taken as constant and real-valued parameters.The number of monomers in the two regions are N S and N W , respectively.
For rationalizing the results of this Frenkel exciton model, we introduce an effective Hamiltonian comprising four states |0⟩, |P⟩, |X S ⟩ and |X W ⟩, where |X S ⟩ and |X W ⟩ indicate strongly coupled (inside the slits) and weakly coupled excitons (in between two slits), respectively.The ground state energy of the system is set to zero.The exciton energies are approximated as E X S = E X W = E SQ − 2J.The collective transition dipole moments and coupling strengths can be calculated as µ X S = √N s µ SQ and µ X W = √N W µ SQ for the transition dipole moments, and V S = √N S V S,n and V W = √N W V W,n for the coupling strengths.
We simulate the dynamics following optical excitation by numerically integrating the master equation in Lindblad form 57,58 ρ ̂̇= − H ̂= H ̂S + H ̂int (t) describes the free evolution of the system via H ̂S and its light-matter interaction is governed by the time-dependent interaction Hamiltonian H ̂int (t) = −μ ̂E(t) (6) which accounts for optical excitation of the system by an external light field in dipole approximation.μ ̂ denotes the transition dipole moment operator of the system.We assume a short and sufficiently weak 5-fs pulse at 1.6 eV.Dephasing and relaxation processes are incorporated through appropriate Lindblad operators L ̂k 58,59 .For more details see Supporting Information section 8.

Simulation of nonlinear signals
To simulate the experimental pump-probe and 2DES maps, we compute the full dynamics of the density matrix ρ ̂ considering a series of interactions with up to three laser pulses.This allows us to calculate the sample polarization P(t) = Tr(μ ̂ρ ̂(t)) and the resulting optical spectra 39 .
For the numerical integration of the master equation in Lindblad form, Eq. ( 5), we use a nonperturbative approach 58,59 .The total electric field, entering the interaction Hamiltonian in Eq. ( 6), comprises of up to three laser pulses (pump 1, pump 2 and probe) with amplitude E 0,n , pulse duration Δt = 10 fs, frequency ω L = 1.6 eV/ℏ, phase ϕ n for phase cycling and a time shift t n ′ .We label the delay between the pump pulses as the coherence time τ = t pu2 ′ − t pu1 ′ and the delay between the second pump pulse and the probe pulse as the waiting time T = t pr ′ − t pu2 ′ .A detection time t of zero corresponds to the arrival time of the probe pulse.The detection time axis coincides with the one used for the numerical integration of Eq. ( 5).
In the simulations, we calculate the linear polarization P (1) (t) without pump, or a total polarization P tot (τ, T, t), including all nonlinear signals that arise from all possible interactions with the three pulses 39 -We then deduce linear and total susceptibilities as a function of the detection energy E det from these polarizations via a Fourier transform along the detection time t χ (1) (E det ) = 1 ε 0 ℱ(P (1) (t))/ℱ(E pr (t)) χ tot (τ, T, E det ) = 1 ε 0 ℱ(P tot (τ, T, t))/ℱ(E pr (t)) where ε 0 denotes the vacuum dielectric constant.
We employ a phase-cycling scheme in our simulations to isolate the desired linear and nonlinear contributions to our signal.For this, we calculate P tot (τ, T, t) for four different phase settings of ] for both pump pulses, with ϕ pr = 0 [60][61][62] .We then perform the average over these four settings to obtain the phase-cycled average χ PC tot (τ, T, E det ).This allows for extracting the contributions to the nonlinear signal corresponding to those that are measured in the partially collinear experimental geometry.We then obtain the nonlinear signal as χ nl (τ, T, E det ) = χ PC tot (τ, T, E det ) − χ (1) (E det ).
(11) Finally, energy-energy 2DES maps for a fixed waiting time T are obtained by taking the real part of the Fourier transform along the coherence time τ S 2D (E ex , T, E det ) = ℜ(∫ θ(τ)S 2D (τ, T, E det )e iE ex τ/ℏ dτ ∞ −∞ ), (12)   yielding the excitation energy axis E ex .Setting τ = 0 in the simulations allows for calculating pump-probe spectra as a function of E det and T.
The exciton-plasmon system is modeled based on the previously introduced effective Hamiltonian and contains, in addition to the ground state, not only the three one-quantum states |P⟩, |X S ⟩ and |X W ⟩, but also six two-quantum states |2P⟩, |P, X S ⟩, |P, X W ⟩, |XX S ⟩, |X S , X W ⟩ and |XX W ⟩. The plasmon and exciton subsystems are each treated as a harmonic oscillator.To account for the exciton nonlinearity, the two-exciton state is blue-shifted by Δ= 5 meV 41,42 .Dipolar interactions of the strongly (V S ) and weakly (V W ) coupled excitons with the plasmon field are introduced by adding two coupling Hamiltonians in rotating wave approximation, following 35 H ̂XW P = V W (b ̂P † b ̂XW + b ̂XW † b ̂P) with b ̂ † and b ̂ being the creation and annihilation operators.Radiative damping and pure dephasing phenomena are included using the Lindblad formalism.See Supporting Information section 9 for more details.

Figure 1 :
Figure 1: Strong coupling of a gold nanoslit array coated with a J-aggregated thin film.a) Experimental geometry.The nanoslit sample is illuminated with a phase-locked pair of pump pulses, separated by the coherence time  at incidence angle .The pumpinduced change in sample reflectivity is monitored by a probe pulse with the same incidence angle and time-delayed by the waiting time T. b) Linear reflectivity of a 10-nm J-aggregated film of ProSQ-C16 squaraine molecules (inset) on a gold surface.c) Experimental, reflective two-dimensional electronic (2DES) spectrum of the film on a flat gold surface at  = 0 fs.Bottom: Cross section along the detection energy for excitation at 1.591 eV showing the characteristic dispersive line shape of the J-aggregate exciton.Inset: One-quantum (|⟩) and two-quantum (|⟩) excitations contributing to the 2DES exciton peak.d) Waiting time dynamics at the maximum (red) and minimum (blue) of the 2DES exciton peak, showing incoherent relaxation dynamics on a 100fs time scale.e) Angle-resolved linear reflectivity spectra reveal the dispersion relations of upper (UP) and lower (LP) polaritons together with an angle-independent peak of "uncoupled" excitons (X).f) Finite-difference time domain (FDTD) simulation of the angle-resolved reflectivity.

Figure 2 :
Figure 2: Experimental 2DES maps of the J-aggregate-nanoslit array for selected incidence angles at waiting time  = 0 fs.Insets: Cross sections at selected excitation and detection energies.a) 2DES map for  = 18° displaying diagonal and cross peaks at the LP (1.571 eV) and UP (1.70 eV) energies.In addition to the dispersive "uncoupled" X peak (1.593 eV), cross peaks between both polaritons and the X transition are observed.The deduced resonance energies are marked as dashed lines.b) The same general features are also seen at  = 21°.The resonance energies of the polariton peaks are shifted according to their dispersion relation.c)For  = 27°, the detuning between LP and X is sufficiently large to resolve the dispersive line shape of all four diagonal and cross peaks of LP and X.At this angle, the intensities of the UP-related peaks are too weak to be seen.d) Dispersion relation of one-quantum (1Q) and two-quantum (2Q) excitations, as deduced from 2DES maps (open circles).The 1Q energies match well the linear dispersion relation (solid lines).2Q excitations of lower polaritons (2LP) and "uncoupled" excitons (XX) are extracted together with mixed LP/X and UP/X peaks, while 2UP and UP/LP peaks are lacking.

Figure 3 :
Figure 3: 2DES dynamics revealing polariton Rabi oscillations.a) Experimental 2DES map at  = 0 fs and  = 27°.b) Simulation of the 2DES map using the effective coupling Hamiltonian introduced in the main text.As in the experiment, diagonal and cross peaks with dispersive line shape involving the LP and X transitions appear, while the UP-related peaks are too weak.c) The waiting time dynamics of the (LP,LP) diagonal and (LP,X) / (X,LP) cross peak reveal pronounced Rabi oscillations with a period matching the peak splitting while such oscillations are lacking at the (X,X) diagonal.The dynamics detected on the positive and negative sides of the dispersive peaks are displayed in the upper and lower panels, respectively.d) Simulations of the waiting time dynamics for all peaks shown in c).e) Angle-resolved pump-probe dynamics detected at the positive peak of the dispersive LP resonance.f) Rabi oscillations periods 2ℏ/ (filled symbols) extracted from the 2DES waiting time dynamics for different incidence angles.UP-LP oscillations (red circle) are only resolved at  = 21°, near the crossing angle.All other periods (blue squares) match the X-LP splitting.The 2DES peak splittings (green triangles) follow the energy differences  between the corresponding UP, LP and X transitions, as deduced from the linear dispersion (open circles and lines).

Figure 4 :
Figure 4:Frenkel exciton simulations of plasmon-driven coherent exciton population oscillations (CPOs).a) A chain of 300 squaraine monomers with disordered site energies and nearest-neighbor coupling forms localized J-aggregate excitons.The monomers are coupled to a delocalized plasmonic mode with a spatially inhomogeneous coupling strength   that decreases in amplitude from 5.0 meV for excitons (  ) inside the narrow slit region to 0.5 meV for those (  ) outside the slits.This leads to the absorption spectrum in (b), displaying a single, delocalized LP mode and several disordered UP transitions together with an "uncoupled" exciton peak X. Excitons and plasmon mode (P) are resonantly excited by a 5-fs pulse.After excitation, the population of P displays oscillations at a period   = 2ℏ/(  −   ), given by the UP-LP splitting.Out-of-phase UP-LP oscillations are seen for the slit excitons.These are superimposed by slower oscillations at   = 2ℏ/(  −   ).These reflect CPOs from outside the slits into the slit region and back, as illustrated by the arrows (a).c) Population dynamics of the   ,   and  states after impulsive excitation, displaying oscillations at   and   .CPOs are seen on both,   and   , but are fully absent for  and   +   .Dynamics for the effective Hamiltonian are shown as dashed lines.d) Fourier transform of the populations displaying CPOs at   , Rabi oscillations at   , and weaker oscillations at   = 2ℏ/(  −   ).