Tuning the Coherent Propagation of Organic Exciton‐Polaritons through the Cavity Q‐factor

Abstract Transport of excitons in organic materials can be enhanced through polariton formation when the interaction strength between these excitons and the confined light modes of an optical resonator exceeds their decay rates. While the polariton lifetime is determined by the Q(uality)‐factor of the optical resonator, the polariton group velocity is not. Instead, the latter is solely determined by the polariton dispersion. Yet, experiments suggest that the Q‐factor also controls the polariton propagation velocity. To understand this observation, the authors perform molecular dynamics simulations of Rhodamine chromophores strongly coupled to Fabry–Pérot cavities with various Q‐factors. The results suggest that propagation in the aforementioned experiments is initially dominated by ballistic motion of upper polariton states at their group velocities, which leads to a rapid expansion of the wavepacket. Cavity decay in combination with non‐adiabatic population transfer into dark states, rapidly depletes these bright states, causing the wavepacket to contract. However, because population transfer is reversible, propagation continues, but as a diffusion process, at lower velocity. By controlling the lifetime of bright states, the Q‐factor determines the duration of the ballistic phase and the diffusion coefficient in the diffusive regime. Thus, polariton propagation in organic microcavities can be effectively tuned through the Q‐factor.

for both the molecular degrees of freedom, 3 and the cavity mode structure: 4

ĤTC =
Here, σ+ j (σ − j ) is the operator that excites (de-excites) molecule j from the electronic ground (excited) state |S j 0 (R j ) (|S j 1 (R j ) ) into the electronic excited (ground) state |S j 1 (R j ) (|S j 0 (R j ) ); R j is the vector of the Cartesian coordinates of all atoms in molecule j, centered at z j ; âkz (â † kz ) is the annihilation (creation) operator of an excitation of a cavity mode with wave-vector k z ; hν j (R j ) is the excitation energy of molecule j, defined as: with V mol S 0 (R j ) and V mol S 1 (R j ) the adiabatic potential energy surfaces of molecule j in the electronic ground (S 0 ) and excited (S 1 ) state, respectively.
The last term in Equation 1 is the total potential energy of the system in the absolute ground state (i.e., with no excitations in neither the molecules nor the cavity modes), defined as the sum of the ground-state potential energies of all molecules in the cavity.The V mol S 0 (R j ) and V mol S 1 (R j ) adiabatic potential energy surfaces are modelled at the QM/MM level of theory, 5,6 as described in the Computational Details section of the main text.
The third term in Equation 1 models the light-matter interaction within the dipolar approximation through g j (k z ): where µ TDM j (R j ) is the transition dipole moment of molecule j that depends on the molecular geometry (R j ); u cav the unit vector in the direction of the electric component of the cavity vacuum field (i.e., |E| = ω cav (k z )/2 0 V cav ), chosen along the y-direction (see Figure S1); 0 the vacuum permittivity; and V cav the cavity mode volume.
Figure S1: One-dimensional (1D) Fabry-Pérot micro-cavity model. 7Two reflecting mirrors located at − 1 2 x and 1 2 x, confine light modes along this direction, while free propagation along the z direction is possible for plane waves with in-plane momentum k z and energy ω cav (k z ).The vacuum field vector (red) points along the y-axis, reaching a maximum amplitude at x = 0 where the N molecules (magenta ellipses) are placed, distributed along the z-axis at positions z j with 1 ≤ j ≤ N .

