Non-Hermitian Hamiltonians for linear and nonlinear optical response: A model for plexcitons Open Access

The coupling between light modes and radiative transitions of matter creates exciton polariton states, i.e., half-light and half-matter states. These were first demonstrated in atoms inside a microwave cavity and were used to study decoherence and quantum logic operations in atoms.^1–4 More recently, polaritonic states have also been observed in molecules in microcavities, having both infrared and visible wavelengths. They open new paradigms to influence chemical processes.^5–16 For example, the delocalization of polaritonic states is advantageous for long-distance energy transfer,^11,17–22 polaritonic chemistry,^23–26 and photocatalysis.^27,28 They can also be used for lasing and to study nonequilibrium condensates and phase transitions.^5,29 Using plasmonic nanoparticles instead of cavities confines the light to subwavelength dimensions. This allows creating billions of polaritonic systems in solution.³⁰ It allows many more interesting arrangements that exhibit strong coupling in the single or few molecule limit.³¹ 

The number of molecules in the interaction volume and the properties of the cavity determine the light–matter coupling strength, while the quality of the cavity and the molecular dissipative processes determine the coherence lifetime. The coupling strength and coherence lifetimes span a large parameter space that can be broken into different regimes, which are also associated with different phenomena.^32,33 For example, for weak coupling, we observe a Purcell enhancement of the fluorescence rate, while for strong coupling, we see a Rabi splitting, which is reflected in a splitting of the absorption into two peaks termed upper and lower polariton branches. However, such an apparent splitting is also possible in the presence of an interference that appears already in the weak coupling regime (also called electromagnetically induced transparency or a Fano interference³³). Finally, the ultrastrong coupling and deep strong coupling regimes are associated with the breaking of the rotating-wave approximation and are difficult to reach using molecules in cavities but have been observed with plasmonic lattices.³⁴ 

It is important to discuss dissipation when considering polariton transitions. While in atoms, the sources of dissipation are few and well understood, an intense experimental and theoretical effort has been necessary to assign the timescale and processes of coherent and incoherent dynamics in molecular polaritons.^30,35–40 Briefly, time-resolved spectroscopy experiments have detected coherent dynamics in the form of Rabi oscillations via transient absorption³⁷ and established the possible relaxation pathways: the upper branch can decay via radiation damping to the ground state,⁴¹ inject into the dark plexciton states,^42,43 or decay directly into the lower branch.⁴⁴ The dark states themselves can decay to the lower polariton branch, with some degree of back-transfer (from lower polariton to dark reservoir) taking place.⁴⁵ The long lifetime of the dark states is responsible for the excited state stabilization with decay times for up to several μs.^42,46,47 In systems, where the cavity mode is confined to a plasmonic nanoparticle, the dynamics are strongly influenced by electron–electron and electron–phonon scattering inside the metal nanoparticle. Notably, a direct coupling between the molecular dark states and metal surface states appears to limit the lifetime of the dark states to that of the electron–electron scattering time.⁴⁰ 

