{'fraction_non_alphanumeric': 0.03374771053349441, 'fraction_numerical': 0.008573321382642921, 'mean_word_length': 4.672413793103448, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Plasmonic Purcell Effect in Organic Molecules

By means of quantum tensor network calculations, we investigate the large Purcell effect experienced by an organic molecule placed in the vicinity of a plasmonic nanostructure. In particular, we consider a donor-{\pi} bridge-acceptor dye at the gap of two Ag nanospheres. Our theoretical approach allows for a realistic description of the continua of both molecular vibrations and optical nanocavity modes. We analyze both the exciton dynamics and the corresponding emission spectrum, showing that these magnitudes are not accurately represented by the simplified models used up to date. By disentangling the molecule coupling to radiative and non-radiative plasmonic modes, we also shed light into the quenching phenomenology taking place in the system.

The Purcell effect [1] lies at the core of quantum electrodynamics, as it reveals that the radiative properties of any quantum emitter are not inherent to it, but also depend on the electromagnetic (EM) vacuum in its surroundings. A few decades ago, this phenomenon, and in particular, the pursuit for the full inhibition of spontaneous emission at EM band gaps, was one of the driving forces behind the development of photonic crystals [2]. More recently, much interest has focused on metallic nanocavities [3,4]. The large and spectrally complex photonic density of states associated to plasmonic resonances allows an unprecedented control over spontaneous emission in these nanostructures [5,6]. In recent years, different light sources have been used to probe plasmonic Purcell enhancement phenomena, such as quantum dots [7,8], solid-state color centers [9,10] or passing electron beams [11,12].
Within the TLS approach, only a few timescales play an important role, and in particular, the relative values of the plasmon-exciton exchange rate and their respective losses can be used to distinguish between the weak-coupling and strong-coupling regimes in a fundamental analysis. In contrast, in molecules, the vibronic coupling and the vibrational dynamics introduce several new timescales, such as the coherent nuclear oscillation period or the rate of vibrational energy dissipation and thermalization time, and the simple dichotomy between weak and strong coupling can be expected to give way to a richer phenomenology. For the case of not too large Purcell factors, radiative decay is typically slower than the vibrational relaxation and thermalization times (on the order of 1 ps or less [31]), and molecular emission can be assumed to proceed from the lowest excited molecular state, with all vibrational modes in thermal equilibrium. This assumption gives rise to a Fermi's golden rule (FGR)-based approach in which the emission spectrum of the organic molecule near a plasmonic structure is given by the product of its free-space spectrum and a frequency-dependent Purcell enhancement [19,20,32]. However, as plasmonic structures can achieve Purcell factors on the order of 10 6 [28], radiative decay rates can be decreased from their free-space values in the nanosecond range to femtoseconds, such that radiative decay does not proceed from vibrationally relaxed molecules, the FGR approximation breaks down, and non-equilibrium effects play an important role. More recently, cavity-QED approaches have been able to incorporate the interplay between the exciton and one (or a few) vibrational modes [33][34][35][36][37]. However, a realistic description of typical organic molecules could require considering hundreds of vibrational modes.
In this Letter we present an accurate theoretical framework that is able to treat on an equal footing the electronic, vibrational and plasmonic degrees of freedom associated with the Purcell effect experienced by organic molecules placed within a metallic nanocavity. Our numerical scheme is based on a quantum tensor network (TN) method [38,39], which helps provide a complete picture of this phenomenon. In addition, it allows for a validation of the simple models introduced above and also for an in-depth study of the influence of the large number of vibrational modes on the spontaneous emission rate and emission spectrum. We show that while the TLS model is able to reproduce the very short timescales of the spontaneous decay, the FGR approach accounts only for the behavior at very long times. On intermediate time scales, a model that incorporates a single vi-brational mode is able to extent the validity of the TLS to slightly longer times. However, we demonstrate that it is mandatory to account for the whole set of vibrational modes to capture most of the time dynamics and the corresponding spectrum.
As a case study, we consider the excited-state dynamics of a single organic molecule placed at the gap center of a silver nanosphere dimer, as sketched in the inset of Fig. 1(a). This plasmonic structure resembles the two geometries in which the largest Purcell factors have been reported, the so-called nanoparticle on mirror [4] and bowtie antenna [16] geometries. In this work we have chosen a donor-π bridge-acceptor (D-π-A) organic dye, labeled CPDT [40], as a prototypical organic molecule because it displays a significant transition dipole moment and a large Stokes shift (> 0.4 eV), yielding excellent absorption and emission capabilities. The Hamiltonian for this hybrid system can be written as [h = 1] where σ − , b ω and a i,ω denote the lowering operator of the molecular electronic transition, annihilation operator for the molecular vibration and i-component of the plasmonic mode at frequency ω. Here the index i is either nr or r, representing the non-radiative and radiative plasmonic modes, respectively. Notice that only the latter can be detected in the far-field. The first line in Eq.
(1) describes the free-standing molecule, where we have used a Holstein-type Hamiltonian [41] accounting for the continuum of vibrational modes, with coupling strengths λ ω . The exciton frequency ω e corresponds to the vertical transition from the vibrational ground state of the electronic ground state. The second part in Eq. (1) describes the continuum of radiative and non-radiative plasmonic modes, and their coupling, weighted by g i,ω , to the molecular exciton. This plasmon-exciton coupling can be encoded in the spectral density J i (ω) = g 2 i,ω , in a similar way as the vibrational spectral density J v (ω) = λ 2 ω does for the vibronic coupling.
We numerically solve both the excited-state dynamics and the emission spectrum with our quantum TN method. To apply this TN approach, the two continua in Eq. (1) are transformed to a chain form in which each continuum (EM and vibrational) is represented by a chain of nearest-neighbor coupled oscillators, with only the first oscillator coupled to the exciton. More details of our TN numerical framework can be found in the Supplemental Material [42]. The first site in the vibrational chain is usually termed the reaction coordinate (RC), whose frequency is given by with λ RC = J v (ω)dω being the coupling strength between the RC and the exciton. If the molecule is initially in its excited state, the coupling to the plasmonic chain leads to its decay whereas the coupling to the vibrational chain dresses the electronic state.
To shed light into the effect of molecular vibrations and uncover the relevant timescales in the spontaneous emission process, we compare our TN numerical results against three simplified models. As the simplest choice, by discarding all the molecular vibrations, the standard TLS model predicts a decay only dictated by the plasmonic environment, is the total plasmonic spectral density. In a second step, by keeping only the RC within the vibrational chain, we can derive an approach, dubbed here as the single vibration mode (SVM) approximation, in which the RC comprises all the vibrational response. As a difference, within the FGR model [42], the spontaneous emission rate is simply calculated as the spectral integral of the product of J p (ω) with the so-called lineshape function, D em (ω), which represents the available optical transitions connecting the ground and the excited electronic states of the molecule under the assumption that the vibrational modes are in thermal equilibrium, leading to As commented above, we consider a silver nanosphere dimer as an example of plasmonic cavity, with a 1 nm gap and 20 nm radius, embedded in a matrix with refractive index n D = 2.1. The dielectric constant for silver is taken from Ref. [43]. As shown in Fig. 1(a), the plasmonic spectral density is characterized by a radiative dipole mode located at around 1.5 eV, a quadrupole mode at 2.1 eV, and a dominant non-radiative mode, the so-called pseudomode [44], emerging at 2.9 eV. A single CPDT molecule is located at the gap center and we assume that its transition dipole moment (modulus µ = 0.1 e·nm) is pointing along the line that connects the two nanospheres. In the numerical calculations, we assume that the initial state originates from a Franck-Condon excitation, i.e., a vertical transition that could result after an ultrashort laser pulse excitation of the electronic ground state.
Electronic structure calculations within density functional theory and its time-dependent version are performed to obtain the vibronic coupling constants of the CPDT dye for the displaced harmonic oscillator model. Both the vibrational spectral density and lineshape function of the CPDT molecule used in our calculations are shown in Fig. 1(b). In order to explore different regions of the plasmonic spectral density while keeping the vibrational structure unchanged, we artificially shift the (vertical) exciton energy of the molecule to different values.
We focus first on the excited-state dynamics. Figure 2 shows the evolution of the exciton population, σ + σ − , evaluated with the TN method at four exciton frequencies. It is compared against the results of the three different models (TLS, SVM and FGR) for the same cases. First, it is evident that the TLS model largely fails to describe the decay dynamics for all the cases depicted in Fig. 2. This failure highlights the importance of going beyond the TLS approach when defining the coupling regimes in the interaction of organic molecules with nanophotonic structures, implying that the simple distinction between weak and strong coupling has to be modified for organic molecules. Still, it is interesting to note that TLS seems to be valid up to times of about τ 1 ≈ 4 fs (grey arrows in Fig. 2). This is because, soon after its generation, the exciton vibrational wavepacket remains in the vicinity of the Franck-Condon region before exploring the potential energy surface of the excited state [45]. Hence, within a very short timescale, vibrations do not play any role yet and the TLS model describes the molecular decay. This time scale, τ 1 , can be estimated as a fraction of the period associated with the RC harmonic motion, T RC = 2π/ω RC , with τ 1 ≈ T RC /6 = 3.8 fs giving a good estimate of this initial timescale.
When the wavepacket initially moves away from the Franck-Condon region, its dynamics is dominated by the RC, as revealed by the accuracy of the SVM description for all ω e in Fig. 2. This regime holds until both dephasing and decay of the RC into other sites in the TN vibrational chain becomes important. We can estimate this second timescale, τ 2 , as the inverse of the coupling between the first (RC) and second sites of the vibrational chain, This estimation gives τ 2 = 1/γ d = 17.6 fs (brown arrows in Fig. 2). Following these arguments, it can be understood why the FGR model, in which all vibrational modes are taken into account, works better than the TLS and SVM models at sufficiently long times, as observed in Fig. 2. However, FGR assumes that the vibrational modes at all times are in thermal equilibrium and thus neglects the initial strongly non-equilibrium state and its coherent wave packet motion. The FGR approach thus fails to capture the short-time dynamics, and in particular cannot represent the weak oscillations observed in the TN calculations.
We next examine the frequency-dependent population of the EM modes, S em (ω) = i a † i (ω)a i (ω) . Notice that the sum extends over both radiative and non-radiative plasmonic modes, so we name this physical magnitude the near-field emission spectrum. Within the FGR approach, this quantity is independent of time and can be written as S em (ω) ∝ D em (ω)J p (ω) [42]. Figure 3 shows the near-field emission spectra calculated for the four exciton frequencies and the four theoretical approaches in  to ω e contribute to the spontaneous emission process. This leads to single-peaked spectra broadened by the plasmonic environment, very different to the spectra obtained within the TN framework. On the other hand, the SVM approximation largely fails to reproduce the TN-spectra for low exciton frequencies (1.5 and 2.0 eV) although it provides a reasonable approximation for higher ω e . This is in accordance with its better accuracy describing the excited-state decay dynamics for those frequencies (see Fig. 2).
In principle, one could have expected that the FGR approach should work better than the other two simplified models, as it incorporates the whole vibration spectrum via the lineshape function of the organic molecule. However, as observed in Fig. 3, this is not the case for the four chosen ω e . The FGR approach only works well when the coupling of the molecular exciton with its EM environment is weak enough such that vibrational thermal equilibrium is reached before emission takes place. However, for large plasmon-exciton couplings leading to huge Purcell factors as those associated with plasmonic fields, the molecular exciton decays so fast that vibrational equilibrium is not reached and the spectrum is significantly different from the stationary limit. In agreement with recent experiments [23,24], our TN calculations show that this fast decay can be utilized to strongly modify the branching ratio of the emission by organic molecules.
Finally, taking advantage of the capability of our TN theoretical framework to separate contributions from plasmonic radiative and non-radiative channels, here we address the interplay between the Purcell effect and the quenching phenomenon. Figure 4 shows the evolution of the population of EM modes, n i = dω a † i (ω)a i (ω) (i being either nr or r), as a function of exciton frequency and evaluated at two different times. At very short time scales (t = 5 fs, panel a), the TLS approximation reproduces the TN result. The non-radiative plasmon population is maximum at exciton frequencies at resonance with the pseudomode (around 2.9 eV). The radiative plasmon population becomes relevant at longer times, as shown in Fig. 4(b). In the TLS approximation, this population is maximum when the exciton frequency coincides with the dipolar mode (1.5 eV). On the contrary, the TN calculations reveal that, the maximum far-field emission takes place when the exciton frequency is slightly bluedetuned. This setup maximizes the spectral overlap between D em (ω) and J r (ω). Within this exciton frequency range, the effective quantum yield, shown in the inset of Fig. 4(b), can be as high 30%, despite the fact that the gap in the plasmonic cavity is 1 nm wide and quenching is expected to be dominant. For higher ω e , non-radiative components take over, completely quenching the far-field emission of photons by the molecule, as expected.
To conclude, we have applied an accurate quantum tensor network method, able to treat electronic, electromagnetic and vibrational degrees of freedom in lightmatter scenarios on an equal footing, to study the Purcell effect occurring when organic molecules interact with plasmonic nanocavities. Using this numerical framework we have tested different simplified models that have been used to analyse this phenomenon lately. We have found that a two-level description of the molecule fairly captures the first steps of the decay process and a simple approach based on the Fermi golden rule can describe its very long time behavior. However, an accurate treatment of the exciton coupling to all the vibrational modes in the organic molecule is mandatory to resolve intermediate timescales. As a consequence, none of the simplified models used up to date are able to provide a precise estimation of near-or far-field emission spectra in hybrid systems involving organic molecules and plasmonic cavities. As an interesting extension of the current work, the influence of non-linear vibrational mode couplings leading to vibrational energy redistribution could be investigated, as they are not included in the current Holstein-type model and are expected to become important on slightly longer timescales than treated here [46,47]. Our findings also reveal that significant effective quantum yield values can be achieved in situations where strong interaction to non-radiative plasmonic modes takes place and quenching is expected to dominate.

I. TENSOR NETWORK CALCULATIONS
Given that the plasmonic environment consists of both radiative and non-radiative parts, one may guess that the decay dynamics depend on the respective values of each components. However, as shown in the following, the excited-state dynamics of the molecules are governed by the total EM spectral density J p (ω) only. To this end, starting from Eq. (1) in the main text, we first introduce two sets of new plasmonic modes a ω = α nr,ω a nr,ω + α r,ω a r,ω , (S1) d ω = α r,ω a nr,ω − α nr,ω a r,ω , where the coefficients α i,ω are defined by α i,ω = g i,ω /g ω , with g ω = g 2 nr,ω + g 2 r,ω . Rewriting Eq. (1) in the main text with the new operators gives Here one can see that d ω acts as a dark mode that is decoupled from the exciton. Moreover, the original EM modes can be expressed in terms of the new mode operators as a nr,ω = α nr,ω a ω + α r,ω d ω , (S4) a r,ω = α r,ω a ω − α r,ω d ω .
Taking the non-radiative component as an example, the population in this branch can be calculated from a † nr,ω a nr,ω = α 2 nr,ω a † ω a ω + α 2 r,ω d † ω d ω +α nr,ω α r,ω ( a † ω d ω + a ω d † ω ). (S6) Since we assume an initial vacuum state for the EM modes and mode d ω is decoupled from the exciton, any expectation value related with d ω vanishes. With this argument we arrive to the following expression for both the radiative and non-radiative components, with the spectral density J i (ω) = g 2 i,ω and J p (ω) = g 2 ω = J nr (ω) + J r (ω). Therefore, the excited-state dynamics of the molecules depends only on the a ω operators with coupling strength g ω , and hence the total EM spectral density J p (ω).
To perform the tensor network (TN) calculations, we first rewrite the system Hamiltonian into a discrete form where we have omitted the d ω operators for simplicity. Note that after the discretization, b k and a l are dimensionless operators, while λ k and g l have the same dimension as frequency. The spectra S em (ω) = a † ω a ω defined in the main text can be obtained from S em (ω) = a † l a l /∆ω, with ∆ω being the frequency interval for the discretization.
To enable the treatment of Eq. (S8) by the density matrix renormalization group methods [S1], a crucial step is applying a chain mapping transformation. It transforms the Hamiltonian into a form containing one-dimensional chains with only nearest-neighbor interactions [S2, S3]. By doing so the system Hamiltonian can be written in an equivalent chain form Therefore, the exciton is now coupled to the vibrational and plasmonic chains. The former accounts for the dressing of the exciton by the molecular vibrations, while the latter leads to the exciton decay. Among many other chain parameters,ω v 1 , t v 1 and t v 2 , denoted by ω RC , λ RC and γ d respectively in the main text, are of vital importance for determining the timescales where the TLS and SVM approximations work.
The Hamiltonian in Eq. (S9) is formidable with a brute-force approach for a large number of modes. By expressing the system's wave function as a TN, i.e., a network of interconnected tensors with fewer coefficients, the computation demand can be greatly reduced [S4]. Therefore, TN methods provides an efficient way to simulate many-body dynamics and allow for the computation of a quasi-exact solution to the time-dependent Schrödinger equation for very complex systems. In this work we calculate the dynamics with the time-dependent variational matrix product states (TDVMPS) technique [S5-S7].
The main parameters for the TN calculations in this work are as follows. For the plasmonic part, the EM modes with frequency range from 0.6 eV to 6.0 eV have been considered, while for the molecular vibrations, all the modes of the CPDT dye have been included, lying within 0.1 to 500 meV. Both plasmonic and vibrational environment are discretized into L = M = 650 modes. Moreover, to perform the TN calculations, a bond dimension D = 22 and a time step ∆t = 5 a.u. have been used. All these parameters have been checked to guarantee the convergence of the results.
With regard to the organic molecule, electronic structure calculations within density functional theory and its time-dependent version were performed to obtain the vibronic coupling constants for the displaced harmonic oscillator model. The long-range corrected CAM-B3LYP functional [S8] was used together with the 6-31G(d) basis set to properly describe the charge-transfer transition. For the most stable molecular conformer of the CPDT dye [S9], ground and excited state harmonic frequencies were calculated and the Huang-Rhys factors were obtained considering Duschinsky rotation effects [S10]. All calculations were performed with the Gaussian 16 package [S11]. Note that the CPDT dye is closely related to the commercial C218, Y123 and Dynamo Red dyes widely employed in dye sensitized solar cells. faster than 1/γ M . This leads to a fixed probability distribution among states {|e, 0, M } during the whole decay process. In this case, the decay of the excited-state population is given by a single exponential as σ + σ − = e −γt , with the decay rate being

II. FERMI'S GOLDEN RULE
The above expression can also be rewritten as with the averaged lineshape function D em (ω) given by In this case, the near-field emission spectrum reads Strictly speaking, as we don't include any phonon redistribution mechanism in our Hamiltonian Eq. (S8), the decay should follow the multi-exponential dynamics discussed in case (i). However, as we will show below, there is no significant difference between both cases for the studied system. To calculate the decay rate in Eq. (S13) and (S18) and hence the lineshape function, detailed knowledge of the vibrational eigenstates and eigenenergies for both the excited and ground state potential energy surfaces is required. Here we focus on the case that the molecular vibrations can be approximately described by displaced harmonic oscillators as in Eq. (S8). With respect to the initial state, we consider a Franck-Condon excitation, the vertical transition from the vibrational ground state of the electronic ground state by ultrashort laser pulses. This leads to an initial distribution Now we compare the two limiting cases mentioned above. Figure S1 shows the comparison between the multi-exponential and single-exponential dynamics for an initial Franck-Condon excitation. To ease the numerical computation, only a few vibrational modes have been included. The two cases show slight differences if we consider one vibrational mode only. By increasing the number of vibrational modes, the calculations tend to converge and the differences between those two cases disappear. Moreover, due to the broadband structure and strong Purcell enhancement effect of the plasmonic pseudomode, it is easier to converge calculations for excitons coupled to this mode, meaning that the decay rates are not sensitive to the lineshape function. In contrast, the relatively narrower structure of the lower-order modes necessitate a more accurate description of the lineshape function, and hence the inclusion of more vibrational modes.
Although the brute-force calculation of the lineshape function based on Eq. (S14) and (S19) works quite well for a few number of vibrational modes, they are prohibitive for hundreds of modes, which is the typical number of vibrational modes for organic molecules. Fortunately, if fast thermalization takes place for the molecular vibrations due to its interaction with a thermal reservoir at temperature T , an analytical expression for the lineshape function can be derived as [S12] D em (ω) = 1 2π dte −i[ω−(ωe−∆re)]t−G(0)+G(t) . (S22) The time-dependent function in the exponent reads T k for each vibrational mode k to simulate the Franck-Condon distribution with a thermal distribution, as long as the factor λ k /ω v k is small, say, less than 0.3. This pseudo temperature is found to be which is obtained by forcing the thermal distribution to have the same vibrational ground state population as the Franck-Condon distribution.
In the main text, the lineshape function and the FGR results are based on Eqs. (S22)-(S24) with the assumption of fast thermalization at T = 300 K. To justify this choice, we compare the decay rates as well as the lineshape function calculated at T = 300 K and pseudo temperature T k in Fig. S2. Excellent agreement for the decay rates can be found. As for the lineshape function, except for the sharp vibronic structures shown in the pseudo temperature case, the main peak positions are also in good agreement.