Multi-mode cavity model
To discretize the cavity dispersion, we follow Michetti and La Rocca, 7 and impose periodic boundary conditions in the z-direction of the 1D cavity.Under these conditions, the wave vectors, k z , adopt discrete values: k z,p = 2πp/L z with p ∈ Z and L z the length of the 1D cavity.After discretization the molecular Tavis-Cummings Hamiltonian in Equation 1 becomes a (N + n mode ) by (N + n mode ) matrix with four blocks: 4 The elements of this matrix are evaluated in the product basis of adiabatic molecular states times cavity mode excitations: for 1 ≤ j ≤ N , and for N < j ≤ N + n mode .In these expressions |00..0 indicates that all Fock states associated with the n mode cavity modes are empty.The basis state |φ 0 is the ground state of the molecule-cavity system with no excitations in neither the molecules nor cavity modes: The upper left block, H mol , is an N × N matrix that contains the single-photon excitations of the molecules.Because we neglect direct excitonic interactions between molecules, this block is diagonal, with elements labeled by the molecule indices j: for 1 ≤ j ≤ N .Each matrix element of H mol thus represents the potential energy of a molecule, j, in the electronic excited state |S j 1 (R j ) while all other molecules, i = j, are in the electronic ground state |S i 0 (R i ) : The lower right block in Equation 4, H cav , is an n mode × n mode matrix (with n mode = n max − n min + 1) containing the single-photon excitations of the cavity modes, and is also diagonal: for n min ≤ p ≤ n max .Here, â † p excites cavity mode p with wave-vector k z,p = 2πp/L z .In these matrix elements, all molecules are in the electronic ground state (S 0 ).The energy is therefore the sum of the cavity energy at k z,p , and the molecular ground state energies: where, ω cav (k z,p ) is the cavity dispersion (dashed-dotted curve in Figure 1b, main text): with ω 0 the energy at k z,0 = 0, n the refractive index of the medium and c the speed of light in vacuum.
The two N × n mode off-diagonal blocks H int and H int † in the multi-mode Tavis-Cummings Hamiltonian (Equation 4) model the light-matter interactions between the molecules and the cavity modes.These matrix elements are approximated as the inner product between the molecular transition dipole moments on the one hand, and the transverse electric field of the cavity modes at the center of the molecules, on the other hand: with eigenenergies E m .The expansion coefficients β m j and α m p reflect , respectively, the contribution of the molecular excitons (|S j 1 (R j ) ) and the cavity mode excitations (|1 p ) to polariton |ψ m .

Ehrenfest molecular dynamics simulations
MD trajectories of all molecules (including environment) were computed by numerically integrating Newton's equations of motion.The multi-mode Tavis-Cummings Hamiltonian (Equation 4) was diagonalized at each time-step of the simulation to obtain the N + n mode (adiabatic) polaritonic eigenstates |ψ m and energies E m .The total polaritonic wavefunction |Ψ(t) was coherently propagated along with the classical degrees of freedom of the N molecules as a time-dependent superposition of the N + n mode time-independent adiabatic polaritonic states: where c m (t) are the time-dependent expansion coefficients of the time-independent eigenstates, |ψ m , defined in Equation 14.A unitary propagator in the local diabatic basis was used to integrate these coefficients, 9 while the nuclear degrees of freedom of the N molecules evolved on the mean-field potential energy surface:

Cavity decay
Radiative loss through the imperfect cavity mirrors was modeled as a first-order decay process into the overall ground state of the system (i.e., no excitation in neither the molecules nor the cavity modes). 3Assuming that the intrinsic decay rates, γ cav , are the same for all modes, the total loss rate was calculated as the product of γ cav and the total photonic weight, Since , changes in the real and imaginary parts of the (complex) expansion coefficients c m (t) due to spontaneous photonic loss through the mirrors of a low-Q cavity, are: Simultaneously, the population of the zero-excitation subspace, or ground state, ρ 0 (t+∆t), increases as: 2 Further simulation details

Rhodamine model
The Rhodamine molecules, one of which is shown schematically in Figure S2, were modelled with the Amber03 force field, 10 using the parameters provided by Luk et al. 3 After a geometry optimization at the force field level, the molecule was placed at the center of a rectangular box and 3,684 TIP3P water molecules, 11 were added.The simulation box, which contained 11,089 atoms, was equilibrated for 2 ns with harmonic restraints on the heavy atoms of Rhodamine (force constant 1000 kJmol −1 nm −1 ).Subsequently, a 200 ns classical MD trajectory was computed at constant temperature (300 K) using a stochastic dynamics integrator with a friction coefficient of 0.1 ps −1 .
The pressure was kept constant at 1 bar using the Berendsen isotropic pressure coupling algorithm 12   with a time constant of 1 ps.The LINCS algorithm was used to constrain bond lengths, 13 while SETTLE was applied to constrain the internal degrees of freedom of water molecules, 14 enabling a time step of 2 fs in the classical MD simulations.A 1.0 nm cut-off was used for Van der Waals' interactions, which were modeled with Lennard-Jones potentials.Coulomb interactions were computed with the smooth particle mesh Ewald method, 15 using a 1.0 nm real space cut-off and a grid spacing of 0.12 nm.The relative tolerance at the real space cut-off was set to 10 −5 .
Figure S2: Rhodamine model used in our simulations.The QM subsystem, shown in ball-and-stick representation, is described at the HF/3-21G level of theory in the electronic ground state (S 0 ), and at the CIS/3-21G level of theory in the first singlet excited state (S 1 ).The MM subsystem, consisting of the atoms shown in stick representation, and the water molecules (not shown), are modelled with the Amber03 force field.
From the second half of the 200 ns MD trajectory, snapshots were extracted and subjected to further equilibration for 10 ps at the QM/MM level.The time step was reduced to 1 fs.As in previous work, 3,4,16 the fused ring system was included in the QM region and described at the RHF/3-21G level, while the rest of the molecule as well as the water solvent, were modelled with the Amber03 force field 10 and TIP3P water model, 11 respectively (Figure S2).The bond connecting the QM and MM subsystems was replaced by a constraint and the QM part was capped with a hydrogen atom.The force on the cap atom was distributed over the two atoms of the bond. 17The QM system experienced the Coulomb field of all MM atoms within a 1.