Simulations of the optical response are crucial to correctly interpret the experiments and more so in time-resolved spectroscopy experiments, such as two-dimensional coherent spectroscopy (2DCS) or transient absorption. Coupled oscillator models work well to understand linear optical absorption.⁴⁴ For transient absorption, they can also work as long as the experiment can be simulated as two linear experiments: the bleach of the ground state and the absorption of an excited state dominated by the signal from the remaining molecules in the ground state. This is true when the main nonlinearity is a Rabi contraction whereby the Rabi splitting of the ground state is reduced because of the molecules that have been excited.³⁷ However, they are not sufficient for simulating a third-order signal as, in general, all orders of the field appear together, and it is not straightforward to disentangle them. Input–output theory has also been used to compute with a very good agreement the experimental two-dimensional infrared spectroscopy (2DIR) signal of molecules in infrared cavities.²⁵ The Rabi contraction and softening of the vibrational mode (for 2DIR) have been proposed as the mechanisms of nonlinearity. Similar physics has been found in polaritons formed from molecules adsorbed on plasmonic lattices.³⁷ In an earlier work, we have simulated the two-dimensional electronic spectroscopy response of molecules coupled to plasmonic nanoparticles using an extension of the response function formalism to non-diagonal non-Hermitian Hamiltonians, successfully reproducing the spectra at early times using a Rabi contraction nonlinearity and at late times using a thermal expansion of the nanoparticle.⁴⁰ Gu et al. also employed a non-Hermitian Hamiltonian to calculate several nonlinear experiments, finding significant modifications of the energy levels as well as the selection rules.⁴⁸ Models of disordered polaritons have suggested the existence of exceptional points in these non-Hermitian Hamiltonians, which vary according to the degree of disorder,⁴⁹ and these have been measured experimentally in infrared cavities.⁵⁰ Non-Hemiticity of the Hamiltonian can arise from a number of different scenarios not only due to energy relaxation but also from particle decay, for example, in photoionization. Moiseyev et al. proposed a more realistic model of molecules, including the ionized continuum states, and studied their modification inside a cavity.⁵¹ 

There are two difficulties related to the description of the dynamics of hybrid plasmonic systems. First, the description of the optical response in terms of non-diagonal non-Hermitian Hamiltonians is not typical in the approach for calculating 2D spectra in molecular systems. Usually an excitonic basis is used, where diagonal fluctuations induced by the bath are added and result in a line shape function after bath-averaging.⁵² Second, plasmonic systems add a rich dimension of dynamics in terms of the non-equilibrium distribution of hot electron–hole pairs, formed inside the metallic band structure after plasmon decay. The dynamics inside the metallic band observed in transient absorption experiments has been beautifully demonstrated in pure metals;⁵³ however, the description has not yet been extended to mixed molecule-metal systems. One complication is the co-existence of coherent effects related to the Rabi oscillation between two discrete transitions along with the incoherent scattering of continuous non-equilibrium distribution of electrons in the metal. The latter phenomena is described by the dielectric constant of an electron gas whose temperature is described by the two-temperature model⁵⁴ and has not been adapted to handle the molecular states.

In this work, we consider the linear and nonlinear optical response of non-diagonal non-Hermitian Hamiltonians. We begin by providing a derivation of the non-Hermitian Hamiltonian and discuss its regime of validity. We then describe the decomposition of the linear and third-order optical response using double-sided Feynman diagrams. We conclude with a discussion of the polariton branch energy structure in the weak and strong coupling regimes and their distinctive spectral signatures.

The Hamiltonian from Eq. (1) can be derived from classical arguments⁵⁶ or obtained from the quantization of the modes of the electromagnetic field sustained by a spherical metallic particle.⁵⁵ In  Appendix A, we show a derivation suitable for the parameters of plexcitons explored earlier,⁴⁰ while in  Appendix B, we show a derivation obtained by projecting out continuum states. However, the effective non-Hermitian operator remains an incomplete description for time-resolved experiments as it cannot capture the dynamics occurring in the implicit degrees of freedom, i.e., the degrees of freedom that when removed induce the non-Hermiticity. In our case, these correspond to the excited electron–hole pairs in the metallic conduction band. The non-Hermitian Hamiltonian can be used for calculating the linear response, but it is not justified for higher-order signals. The dephasing of the optical coherence can give a width to the line shape; however, when this width arises, for instance, from decay into a different excitation (i.e., electron–hole states), these will contribute to the excited-state absorption in a way that cannot be handled by non-Hermiticity alone. We provide a justification for its use under restricted conditions in the next section [optical response, Eq. (8)] and derive the expressions applicable to the general case.