Initial conditions
The initial excitation of a polariton wavepacket was created by assigning to the expansion coefficients c m (t = 0) values of a Gaussian distribution centered at k-vector k c = 80 2π Lz = 10.05 µm −1 and covering the whole UP branch (i.e., excluding LP and dark states in the initial wavepacket): 8 where ζ = 10 −14 m 2 is a coefficient characterising the shape of the wavepacket and k m z the expectation value of the in-plane momentum of polariton |ψ m , evaluated as: with k z,p = 2πp/L z the discrete wave vector in a periodic 1D cavity of length L z (see subsection 1.2).

Populations of the lower polariton, upper polariton, and dark states
The time evolution of the populations in the LP, UP, and dark states (plotted in panels g, h, i of

Wavepacket propagation
To monitor the propagation of the wavepackets, we plotted the probability density of the total time-dependent wave function |Ψ(t)| 2 at the positions of the molecules, z j , as a function of time (panels a-f in figure 2 in the main text).We thus represent the density as a discrete distribution at grid points that correspond to the molecular positions, rather than as a continuous distribution.
The probability density of the total time-dependent wave function |Ψ(t)| 2 was calculated as the sum of the probability densities of the molecular |Ψ mol (t)| 2 and photonic |Ψ pho (t)| 2 contributions.
The amplitude of |Ψ mol (t) at position z j in the 1D cavity (with z j = (j − 1)L z /N for 1 ≤ j ≤ N ) was obtained by projecting the excitonic basis state in which molecule j at position z j is excited, onto the total wave function (Equation 15): with β m j the expansion coefficient of the excitonic basis state σ + j |φ 0 in polaritonic state |ψ m (Equation 14), c m (t) the time-dependent expansion coefficients of the total wavefunction |Ψ(t) (Equation 15), and |φ 0 the ground state of the molecule-cavity system with no excitations of neither the molecules nor cavity modes (Equation 7).
The cavity mode excitations are described as plane waves that are delocalized in real space.We therefore obtained the amplitude of the cavity mode excitations in polaritonic eigenstate |ψ m at position z j by Fourier transforming the projection of the cavity mode Fock states, in which cavity mode p is excited, onto |ψ m : where α m p is the expansion coefficient of the cavity mode excitation α † p |φ 0 in polaritonic state |ψ m (Equation 14) and we normalized by 1/ √ N rather than 1/ √ L z , as we represent the density on the grid of molecular positions.The total contribution of the cavity mode excitations to the wavepacket at position z j at time t was then obtained as the weighted sum over the Fourier transforms: with c m (t) the time-dependent expansion coefficient of the adiabatic polaritonic state |ψ m in the total wave function |Ψ(t) (Equation 15).

Transient transmission
Under the assumption that we can neglect reflection, the transmission of light through the system, T , is related to absorbance, A, via the Lambert-Beer law: with ε a the absorption coefficient, C the concentration of absorbers (in our case the Rhodamines), and d the length of the path through which the light passes.Based on the article of Pandya and co-workers, 21 we assume that they probed how the transmission of photons with an energy equal to the excitation from the total ground state |φ 0 (Equation 7, i.e., no excitation in neither the molecules, nor cavity) into the LP branch at 640 nm for their BODIPY-R cavity systems, changes as a function of time (t) and position (here z) after the interaction of the molecule-cavity system with the pump pulse.Assuming furthermore that excitation from the LP or UP into the two-photon manifold is negligible due to absence of resonant transitions at that wavelength, the absorbance of the sample at position z and time t is proportional to the concentration of unexcited molecules: with C 0 the concentration of the Rhodamine molecules in the cavity before interaction with the pump-pulse.We assume that C 0 is uniform and homogeneous.With these approximations, the transient normalized differential transmission in our simulations was calculated as: In analogy to Equation S4 in the Supporting Information of the work by Balasubrahmanyam et al. 22 , the mean squared displacement (MSD) of the transient transmisson signal ∆T /T 0 was calculated as: Here, we treated ε a d as a single parameter between 0 and 1.In Figure S3, we plot the MSD T of the transmitted signal for various values of ε a d.Based on the similarity of the plots, we conclude that the results of the analysis are not very sensitive to the choice of this parameter.

Estimation of the duration and propagation velocity of the ballistic phase
The propagation velocity, υ bal , and duration, τ bal , of the ballistic phase were obtained by fitting the model of Pandya et al. (Equation 5 in their paper 21 ) without σ 0 (which is zero in our simulations) to the MSD T of ∆T (z, t)/T 0 from t = 0 to t = t MSD max T , when the MSD T reaches its maximum: In Figure S4a these fits are plotted as dashed lines.The good agreement between the fit and the MSD T suggests that this function captures the initial MSD T as a function of time for all cavities.
The propagation velocity, υ bal , only weakly depends on the cavity lifetime, τ cav (Figure S4b), because transport in the ballistic regime is primarily governed by the group velocities of UP states.
The duration of the ballistic regime, τ bal , however, depends more strongly on the cavity lifetime (Figure 4b in the main text), suggesting that the increase of propagation distance with cavity lifetime is mainly due to a longer duration of the ballistic phase.28).
The error bars in the plots of υ bal (Figure S4b) and τ bal (Figure 4b in the main text) correspond to the standard deviation σ x of the S = 5 trajectories for each cavity mode lifetime: where x i and x are, respectively, the values of υ gr or τ bal , and their averages, respectively.