In general, the effective operator H_eff(z) is nonlinear in that it depends on the frequency parameter z. The physical meaning of this is that the particle density can transfer from the P partition to the Q partition (and back), which can play a role in the dynamics. This goes beyond the scope of this article, although the description in the wideband approximation will be valid at short enough times. Before understanding the spectral signatures depending on the non-Hermiticity, we isolate the important regimes of Eq. (1), which can correspond to different symmetries of the eigenvalues. As mentioned, non-Hermiticity can arise in a number of settings and is usually connected with a manifold of states that are not explicitly described. The physical effect of the imaginary part of the energy is to destroy (when negative) or create (when positive) particle density. The linear response of a system is not sensitive to what happens to the particle density that the non-Hermitian part of the Hamiltonian destroys. The fate of the particle density only becomes important in nonlinear response (most easily accessible in time-resolved experiments) when the particle density that leaves the system can come back into it. This can be captured exactly by Eqs. (6) and (7), provided that the effective operators can be calculated exactly.

The behavior of polaritonic states can be categorized in regimes according to the strength of the light–matter coupling, the magnitude of the dephasing rates, and detuning.³² We limit ourselves to coupling strengths for which the rotating-wave approximation is valid. For isoenergetic plasmonic and molecular transitions (ω_P = ω_J), the non-Hermitian Hamiltonian [see Eq. (A1)] has an exceptional point at S = 2g/(γ_P − γ_J) = 1,⁵⁶ taken as a boundary between the weak (S < 1) and strong (S > 1) coupling regimes. An additional dimensionless parameter, the cooperativity C = g²/γ_Pγ_J, has been proposed to separate the region of weak coupling C < 1 and the region where interference processes are present C > 1.^33,60 We will systematically analyze the first and third-order optical response above and below the exceptional point with six parameter families cases marked in Fig. 2(a). Point A (red) corresponds to the weak coupling regime for C = 0.20 and S = 0.21, point B (green) corresponds to the interference regime with C = 1.8 and S = 0.63, and point C (blue) corresponds to the strong coupling with C = 11.25 and S = 1.58. Points D (C = 0.25, S = 0.67), E (C = 0.44, S = 1.6), and F (C = 1.78, S = 3.2) are meant to explore a small region where we can have condition S > 1 but C < 1 and is outlined in  Appendix C. In all cases, the molecular dephasing is set at γ_J/ω₀ = 0.01. As mentioned previously, we consider two limiting cases where the plasmonic transition couples predominantly to the far field (i.e., μ_J = 0) and where the molecular transition couples to the far field (μ_P = 0).

As has been pointed out, for example, in Ref. 56, the behavior is different on either side of the exceptional point. We make this explicit by plotting the eigenvalues of Hamiltonian (A10) and eigenvectors as a function of the coupling strength (Fig. 3). For S < 1, the real part of the eigenvalues (i.e., the transition frequencies) is equal, while the imaginary part (i.e., the dissipation) is not. The two polariton branches are either predominantly plasmonic or predominantly molecular, and the distinction is carried by the dephasing rate of the branch (imaginary part of the eigenvalue). At the exceptional point, both eigenvalues and eigenvectors coalesce and become indistinguishable. For S > 1, the real part of the eigenvalues begins to split (i.e., Rabi splitting), while the imaginary part is now identical for both branches. The two polariton branches are now (and for all stronger couplings) an equal mixture of plasmon and exciton. We analyze in more detail the first and third order optical response of the non-Hermitian Hamiltonian from Eq. (A1) using the double-sided Feynman diagram formalism⁵⁸ in order to analyze the expected spectral signatures.

We consider the nth order response of the density matrix using Eq. (4) and calculate the signal of two-dimensional electronic spectroscopy and linear absorption in the impulsive limit.⁵⁷ 