Estimation of the diffusion coefficient in the diffusive phase
As explained in the main text, we extracted the diffusion coefficient from the wavepacket MSD.We consider only the part of the wavepacket that is moving slower than the maximum group velocity of the LP (υ max LP , Figure S5) and performed a linear fit to the last 100 fs of the trajectories: These diffusion coefficients are plotted as a function of cavity lifetime in Figure 4c of the main text.Transport of excitons in organic materials can be enhanced through polariton formation when the interaction strength between these excitons and the confined light modes of an optical resonator exceeds their decay rates.While the polariton lifetime is determined by the Q(uality)-factor of the optical resonator, the polariton group velocity is not.Instead, the latter is solely determined by the polariton dispersion.Yet, experiments suggest that the Q-factor also controls the polariton propagation velocity.To understand this observation, we performed molecular dynamics simulations of Rhodamine chromophores strongly coupled to Fabry-Pérot cavities with various Q-factors.Our results suggest that propagation in the aforementioned experiments is initially dominated by ballistic motion of upper polariton states at their group velocities, which leads to a rapid expansion of the wavepacket.Cavity decay in combination with non-adiabatic population transfer into dark states, rapidly depletes these bright states, causing the wavepacket to contract.However, because population transfer is reversible, propagation continues, but as a diffusion process, at lower velocity.By controlling the lifetime of bright states, the Q-factor determines the duration of the ballistic phase and the diffusion coefficient in the diffusive regime.Thus, polariton propagation in organic microcavities can be effectively tuned through the Q-factor.