We can decompose the absorption spectra into each polariton contribution, i.e., plot each term of the sum in Eq. (9) separately. In Fig. 4, we show the total absorption (green, solid) decomposed into the upper (blue, dotted) and lower (yellow, dashed) polaritons for different values of the coupling corresponding to points A, B, and C of Fig. 2. We consider two limiting cases, one where the plasmon predominantly couples to the far field (top row) and one where the J-aggregate predominantly couples to the far field (bottom row). The Euler plane of the insets show the transition dipole moments for each polariton branch, μ₋, μ₊. We first consider the case where the plasmon couples to the far field (top row). We can see that for all parameters, there is the appearance of a peak splitting; however, below the exceptional point, it arises from the destructive interference of a broad transition and a narrow transition, while above the exceptional point, it arises from a Rabi splitting. Below the exceptional point (A and B), μ₊ is purely real, while μ₋ is purely imaginary. The branch with the purely imaginary transition dipole moment will have a negative contribution to the total absorption and cause an interference dip. We note that the individual contributions of the branches are not physical by themselves so that a negative contribution does not mean a stimulated emission but rather intuitively describes the interference process. The dip is clearly seen for spectra (a) and (b) where the yellow (dotted) line cancels the plasmonic absorptive contribution. This results in a spectrum that seems to have two apparently separate peaks well below the condition for Rabi splitting. However, while the spectra might look similar above and below the exceptional point, the nature of the excitation and consequently the excited state dynamics are expected to be very different. As the coupling increases, the transition dipole moments are no longer purely real or purely imaginary but have both real and imaginary components (see the Euler plane insets of Fig. 4). The contribution of each branch now amounts to two true peaks separated by Rabi splitting. While qualitatively the final spectra look similar, the nature of the resonances is not. For example, exciting on the right or left of the dip in the interference regime excites the same homogeneously broadened excited state, while in the Rabi splitting regime, exciting on the right or left addresses a different excited state (upper polariton or lower polariton). This has been recognized theoretically and experimentally in the fluorescence signatures of quantum dots coupled to plasmons.⁶¹ The case where the J-aggregate couples predominantly to the far field is strikingly different. Below the exceptional point, we see no evidence of peak splitting since the negative contribution is too broad. Above the exceptional point, we see the expected Rabi splitting although with less contrast than when the plasmon couples to the far-field [Figs. 4(d) and 4(e)]. The interference dip is only observed if it is the plasmon that couples predominantly to the far-field (as is usually the case). The depth of the minimum is modulated by the molecular dephasing rate and the degree to which the molecular aggregate also couples to the far-field. The set of points D, E, and F show that having the condition S > 1 is not enough to observe a peak splitting. In the region where S > 1 but C < 1, the system where the molecular aggregate couples more strongly to the field does not show a Rabi splitting (see Appendix E, Fig. 8). While the features of destructive interference are more easily identified, there will also be regions of constructive interference away from the dip. Thus, both the cooperativity and strength parameters are important for classifying the features of the spectra. The difference between interference and Rabi splitting has also been observed in simulations with disordered molecules.⁴⁹ There, the narrower cavity absorption generates an interference pattern in the heterogeneously broadened molecular absorption, while in the case presented here, the narrow molecular absorption generates an interference pattern in the homogeneously broadened plasmonic absorption.

It is appropriate to clarify an important point concerning the interference process. This dip is often referred to either as a Fano interference or an electromagnetically induced transparency (EIT), and both processes denote different physics and spectral signatures. The simplest system where EIT is observed is a three-level Λ system consisting of a ground state manifold with two levels and one discrete excited state. Its spectral signature is a frequency region of suppressed absorption. The Fano interference appears in a structure akin to the Λ system where one of the levels becomes a continuous manifold of levels and the Fano line shape is characterized by a distinctive asymmetric line shape which also includes a frequency region of suppressed absorption.⁶² The plexciton Hamiltonian is more closely related to a Λ system where instead of having two light-fields coupling each ground state with the excited state, we have an external light field coupling the ground state with the plasmon excitation and radiative coupling connecting the excited plasmon with the excited J-aggregate [Fig. 5(a)].

Figure 5 shows the absorption spectrum for a Λ system and a Fano system as a function of the asymmetry parameter, which is related to the ratio of the different coupling elements.⁶² For Fano, the asymmetry parameter is $q=μJπμPg$⁠. We choose the same definition for an effective asymmetry parameter for the Λ system although its physical meaning is not exactly the same.⁵⁸ For the Fano model, the limiting cases μ_J = 0 (q = 0) and μ_J → ∞ (q → ∞) correspond to anti-Lorentzian and Lorentzian line shapes, respectively. For finite values of μ_J (or q), the Fano system shows distinctive asymmetries that reverse sign as the sign of q is reversed (middle panel) while shifting the condition of destructive interference to the left or right of the zero detuning condition. The Λ system [Fig. 5(a)] shows a similar behavior for the limiting cases q = 0 and q = ∞ if we allow that far away from the resonance condition, the absorbance vanishes instead of going to a finite value as in the Fano model. However, for finite values of the effective asymmetry parameter and unlike the Fano model, we do not find an asymmetry. For finite q, the contrast at the point of destructive interference is reduced, and the plots are identical for q and −q. This point will become important in the analysis of nonlinear signals. We mention that real systems often do show asymmetries, although this is due more to the detuning of the plasmonic transition with respect to the molecular bright transition, and not a true Fano asymmetry.

We now turn to a more detailed analysis of the different features that appear in the spectrum of third-order spectroscopies of non-Hermitian Hamiltonians.

Figure 7 shows the spectrum for weak 7(a) and strong 7(b) coupling along with their decomposition into diagonal elements ((UP,UP) and (LP,LP)) and cross-terms ((UP,LP) and (LP,UP)). For the weak coupling case, the diagonal elements are positive while it is the cross-terms that provide the negative signal and split the main plasmonic resonance into four regions. That the negative features arise from the cross-terms between upper and lower branch is another indication that the dips are a result of an interference process between the plasmonic resonance and the molecular resonance and not a splitting of the energies. Having transition dipole moment pairs μ_±g, μ_g± which are not complex conjugates of each other introduces phase shifts in the resonances. As a result they acquire dispersive character along with the negative features. These negative features do not appear in the GSB and SE contributions in molecular systems where non-diagonal non-Hermitian Hamiltonians are not common. However, they should not be interpreted as ESA features. In the linear response, the negative features are always compensated by positive ones so that absorption is always positive. Instead, our physical picture is that they are the inevitable expression of destructive interferences when we decompose the spectrum into some contributions that are strictly positive, and so require some negative contributions to suppress the absorption (this is also more clearly seen in the linear optical absorption case).

The strong coupling regime is qualitatively different. Each of the four contributions now appear at a different position in the spectrum because of the Rabi splitting. A notable difference between the GSB signal between the two regimes is that the interference regime has a D₄ fourfold reflection symmetry, while the strong coupling regime possesses a C₂ rotational symmetry (as we show later when both μ_J ≠ 0 and μ_P ≠ 0 some symmetries are lost).

We now analyze the effect of having both μ_P ≠ 0 and μ_J ≠ 0. Figure 9 shows the GSB for different values of the coupling strength g (columns) and also different values of the ratio μ_J/μ_P (rows). Along the columns, we observe the same qualitative difference between weak and strong coupling, namely, the transition in the symmetry of the line shape. In this case when μ_J ≠ 0, we lose some symmetries. Strikingly, while the linear response for different values of μ_J/μ_P does not show an asymmetry as the effective asymmetry parameter changes, the third-order response shows a clear asymmetry along the detection dimension. Projections along the excitation dimension (lower row of Fig. 9) show that this asymmetry will be visible in transient absorption. This suggests that differences in the linear spectrum and the 2D projection onto the detection dimension can reveal information on the Hamiltonian. Asymmetries in the linear response will denote different transition frequencies between the plasmonic and molecular parts, while additional asymmetries in the transient absorption will indicate finite values of μ_J.

The description of physical phenomena using non-Hermitian Hamiltonians is pervasive and extends far beyond the physical systems described here.⁶⁶ These can exist in metamaterials, engineered Floquet states, and gain/loss media, in general. The treatment described in the previous section (Third order response) applies to all of these systems as long as they have a non-zero nonlinear response. Our findings suggest that it is interesting to consider two-dimensional spectroscopy (or more generally higher orders) as tools that can distinguish between different symmetries of the eigenvalues by encoding them in symmetries of the signal.

Rigorously, a Lindblad operator should be used instead of a non-Hermitian Hamiltonian. However, the quantum jump operator that restores the particle density back to the ground state can be encoded by a judicious choice of Feynman diagrams.