Introduction
Achieving long-range energy transfer in organic media is a key requirement for enhancing the efficiency of opto-electronic devices, such as organic diodes or solar cells, in which energy transport is limited by the incoherent diffusion mechanism that governs the motion of Frenkel excitons through materials.Recent experiments suggest that strongly coupling such excitons to the confined, but delocalized, modes of an optical resonator (called a cavity in what follows) can enhance transport through hybridization of the molecular excitons with the confined light modes into polaritons [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15].Polaritons are coherent superpositions of molecular and cavity mode excitations that form when the interaction (g) between molecular excitons and cavity modes exceeds their decay rates (κ mol and γ cav , respectively) [16,17,18].The vast majority of these light-matter hybrid states lack a strong contribution from the cavity mode excitations and are therefore "dark".The few remaining states are bright and dispersive owing to their cavity mode contributions (Figure 1).Therefore, these bright states behave as quasi-particles with low effective mass and large group velocity [19], which can be exploited for controlled and long-ranged in-plane energy transport [20].Indeed, in-plane polariton propagation has been observed in a variety of excitonic materials coupled to the confined light modes of Fabry-Pérot cavities [7,15], Bloch Surface Waves [5,11,21], Surface Lattice Resonances [14], and resonances arising from a dielectric constant mismatch between the excitonic medium and the surrounding environment [12].While these observations are in line with theoretical predictions [19,22,20,23,24], the propagation velocity observed in these experiments, is significantly lower than the group velocities inferred from the polariton dispersion (v g = ∂ω/∂k z ).In previous work [25] we used multi-scale molecular dynamics (MD) simulations to resolve this discrepancy, and showed that on long timescales (> 100 fs) polariton propagation is a diffusion process.This diffusion is due to reversible population transfer between the stationary dark state manifold and the highly mobile bright polaritonic states, which renders the propagation speed much slower than the polariton group velocities.While we observed that cavity loss, caused by photon leakage through imperfect mirrors, reduces the distance over which polaritons propagate, we had not systematically investigated the effect of the cavity mode lifetime, τ cav = γ −1 cav , which is related to the quality factor (Q-factor) via Q = ω cav τ cav .The cavity mode lifetime, in combination with the molecular dephasing rate (κ mol ), determines how strong the light-matter interaction (g) needs to be for the molecule-cavity system to enter the strong coupling regime (for which various criteria are commonly employed [18] cav + κ 2 mol )/2; or (iv) g ≥ (γ cav + κ mol )/2).When strong coupling is achieved, the Qfactor only influences the lifetime of organic polaritons, but not the light-matter coupling strength [26].Therefore, the Q-factor neither affects the Rabi splitting between the lower (LP) and upper (UP) polariton branches (Ω Rabi = 2g √ N with N the number of molecules collectively coupled to the confined light modes, Figure 1), nor the group velocity of the bright polariton states.Yet, recent femtosecond transient absorption microscopy (fs-TAM) experiments by Pandya et al. [15] on BODIPY-R dyes in Fabry-Pérot cavities with varying Q-factors, suggest a correlation between the observed polariton velocity and the cavity Q-factor.To address this controversy and determine the effect of the cavity Q-factor on the propagation of organic polaritons, we performed atomistic MD simulations of Rhodamine chromophores strongly coupled to the confined light modes of one-dimensional (1D) uni-directional Fabry-Pérot cavities [27,28] with three different cavity mode lifetimes: τ cav = 15, 30, and 60 fs.As before, the hydrated Rhodamines were modeled at the hybrid Quantum Mechanics / Molecular Mechanics (QM/MM) level [29,30].We calculated mean-field semi-classical MD trajectories of 512 molecules, including their solvent environment, strongly coupled to the 160 confined light modes of a red-detuned cavity (370 meV below the excitation energy of Rhodamine, which is 4.18 eV at the CIS/3-21G//Amber03 level of theory employed here, see Computational Details and Supporting Information (SI) for details).Because in the fs-TAM measurements of Pandya et al. [15] the 10 fs broad-band pump pulses populate mostly UP states, we modeled the initial excitation by preparing a Gaussian wavepacket of UP states centered at hω= 4.41 eV with a bandwidth of σ = 7.07 µm −1 [19].The energy range of the states excited initially in this superposition is indicated by the magenta box in Figure 1b.

Results and Discussion
In Figure 2 we show the time evolution of the probability density of the polaritonic wave function (|Ψ(t)| 2 , Equation 3), after instantaneous excitation of a Gaussian wavepacket of UP states in three Fabry-Pérot microcavities supporting cavity modes with 15, 30, and 60 fs lifetimes, and containing 512 Rhodamine molecules.Animations of the propagation of the total, molecular and photonic wavepackets are provided as SI.In all cavities the wavepacket initially broadens due to the wide range of UP group velocities.Around 30 fs, however, the wavepacket splits into (i) a faster component with a short lifetime that depends on the Q-factor, and (ii) a slower component that is long-lived, but almost stationary.While the lifetime of the slower component is hardly affected by the cavity lifetime, its broadening is Q-factor dependent (Figure 2a-f).The long lifetime of the slower part suggests that it is composed mostly of dark states that lack group velocity, and into which some population of the initially excited UP states has relaxed [31].Nevertheless, due to thermally driven population transfer from these dark states back into propagating polaritons, the slower part still propagates.Because this transfer process is reversible and leads to transient occupation of polaritonic states over a wide range of k z -vectors in both LP and UP branches, propagation occurs in a diffusive manner [25].In contrast, the faster component of the wavepacket is mainly composed of the higher-energy UP states, which have high group velocity.Because the rate at which population transfers from these UP states into the dark state manifold is inversely proportional to the energy gap [32], the main decay channel for these states is radiative emission through the imperfect cavity mirrors.Thus, the lifetime and hence propagation distance of the faster wavepacket component is Q-factor dependent, which is reflected by a faster rise of ground-state population when the cavity mode lifetime decreases (green dashed lines in Figure 2g-i).After the rapid initial expansion of the total wavepacket due to the population in the UP states (blue in Figure 2g-i), transfer into the dark states (black), in combination with irreversible radiative decay from states with the highest group velocity, causes the wavepacket to contract.The extent of this contraction as well as the moment at which it takes place, depends on the cavity mode lifetime, as indicated by both the position, z , and Mean Squared Displacement (MSD) of the wavepackets in Figure 3. Whereas during the expansion phase propagation is dominated by ballistic motion of fast UP states that reach longer distances for higher Q-factors (or equivalently, higher cavity mode lifetimes τ cav ), as indicated by the maximum of the MSD (∼ 68, 23 and 7 µm 2 ), after contraction, propagation continues as diffusion, which is indicated by the linearity of the MSD at the end of the simulations (Figure 3).Diffusion emerges as a consequence of reversible population transfers between stationary dark states and mobile bright states at all k z -vectors in both the UP and LP branches [25].The turnover from ballistic propagation into diffusion is Q-factor dependent and occurs later when the cavity mode lifetime is higher (Figure 3).While simulations provide detailed mechanistic insights into polariton propagation, direct observation of such details is challenging experimentally, in particular because the multiple contributions to a single transient spectral signal of a molecule-cavity system cannot always be unambiguously disentangled [33].In their fs-TAM experiments, Pandya et al. [15] monitored the propagation of the wavepacket, Ψ(z, t), by probing transient changes in cavity transmission at a wavelength that is sensitive to LP absorption.As explained in the SI, to mimic such pump-probe conditions in our simulations, we extracted positiondependent transient changes in the transmission from our trajectories as follows: with ∆T (z, t) = T (z, t)−T 0 the difference between T (z, t), the transmission at position z and time t after excitation, and T 0 = T (z, 0), the transmission before excitation.The variable ε a is the absorption coefficient and d the path length.Because the value of ε a cannot be derived directly from MD simulations, we treated it together with d as a single parameter.Here, we used ε a d= 0.5, but, as we show in SI, varying this parameter does not change the results qualitatively.As was done in experiments [7,15,21], we characterize the propagation of the total wavepacket by the MSD of the transient signal, in our case of the transient transmission (∆T /T 0 , Equation 1): with z 0 the expectation value of the position of the wavepacket at the start of the simulation (t = 0) and the sum is over the positions z i of the N = 512 molecules.Full details of this analysis are provided in SI.
In Figure 4a we plot the MSD T of the transient differential transmission for our cavity systems.As in the experiments (Figure 2c in Pandya et al. [15]), we observe that after a rapid initial increase, the MSD T of the signal decreases.Based on our simulations we attribute this observation to the fast expansion of the wavepacket followed by the contraction.Because two propagation regimes were observed in our simulations, we analysed these regimes separately.In contrast, Pandya et al. assumed a single ballistic phase, and extracted the velocity and duration of that phase from a global fit to the full MSD of the measured ∆T /T 0 signal.Because in the initial stages of the ballistic regime (t < τ cav ) propagation is dominated by the population in UP states with well-defined dispersion, the propagation speed is independent of the Q-factor and determined solely by the UP group velocity (Figure 1c) in all cavities (Figure S2b in SI).However, the duration of this ballistic regime, τ bal , extracted from the ∆T /T 0 signal by fitting the same function as Pandya et al. to the initial rise of the MSD T (SI), depends on the cavity lifetime, and lasts longer if the cavity Q-factor is higher, as shown in Figure 4b.Therefore, as in the MSD plots of the total wavepacket in Figure 3, the MSD T of the ∆T /T 0 signal also reaches the highest value in the cavity with highest Qfactor (or equivalently, the longest cavity mode lifetime τ cav ), in line with the fs-TAM measurements.Whereas in their model Pandya et al. consider only ballistic propagation on a sub-ps timescale, our simulations suggest that also diffusion contributes to propagation on those timescales, when solely the slower part of the wavepacket remains.Therefore, to characterize also this regime, we calculated the diffusion coefficient by fitting the linear regime of the MSD (SI).However, because in the MSD T of the transient transmission (Figure 4a), the linear regime is difficult to discern, we performed the linear fit to the MSD associated with the slower component of the wavepacket at the end of the trajectories (Figure S5 in SI).
In Figure 4c we plot the diffusion coefficients as a function of cavity mode lifetime.Since the Q-factor determines the lifetime of the population in the propagating bright states, the diffusion constant increases with the cavity lifetime.
Because in our simulations we cannot couple as many molecules to the cavity as in experiment (i.e., 10 5 -10 8 molecules [34,35,36]), we overestimate the diffusion coefficient.As we could show previously [32], the rate of population transfer from dark to bright states is inversely proportional to N , whereas the rate in the opposite direction is independent of N .Therefore, the population in the bright propagating states is overestimated when only 512 molecules are coupled to the cavity, leading to a faster diffusion.The overestimation of the diffusion coefficient thus leads to a much more pronounced increase of the wavepacket MSD than in experiment, where the total population residing in the propagating states is significantly lower [21], and diffusion would be hardly observable on sub-ps timescales.Nevertheless, despite these quantitative differences, our simulations provide a qualitative picture that is in line with experimental observations [15].The results of our simulations suggest that the cavity lifetime controls both the duration and length of the initial ballistic phase (Figure 4b) as well as the diffusion constant in the diffusive regime (Figure 4c).Thus, without affecting polariton group velocity, the cavity Q-factor provides an effective means to tune energy transport in the strong coupling regime.

Conclusion
To summarize, we have investigated the effect of the cavity Q-factor on polariton propagation by means of atomistic MD simulations.In line with experiments, we find that the Q-factor determines the propagation velocity and distance of organic polaritons via their lifetimes without affecting group velocities.
Our findings therefore resolve the unexpected correlation between Q-factor and propagation velocity reported by Pandya et al. [15].Our results furthermore underscore that to understand the mechanism of polariton propagation and interpret experiments, it is necessary to include: (i) atomic details for the material; (ii) multiple modes for cavity dispersion; (iii) cavity decay; and (iv) sufficiently many molecules to have dark states providing an exciton reservoir.In particular, treating the molecular degrees of freedom of many molecules is essential for observing wavepacket contraction that is caused by cavity loss in combination with reversible non-adiabatic population transfer between propagating bright states and the stationary long-lived dark state manifold.Our work suggests that an ab initio description of molecules in multi-mode cavities could pave the way to systematically design or optimize polariton-based devices for enhanced energy transport.

Computational Details
We performed mean-field semi-classical [37] MD simulations of 512 Rhodamine chromophores with their solvent environment, strongly coupled to 1D Fabry-Pérot cavities with different radiative lifetimes: τ cav = 15 fs, 30 fs, and 60 fs.To model the interactions between the molecules and the confined light modes of the cavity, we used a Tavis-Cummings Hamiltonian, in which the molecular degrees of freedom are included [38,28].A brief description of our multi-scale cavity MD approach is provided as SI.
In our simulations the Rhodamine molecules were modelled at the QM/MM level, with the QM region containing the fused ring system of the molecule (Figure S2 in SI).The ground-state electronic structure of the QM subsystem was described at the restricted Hartree-Fock (HF) method in combination with the 3-21G basis set [39], while the excited-state electronic structure was modeled with Configuration Interaction, truncated at single electron excitations (CIS/3-21G).The MM region, which contains the rest of the chromophore as well as the solvent (3684 water molecules), was modeled with the Amber03 force field [40] in combination with the TIP3P water model [41].At this level of QM/MM theory, the excitation energy of the Rhodamine molecules is 4.18 eV [38].In previous work, we showed that despite the overestimation of the vertical excitation energy, the topology of the potential energy surfaces is not very sensitive to the level of theory for Rhodamine [31].
The uni-directional 1D cavity with a length of L z = 50 µm, with z indicating the in-plane direction (L x = 163 nm is the distance between the mirrors and x thus indicates the out-of-plane direction, see Figure S1 in the SI), was red-detuned by 370 meV with respect to the molecular excitation energy (4.18 eV at the CIS/3-21G//Amber03 level of theory, dashed line in Figure 1b), such that at wave vector k z = 0, the cavity resonance is hω 0 = 3.81 eV.The cavity dispersion, ω cav (k z ) = ω 2 0 + c 2 k 2 z /n 2 , was modelled with 160 modes (0 ≤ p ≤ 159 for k z = 2πp/L z ), with c the speed of light and n the refractive index.Here, we used n = 1.See SI for further details on the cavity model.The Rhodamine molecules were placed with equal inter-molecular distances on the z-axis of the cavity.To maximize the collective light-matter coupling strength, the transition dipole moments of the Rhodamine molecules were aligned to the vacuum field at the start of the simulation.The same starting coordinates were used for all Rhodamines, but different initial velocities were selected randomly from a Maxwell-Boltzmann distribution at 300 K.With a cavity vacuum field strength of 0.36 MVcm −1 (0.0000707 au), the Rabi splitting, defined as the energy difference between the bright lower (LP) and upper polariton (UP) branches at the wave-vector k res z where the cavity dispersion matches the molecular excitation energy (Figure 1b), is ∼ 325 meV for all cavities (τ cav = 15 fs, 30 fs, and 60 fs).While the choice for a 1D cavity model with only positive k z vectors was motivated by the necessity to keep our simulations computationally tractable, it precludes the observation of elastic scattering events that would change the direction (i.e., in-plane momentum, hk) of propagation.Furthermore, with only positive k z vectors, polariton motion is restricted to the +z direction, but we could show previously [25] that this assumption does not affect the mechanism of the propagation process.Ehrenfest MD trajectories were computed by numerically integrating Newton's equations of motion using a leap-frog algorithm with a 0.1 fs timestep.The multi-mode Tavis-Cummings Hamiltonian (See SI) was diagonalized at each time-step to obtain the N + n mode (adiabatic) polaritonic eigenstates |ψ m and energies E m .The total polaritonic wavefunction |Ψ(t) was coherently propagated along with the classical degrees of freedom of the N molecules as a time-dependent superposition of the N + n mode timeindependent adiabatic polaritonic states: where c m (t) are the time-dependent expansion coefficients of the time-independent polaritonic eigenstates |ψ m (SI).A unitary propagator in the local diabatic basis was used to integrate these coefficients [42], while the nuclear degrees of freedom of the N molecules evolve on the mean-field potential energy surface.Results reported in this work were obtained as averages over five trajectories for each cavity lifetime.For all simulations we used Gromacs 4.5.3[43], in which the multi-mode Tavis-Cummings QM/MM model was implemented [28], in combination with Gaussian16 [44].Further details of the simulations are provided in the SI.

( 13 )
for 1 ≤ j ≤ N and n min ≤ p ≤ n max .Diagonalization of the multi-scale Tavis-Cummings Hamiltonian in Equation 4 yields the N + n mode hybrid light-matter states |ψ m :

Figure 2
Figure 2 in the main text) were obtained by summing over expansion coefficients that belong to LP, UP or dark states: i.e., m |c m (t)| 2 , with m ∈ LP (160 low-energy states), UP (160 higher-energy states), or dark states (the remaining states), respectively.With a Rabi splitting of ∼ 325 meV, the LP and UP branches are sufficiently separated from the dark state manifold for a correct assignment of all states in our simulations.

Figure S3 :
Figure S3: Mean squared displacement (MSD T ) of the ∆T (z, t)/T 0 at time t after instantaneous excitation of a Gaussian wavepacket of UP states in a cavity with 512 molecules and a lifetime of τ cav = 60 fs, plotted for various values of a d.

Figure S4 :
Figure S4: Panel a: MSD T of the transient differential transmission after excitation in the perfect cavity (purple) and in cavities with the decay rates τ cav of 60 fs (red), 30 fs (green), and 15 fs (blue).Panel b: Average polariton group velocity υ gr as a function of the cavity lifetime, plotted as a fraction of the group velocity υ 0gr in the perfect cavity, extracted from a quadratic fit to MSD T s (Equation28).

Figure S5 :
Figure S5: Mean squared displacement of the diffusive part of the total polariton wave function |Ψ(t)| 2 as a function of time in cavities with decay rates τ cav of 60 fs (red), 30 fs (green), and 15 fs (black).The dashed lines are linear fits to the last 100 fs of the MSDs.

Figure 1 :
Figure 1: Panel a: Schematic illustration of an optical Fabry-Pérot micro-cavity filled with Rhodamine chromophores (not to scale).Panel b: Normalised angle-resolved absorption spectrum of the cavity, showing Rabi splitting between the lower polariton (LP, red line) and the upper polariton (UP, blue line) branches.The cavity dispersion and excitation energy of the molecules (4.18 eV at the CIS/3-21G//Amber03 level of theory) are plotted by point-dashed and dashed lines, respectively.The purple frame encloses the range of polaritonic states excited instantaneously by the broad-band pump pulse.Panel c: Group velocity of the LP (red) and UP (blue), defined as ∂ω(k z )/∂k z .9

Figure 2 :
Figure 2: Polariton propagation after resonantly exciting a wavepacket of states in the UP branch centered at z=10 µm.Panels a, b and c: probability density of the total wave function, |Ψ(t)| 2 , as a function of distance (horizontal-axis) and time (vertical-axis) in cavities with different Q-factors (i.e., τ cav = 60, 30 and 15 fs, respectively).Colored arrows in panel a correspond to the time points of the 1D projection in panels d, e and f.The dashed purple and yellow lines indicate propagation at the maximum group velocity of the LP (68 µm ps −1 ) and UP (212 µmps −1 ) branches, respectively.Panels d, e and f: probability density of the total polariton wave function, |Ψ(t)| 2 , at different time points as a function of distance.Panels g, h, and i: populations of the UP (blue), LP (red), and dark (DS, black) states, as well as of the ground state (GS, green dashed line) as functions of time.

Figure 4 :
Figure 4: Panel a: Mean squared displacement of the transmission signal (MSD T ) for different cavities mode lifetimes: τ cav = 60 (red), 30 (green) and 15 fs (black).Circles represent data points for individual runs, while the curves show the averages over all trajectories (five for each τ cav ).Panel b: the duration of the ballistic phase as a function of cavity mode lifetime.Panel c: the diffusion coefficient in the diffusion regime as a function of cavity mode lifetime.