There exist many models to describe the response of emitters coupled to plasmonic resonances, which can fit absorption spectroscopy, scattering, or photoluminescence data well. For example, the coupled oscillators model can explain the regimes of interference and Rabi splitting and reproduce the measured line shapes. However, the formalism is not apt to describe nonlinear spectroscopy because of the impossibility of isolating the signals to a certain order in the field. Ab initio approaches involving the solution of Maxwell’s equations around a metallic nanoparticle result in Fano-like equations with slightly modified detunings or in non-Hermitian Hamiltonians for molecules coupled to plasmonic resonances. These expressions are sometimes reflective of scattered signals and not absorptive signals,^33,67,68 as such they are suitable for linear response but not nonlinear signals. The use of non-Hermitian Hamiltonians as justified above solves this problem and allows the use of the response function formalism and double-sided Feynman diagrams to describe interference phenomena and, in particular, to suggest its use for discerning between interference processes and Rabi splitting. We mention that it is not the only option, and previous works that compare scattering and photoluminescence (PL) show a marked difference in the Rabi splitting with significant reductions in the PL spectra compared to the scattered light.³³ For regions of interference, the peak splitting entirely disappears from the PL spectrum, while it is visible in the scattered spectrum,^61,69 resulting in a diagnostic tool, which can be applied in the case of few emitters coupled to a plasmonic nanoparticle.

One of the findings of this paper is that the third-order signal is qualitatively sensitive to changes in the symmetry of the eigenvalues, unlike the first order signal. The symmetry of the linear and nonlinear signal with respect to reflection along ω_τ = ω₀ and ω_t = ω₀ can be analytically obtained by the symmetries in the eigenvalues and eigenvectors below and above the exceptional point. We define the detunings with respect to the transition frequencies as δ_i = ω_i − ω₀ and the reflection operations that take δ_τ → −δ_τ and δ_t → −δ_t as $Oτ()$ and $Ot()$⁠, respectively. Below the exceptional point, we have the following properties of the eigenvalues and transition dipole moments: δ₊ = δ₋ = 0, μ₊ is purely real, and μ₋ purely imaginary. It is straightforward to verify that $Ot(S(1))=S(1)$⁠. In the strong coupling limit, we have that $μ−2=(μ+2)*$ so that it is also verified that because δ₋ = −δ₊, $Ot(S(1))=S(1)$⁠, and the signal is symmetric with respect to the central transition frequency ω₀.

In a similar fashion, we can calculate the symmetries of the third order signal described by Eq. (11) using Eqs. (2) and (3). Below the exceptional point, we have that $Ot(S(3))=Oτ(S(3))=S(3)$⁠, and above it, we recover that $Ot(Oτ(S(3)))=S(3)$⁠. This work raises the possibility of using the symmetry of the line shape of nonlinear spectroscopy to study the symmetry of the eigenvalues of the Hamiltonian. In real systems, we expect to have deviations from these idealized Hamiltonian; however, we do expect the symmetries of the spectra to be recognizable.

We have derived a non-Hermitian Hamiltonian for the calculation of the linear and nonlinear optical response. Decomposing the Hamiltonian into bright and dark partitions and expressing the linear and nonlinear optical response in terms of effective operators provide a new framework for systems with large manifolds of dark states. An analysis of the linear optical absorption reveals that the exceptional point separating weak from strong coupling regime is clearly illustrated when decomposing the total signal into individual Feynman diagram contributions. Below the exceptional point interference effects are indicated by negative contributions to the absorbance, while the Rabi splitting appears above the exceptional point. The decomposition of the nonlinear signal also provides additional insight. We find that interference and Rabi splitting regimes have different symmetries in the two-dimensional maps, an indication that cannot be obtained from linear absorption alone. The connections outlined in this article between the symmetry of the spectral signatures and the nonlinear response can open new ways of thinking for classifying symmetries in the eigenvalue/exciton structure of complex materials.

The authors have no conflicts to disclose.

Daniel Finkelstein-Shapiro: Conceptualization (lead); Formal analysis (lead); Funding acquisition (equal); Writing – original draft (lead). Pierre-Adrien Mante: Conceptualization (supporting); Formal analysis (supporting); Funding acquisition (equal); Investigation (supporting); Writing – original draft (supporting). Sinan Balci: Conceptualization (supporting). Donatas Zigmantas: Conceptualization (supporting); Formal analysis (supporting); Investigation (supporting); Writing – original draft (supporting). Tõnu Pullerits: Conceptualization (supporting); Formal analysis (supporting); Funding acquisition (equal); Investigation (supporting); Project administration (supporting); Writing – original draft (supporting).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

There are two main line-broadening mechanisms. On the one hand, the optical coherence of the molecular component dephases with fluctuations of the transition frequency (with a time constant of 40 fs⁴⁰ for TDBC, which is much faster than the decay induced dephasing into molecular dark states). To capture this dephasing mechanism, it is enough to consider the stochastic fluctuation of the transition frequency of the molecules. The plasmon coherence, on the other hand, dephases via Landau damping whereby the plasmon decays into neutral hot electron–hole states inside the metal. For this process, we can consider the coupling of the plasmon transition to a continuum of electron–hole states via a plasmon–electron coupling element.⁷⁰ These two are very different mechanisms. The lifetime induced broadening due to decay of the plasmon into hot electron–hole states is not a thermal process and has negligible temperature dependence, in stark contrast to line broadening in molecules induced by fluctuations of the bath, which have strong temperature dependence and can be significantly reduced at low temperatures.

With Eq. (A1), we model the molecular aggregate as N two-level systems, which are coupled to their nearest neighbor with strength K. The plasmon is represented by a boson coupled to each two-level system with a coupling strength g₀. The plasmon excitation can decay into hot electron–hole pairs that form a continuum we label by its momentum k. The decay is mediated by the plasmon–electron interaction V_k.⁷⁰ Both the molecules and the plasmon couple to the far-field with transition dipole moments μ₀ and μ_P, respectively (see Fig. 1). These couplings to the light field are analyzed in the main text for the two limiting cases when μ₀ = 0 and when μ_P = 0.

There are three energy scales of the problem: the coupling strength $g=g0N$⁠, the nearest neighbor coupling K, and the transition energy fluctuation amplitudes Δω_ii. To derive the current effective Hamiltonian, we assume the limit Δω_ii ≫ g, K, which is appropriate for the study of interference effects. For simplicity, we also assume here the limit of vanishing intermolecular coupling K. The approach is to diagonalize the Hamiltonian for the molecular component, average over fluctuations over the bath, couple to the plasmonic part, and project onto the bright states.

The Hamiltonian [Eq. (A6)] can be already used to calculate the optical response of the plexcitons; however, the resulting expression for the signal will contain an explicit description of dark states that do not couple to the far-field and so is inefficient. In order to transform to a fully excitonic picture, we must also diagonalize the plasmon-metal Hamiltonian. However, it is much more preferable and transparent to obtain effective operators for the bright transitions.

In this appendix, we explore the region marked by points D, E, and F of Fig. 2 that are in the strong coupling regime (S > 1) but not in the weak coupling regime with respect to the interference parameter (C < 1). Figure 10 shows that point E shows the presence of Rabi splitting for the case where the plasmon couples to the far field but not for the case where the molecular aggregate couples to the far field. We see here the value of the cooperativity parameter. In the case where the molecular component couples predominantly to the far field, the Rabi splitting is absent even though we are in the strong coupling limit if we are in the low cooperativity regime. It takes both parameters to be larger than one to observe some the Rabi splitting.

We can simulate the full spectrum by choosing a model for the excited state absorption (ESA). If the excited state absorption is dominated by the molecules remaining in the ground state after the pump pulse pair, we can simulate the full spectrum by the Feynman diagrams shown in Fig. 6. The differences found with experimental measurements can be ascribed to the finite bandwidth of the pulses.