Multi-scale molecular dynamics simulations of enhanced energy transfer in organic molecules under strong coupling

Exciton transport can be enhanced in the strong coupling regime where excitons hybridize with confined light modes to form polaritons. Because polaritons have group velocity, their propagation should be ballistic and long-ranged. However, experiments indicate that organic polaritons propagate in a diffusive manner and more slowly than their group velocity. Here, we resolve this controversy by means of molecular dynamics simulations of Rhodamine molecules in a Fabry-Pérot cavity. Our results suggest that polariton propagation is limited by the cavity lifetime and appears diffusive due to reversible population transfers between polaritonic states that propagate ballistically at their group velocity, and dark states that are stationary. Furthermore, because long-lived dark states transiently trap the excitation, propagation is observed on timescales beyond the intrinsic polariton lifetime. These insights not only help to better understand and interpret experimental observations, but also pave the way towards rational design of molecule-cavity systems for coherent exciton transport.


Introduction
Solar cells based on organic molecules are promising alternatives to the silicon-based technologies that dominate today's market, mostly because organic photovoltaics (OPV) are cheaper to massproduce, lighter, more flexible and easier to dispose of. A key step in light harvesting is transport of excitons from where photons are absorbed to where this energy is needed for initiating a photochemical process [1], usually deeper inside the material of the solar cell. Because excitons r in organic materials are predominantly localized onto single molecules, exciton transport proceeds via incoherent hops [2]. Such random-walk diffusion is, however, too slow to compete with ultra-fast deactivation processes of singlet excitons, such as radiative and non-radiative decay. As exciton diffusion is furthermore hindered by thermal disorder, propagation distances in organic materials typically remain below 10 nm [2]. Such short diffusion lengths limit the efficiency of solar energy harvesting and require complex morphologies of active layers into nanometer sized domains, e.g., bulk heterojunctions in OPVs, which not only complicates device fabrication, but also reduces device stability [3,4].
Distances of hundreds of nanometers have been observed for the diffusion of longer-lived triplet states [5], but because not all organic materials can undergo efficient inter-system crossing or singlet fission, it may be difficult to exploit triplet diffusion in general. Exciton mobility can also be increased through transient exciton delocalization [6][7][8], but as the direct excitonic interactions are weak in most organic materials, molecules need to be ordered to reach this enhanced transport regime.
Alternatively, permanent delocalization over large numbers of molecules can be achieved by strongly coupling the excitons in the material to the confined light modes of optical cavities, such as Fabry-Pérot resonators (Figure 1a) or nano-structured devices [9][10][11]. In this strong light-matter coupling regime the rate of energy exchange between molecular excitons and confined light modes exceeds the intrinsic decay rates of both the excitons and the confined modes leading to the formation of new coherent light-matter states, called polaritons [12][13][14][15][16][17][18][19][20].
The majority of hybrid states in realistic molecule-cavity systems are dark [21][22][23], meaning that they have negligible contributions from the cavity modes. In contrast, the few states with such contributions, are the bright polaritonic states that have dispersion and hence group velocity, defined as the derivative of the polariton energy with respect to in-plane momentum (i.e., k z in Figure 1b). In the out-of-plane cavity direction (i.e., perpendicular to the mirrors), these states are delocalised over the molecules inside the mode volume, while in the in-plane direction (i.e., parallel to the mirrors) they behave as quasi-particles with a low effective mass and large group velocity (i.e., fractions of the speed of light). These polaritonic properties can be exploited for both out-of-plane [9][10][11][24][25][26][27][28][29][30][31][32][33], and in-plane energy transport .
Indeed, at cryogenic temperatures in-plane ballistic propagation at the group velocity of polaritons was observed for polariton wavepackets in a Fabry-Pérot microcavity containing an In 0.05 Ga 0.95 As quantum well [34]. Ballistic propagation was also observed for polaritons formed between organic molecules and Bloch surface waves [38,42,51], while a combination of ballistic transport on an ultrashort timescale (sub-50 fs) and diffusive motion on longer timescales was observed for cavity-free polaritons [43], for which strong coupling was achieved through a mismatch of the refractive indices between thin layers of densely-packed organic molecules and a host material [55]. In contrast, experiments on strongly coupled organic J-aggregates in metallic micro-cavities suggest that molecular polaritons propagate in a diffusive manner and much more slowly than their group velocities [40]. Furthermore, despite a low cavity lifetime in the order of tens of femtoseconds in these experiments, propagation was observed over several picoseconds, which was attributed to a long lifetime of the lower polariton (LP) [17,40].
To address these controversies and acquire atomistic insights into polariton propagation, we performed multi-scale molecular dynamics (MD) simulations [56, 57] of solvated Rhodamine molecules strongly coupled to the confined light modes of a one-dimensional (1D) Fabry-Pérot microcavity, shown in Figure 1a [58]. As in previous work [56], the electronic ground state (S 0 ) of the molecules was modeled at the hybrid Quantum Mechanics / Molecular Mechanics (QM/MM) level [59], using the restricted Hartree-Fock (HF) method for the QM subsystem, which contains the fused rings, in combination with the 3-21G basis set [60]. The MM subsystem, consisting of the rest of the Rhodamine molecule and the water, was modeled with the Amber03 force field [61]. The first electronic excited state (S 1 ) of the QM region was modelled with Configuration Interaction, truncated at single electron excitations (CIS/3-21G//Amber03). At this level of theory, the excitation energy of Rhodamine is 4.18 eV, which is significantly overestimated with respect to experiments. This discrepancy is due to the limited size of the basis set and the neglect of electron-electron correlation in the ab initio methods. While including electron-electron correlation into the description of the QM region improves the vertical excitation energy, we show in the Supporting Information (SI) that this does not significantly change the topology of the relevant potential energy surfaces, which determines the molecular dynamics ( Figure S3).
We computed semi-classical Ehrenfest [62] MD trajectories of 1024 Rhodamine molecules inside a 1D cavity of length 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). The cavity 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 ℏω 0 = 3.81 eV. The dispersion of this cavity, ω 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) [35]. With a cavity vacuum field strength of 0.26 MVcm −1 , 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), was ∼ 325 meV. 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, ℏk) of propagation. Furthermore, with only positive k z vectors, polariton motion is restricted to the +z direction, but we show in the SI ( Figure S14) that this assumption does not affect our conclusions about the transport mechanism.
Newton's equations of motion were integrated numerically with a 0.1 fs time step using forces derived on-the-fly from the mean-field potential energy surface provided by the total time-dependent polaritonic wave function, |Ψ(t)⟩, which was expanded in the basis of the timeindependent adiabatic eigenstates of the cavitymolecule Hamiltonian (SI) [35,57,58,63]. The total wavefunction was evolved along with the r classical MD trajectory by unitary propagation in the local diabatic basis [64]. A complete description of the methods employed in this work, including details of the multi-scale MD model for strongly coupled molecules, is provided as SI.
In these experiments the strongly-coupled systems were excited either resonantly into the bright polaritonic states [43,47], or off-resonantly into an uncoupled molecular electronic state [38,40,42,46,51]. To understand the effect of the excitation on polariton-mediated transport, we performed simulations for both initial conditions.
Resonant excitation into the LP branch by a short broad-band laser pulse, typically used in time-resolved experiments [34,43,47] was modeled by preparing a Gaussian wavepacket of LP states centered at ℏω = 3.94 eV where the group velocity of the LP branch, defined as v LP gr (k z ) = ∂ω LP (k z )/∂k z , is highest, and with a bandwidth of σ = 0.707 µm −1 [35]. Off-resonant excitation in a molecule-cavity system is usually achieved by optically pumping a higher-energy electronic state of the molecules [38,40,42,51], which then rapidly relaxes into the lowest energy excited state (S 1 ) according to Kasha's rule [65]. We therefore modelled off-resonant photo-excitation by starting the simulations directly in the S 1 state of a single molecule, located at z = 5 µm in the cavity (SI). We assume that the intensity of the excitation pulse in both cases is sufficiently weak for the system to remain within the single-excitation subspace. We thus exclude multi-photon absorption and model the interaction with the pump pulse as an instantaneous absorption of a single photon.
Because the light-confining structures used in previous experiments (e.g., Fabry-Pérot cavities [34,40,47,48], Bloch surface waves [38,42,51], or plasmonic lattices [41,46,54]) span a wide range of quality factors (Q-factors), we also investigated the effect of the cavity mode lifetime on the transport by performing simulations in an ideal lossless cavity with no photon decay (i.e., γ cav = 0 ps −1 ), and a lossy cavity with decay rate of 66.7 ps −1 . This decay rate corresponds to a lifetime of 15 fs, which is in the same order of magnitude as the 2 -15 fs lifetimes reported for metallic Fabry-Pérot cavities in experiments [40,[66][67][68]. In addition to cavity loss, also internal conversion via the conical intersection seam between the S 1 and S 0 potential energy surfaces [69], can provide a decay channel for the excitation. However, because in our Rhodamine model, the minimum energy conical intersection is 1.3 eV higher in energy than the vertical excitation (SI), and is therefore unlikely to be reached on the timescale of our simulations, we neglect internal conversion processes altogether.

Resonant excitation
First, we explore how polaritons propagate after resonant excitation of a Gaussian wavepacket of LP states with a broad-band laser pulse. In Figure 2, we show the time evolution of the probability density of the polaritonic wave function, |Ψ(t)| 2 after such excitation in both a perfect lossless cavity with an infinite Q-factor (γ cav = 0 ps −1 , top panels) and a lossy cavity with a low Q-factor (γ cav = 66.7 ps −1 , bottom panels) containing 1024 Rhodamine molecules. Plots of wavepacket propagation in systems with 256 and 512 molecules are provided as SI ( Figures S4-S5), as well as animations of the wavepackets for all system sizes.

Lossless cavity
In the perfect lossless cavity the total wavepacket |Ψ(t)| 2 initially propagates ballistically close to the maximum group velocity of the LP branch (v LP,max gr = 68 µmps −1 , Figure 1c), until around 100 fs (see animations in the SI), when it slows down as evidenced by a decrease in the slope of the expectation value of the position of the wavepacket ⟨z⟩ in Figure 3a. The change from a quadratic to a linear time-dependence of the Mean Square Displacement (Figure 3c) at t = 100 fs furthermore suggests a transition from ballistic to diffusive motion.
During propagation, the wavepacket broadens and sharp features appear, visible as vertical lines in both the total and molecular wavepackets in Fig. 2 Polariton propagation after on-resonant excitation of a wavepacket in the LP branch centered at z = 5 µm. Panels a, b and c: total probability density |Ψ(t)| 2 , probability density of the molecular excitons |Ψ exc (t)| 2 and of the cavity mode excitations |Ψ pho (t)| 2 , respectively, as a function of distance (horizontal axis) and time (vertical axis), in a cavity with perfect mirrors (i.e., γ cav = 0 ps −1 ). The magenta dashed line indicates propagation at the maximum group velocity of the LP (68 µmps −1 ). Panel d: Contributions of molecular excitons (black) and cavity mode excitations (red) to |Ψ(t)| 2 as a function of time in the perfect cavity. Without cavity losses, no ground state population (blue) can build up. Panels e, f, and g: |Ψ(t)| 2 , |Ψ exc (t)| 2 and |Ψ pho (t)| 2 , respectively, as a function of distance (horizontal axis) and time (vertical axis), in a lossy cavity (γ cav = 66.7 ps −1 ). Panel h: Contributions of the molecular excitons (black) and cavity mode excitations (red) to |Ψ(t)| 2 as a function of time. The population in the ground state, created by radiative decay through the imperfect mirrors, is plotted in blue. Figure 2a-b and as peaks in the wavepacket animations provided as SI. These peaks coincide with the z positions of molecules that contribute to the wavepacket with their excitations during propagation. Such peaks are not observed if there is no disorder and the molecular degrees of freedom are frozen ( Figure S15), but appear already at the start of the simulation when the initial configurations of the molecules are all different ( Figure S21). Similar observations were made by Agranovich and Gartstein [35], who attributed these peaks to energetic disorder among the molecular excitons. We therefore also assign these peaks to a partial localization of the wavepacket at the molecules due to structural disorder that alters their contribution to the wavepacket. In contrast, because the cavity modes are delocalized in space, the photonic wavepacket remains smooth throughout the propagation (Figure 2c).
The transition from ballistic propagation to diffusion around 100 fs coincides with the onset of the molecular excitons dominating the polaritonic wavepacket, as shown in Figure 2d, in which we plot the contributions of the molecular excitons (black line) and cavity mode excitations (red line) to the total wave function (see SI for details of this analysis). Because in the perfect cavity, photon leakage through the mirrors is absent (i.e., γ cav = 0 ps −1 ), the decrease of cavity mode excitations is due to population transfer from bright LP states into the dark state manifold ( Figure S19b) [70][71][72]. Thus, while resonant excitation of LP states initially leads to ballistic motion with the central group velocity of the wavepacket, as evidenced by the quadratic dependence of the Mean Squared Displacement on time (Figure 3c), population transfer into dark states turns the propagation into a diffusion process, r as evidenced by a linear time-dependence of the Mean Square Displacement after ∼100 fs.
Since dark states lack group velocity, and are therefore stationary, while excitonic couplings between molecules are neglected in our model (see SI), propagation in the diffusive regime must still involve bright polariton states. Our simulations therefore suggest that while, initially, molecular vibrations drive population transfer from the propagating bright states into the stationary dark states [73], this process is reversible, causing new wavepackets to form continuously within the full range of LP group velocities. Likewise, the propagation of transiently occupied bright states is continuously interrupted by transfers into dark states, and re-started with different group velocities. This re-spawning process leads to the diffusive propagation of the excitation observed in Figure 2, with an increasing wavepacket width (Figure 3a), in line with experimental observations [40,42,43,51]. Top panels: Expectation value of the position of the total-time dependent wavefunction ⟨Ψ(t)|ẑ|Ψ(t)⟩/⟨Ψ(t)|Ψ(t)⟩ after on-resonant excitation in an ideal cavity (a) and a lossy cavity (b). The black lines represent ⟨z⟩ while the shaded area around the lines represents the root mean squared deviation (RMSD, i.e., ⟨(z(t) − ⟨z(t)⟩) 2 ⟩). Bottom panels: Mean square displacement (MSD, i.e., ⟨(z(t) − ⟨z(0)⟩) 2 ⟩) in the ideal lossless cavity (c) and the lossy cavity (d). Magenta lines are quadratic fits to the MSD and cyan lines are linear fits.

Lossy cavity
Including a competing radiative decay channel by adding photon losses through the cavity mirrors at a rate of γ cav = 66.7 ps −1 , leads to a rapid depletion of the polariton population (Figure 2h), but does not affect the overall transport mechanism: the wavepacket still propagates in two phases, with a fast ballistic regime followed by slower diffusion. However, in contrast to the propagation in the ideal lossless cavity, we observe that the wavepacket temporarily contracts. This contraction is visible as a reduction of both the expectation value of ⟨z⟩ and the Mean Squared Displacement between 60 to 130 fs in the right panels of Figure 3.
Initially the propagation of the wavepacket is dominated by ballistic motion of the population in the bright polaritonic states moving at the maximum group velocity of the LP branch. However, due to non-adiabatic coupling [73], some of that population is transferred into dark states that are stationary. Because non-adiabatic population transfer is reversible, the wavepacket propagation undergoes a transition into a diffusion regime, which is significantly slower, as also observed in the ideal cavity ( Figure 3c).
In addition to these non-adiabatic transitions, radiative decay further depletes population from the propagating bright polaritonic states. Because before decay, this population has moved much further than the population that got trapped in the dark states, the expectation value of ⟨z⟩, as well as the Mean Square displacement, which were dominated initially by the fast-moving population, decrease until the slower diffusion process catches up and reaches the same distance around 130 fs (right panels in Figure 3). Such contraction of the wavepacket in a lossy cavity is consistent with the measurements of Musser and co-workers, who also observe such contraction after on-resonant excitation of UP states [47].
Because of the contraction, it is difficult to see where the transition between ballistic and diffusion regimes occurs in Figure 3d. We therefore instead extrapolated the linear regime, and estimate the turn-over at 30 fs, where the quadratic fit to the ballistic regime intersects the extrapolated fit to the diffusion regime. As in the perfect lossless cavity, the transition between ballistic and diffusion regimes occurs when the population of molecular excitons exceeds the population of cavity mode excitations ( Figure 2h). However, due to the radiative decay of the latter, this turnover already happens around 30 fs in the lossy cavity simulations.
Owing to the short cavity mode lifetime (15 fs), most of the excitation has already decayed into the ground state at 100 fs, with a small remainder "surviving" in dark states (Figure 2h) that lack mobility. Because cavity losses restrict the lifetime of bright LP states, the distance a wavepacket can reach is limited due to (i) the shortening of the ballistic phase, and (ii) the reduction of the diffusion coefficient (i.e., the slope of ⟨z⟩, Figure 3c) in the second phase. Therefore, the overall velocity is significantly lower than in the perfect cavity, suggesting a connection between cavity Q-factor and propagation velocity [47], while also the broadening of the wavepacket is reduced (Figure 3b). Furthermore, because the rate of population transfer is inversely proportional to the energy gap [73], and hence highest when the LP and dark states overlap [71], we speculate that the turn-over between the ballistic and diffusion regimes depends on the overlap between the absorption line width of the molecules and the polaritonic branches, and can hence be controlled by tuning the excitation energy to move the center of the initial polaritonic wavepacket along the LP branch. In addition, the direction of ballistic propagation can be controlled by varying the incidence angle of the on-resonant excitation pulse.

Comparison to experiments
Our observations are in line with transient microscopy experiments, in which broad-band excitation pulses were used to initiate polariton propagation. At low temperatures Freixanet et al. observed ballistic wavepacket propagation for a strongly coupled quantum dot [34]. If we suppress vibrations that drive population transfer by freezing the nuclear degrees of freedom, we also observe such purely ballistic motion ( Figure S15). In contrast, in room temperature experiments on cavityfree molecular polaritons, Pandya et al. identified two transport regimes: a short ballistic phase followed by diffusion [43]. Based on the results of our simulations, we attribute the first phase to purely ballistic wavepacket propagation of photoexcited LP states. The slow-down of the transport in the second phase is attributed to reversible trapping of population inside the stationary dark state manifold. Owing to the reversible transfer of population between these dark states and the LP states, propagation continues diffusively at time scales exceeding the polariton lifetime, in line with experiment [40,43].

Off-resonant excitation
Next, we investigate polariton propagation after an off-resonant excitation into the S 1 electronic state of a single Rhodamine molecule, located at z = 5 µm. In Figure 4 we show the time evolution of the probability density of the total polaritonic wave function, |Ψ(t)| 2 , after such excitation in both a perfect lossless cavity with an infinite Qfactor (γ cav = 0 ps −1 , top panels) and a lossy cavity with a low Q-factor (γ cav = 66.7 ps −1 , bottom panels) containing 1024 Rhodamine molecules. Plots of the wavepacket propagation in systems with 256 and 512 molecules are provided as SI ( Figures S8-S9), as well as animations of the wavepackets for all system sizes.

Lossless cavity
In the lossless cavity with perfect mirrors, the excitation, initially localised at a single molecule, rapidly spreads to other molecules (see animation in the SI). In contrast to the ballistic movement observed for on-resonant excitation, the wavepacket spreads out instead, with the front of the wavepacket propagating at a velocity that closely matches the maximum group velocity of the LP branch (68 µmps −1 , Figure 1c), while the expectation value of the wavepacket position (⟨z⟩, Figure 5a) moves at a lower pace (∼10 µmps −1 ).
Because we do not include negative k z -vectors in our cavity model, propagation can only occur in the positive z direction. With negative k z -vectors, propagation in the opposite direction cancels such motion leading to ⟨z⟩ ≈ 0 ( Figure S14a). Nevertheless, since the Mean Square Displacement is not affected by breaking the symmetry of the 1D cavity, and increases linearly with time in both uni-and bi-directional cavities (Figures 5c and S14b), we consider it reasonable to assume that the mechanism underlying the propagation process is identical.
Because the population of dark states dominates throughout these simulations (Figure 4d), Fig. 4 Polariton propagation after off-resonant excitation of a single molecule located at z = 5 µm into S 1 . Panels a, b and c: total probability density |Ψ(t)| 2 , probability density of the molecular excitons |Ψ exc (t)| 2 and of the cavity mode excitations |Ψ pho (t)| 2 , respectively, as a function of distance (horizontal axis) and time (vertical axis), in a cavity with perfect mirrors (i.e., γ cav = 0 ps −1 ). The magenta and yellow dashed lines indicate propagation at the maximum group velocity of the LP (68 µmps −1 ) and UP (212 µmps −1 ), respectively. Panel d: Contributions of the molecular excitons (black) and cavity mode excitations (red) to |Ψ(t)| 2 as a function of time in the perfect cavity. Without cavity decay, there is no build-up of ground state population (blue). Panels e, f, g: |Ψ(t)| 2 , |Ψ exc (t)| 2 and |Ψ pho (t)| 2 , respectively, as a function of distance (horizontal axis) and time (vertical axis), in a lossy cavity (i.e., γ cav = 66.7 ps −1 ). Panel h: Contributions of the molecular excitons (black), and cavity mode excitations (red) to |Ψ(t)| 2 as a function of time in the lossy cavity. The population in the ground state, created by radiative decay through the imperfect mirrors, is plotted in blue. and direct excitonic couplings are not accounted for in our model (SI), the observed propagation must again involve bright polariton states. Since the initial state, with one molecule excited, is not an eigenstate of the molecule-cavity system, population exchange from this state into the propagating bright states is not only due to displacements along vibrational modes that are overlapping with the non-adiabatic coupling vector [73], but also due to Rabi oscillations, in particular at the start of the simulation.
To quantify to what extent the overall propagation is driven by population transfers due to the molecular displacements, we performed additional simulations at 0 K with all nuclear degrees of freedom frozen. As shown in Figure S17, the propagation is reduced at 0 K, and the wavepacket remains more localized on the molecule that was initially excited, than at 300 K. A quadratic timedependence of the Mean Square Displacement of the cavity mode contributions to the wavepacket ( Figure S18f) furthermore suggest that the mobility at 0 K is driven by the constructive and destructive interferences of the bright polaritonic states, which evolve with different phases (i.e., e −iE m t/ℏ ).
The reduced mobility of the wavepacket at 0 K compared to 300 K ( Figure S18) confirms that thermally activated displacements of nuclear coordinates, which are absent at 0 K, are essential to drive population into the bright states and sustain the propagation of the polariton wavepacket. Thus, as during the diffusion phase observed for on-resonant excitation, ballistic motion of bright states is continuously interrupted and restarted with different group velocities, which makes the overall propagation appear diffusive with a Mean Square Displacement that depends linearly on time (Figure 5c), in line with experimental observations [40,42].
In the perfect cavity, propagation and broadening continue indefinitely due to the long-range ballistic motion of states with higher group velocities. Indeed, a small fraction at the front of the wavepacket, which moves even faster than the maximum group velocity of the LP (indicated by a magenta dashed line in Figure 4), is mostly composed of higher-energy UP states. These states not only have the highest in-plane momenta, but also decay most slowly into the dark state manifod of the perfect cavity due to the inverse dependence of the non-adiabatic coupling on the energy gap [57]. Momentum-resolved photo-luminenscence spectra at two distances from the initial excitation spot ( Figure S23, SI) confirm that the front of the wavepacket is indeed composed of UP states: at short distances (z = 10µm) from the excitation spot (z = 5 µm), the emission spectrum, accumulated over 100 fs simulation time, closely matches the full polariton dispersion of Figure 1b, displaying both the LP and UP branches. In contrast, further away from the excitation spot (z = 20 µm), the emission exclusively originates from the higher energy UP states, suggesting that only these states can reach the longer distance within 100 fs.

Lossy cavity
Adding a radiative decay channel for the cavity mode excitations (γ cav = 66.7 ps −1 ) restricts the distance over which polaritons propagate (Figure 4e-g), but does not affect the overall transport mechanism, as we also observe a linear increase of the MSD with time ( Figure 5d). While the propagation in the lossy cavity initially is very similar to that in the ideal lossless cavity, radiative decay selectively depletes population from the propagating bright states and the wavepacket slows down, as evidenced by the expectation value of the displacement, ⟨z⟩, levelling off in Figure 5b. In addition, since the maximum distance a wavepacket can travel in lossy cavities is determined by the cavity lifetime in combination with the group velocity [38], the broadening of the wavepacket is also more limited when cavity losses are included ( Figure 5b). Furthermore, even if the dark states do not have a significant contribution from the cavity mode excitations, the reversible transfer of population between the dark state manifold and the decaying bright polaritonic states, also leads to a significant reduction of dark state population in the lossy cavity as compared to the ideal lossless cavity (Figure 4d and f). Nevertheless, dark states still provide "protection" from cavity losses as the overall lifetime of the photo-excited moleculecavity system (> 150 fs) significantly exceeds that of the cavity modes (15 fs).

Comparison to experiments
In microscopy experiments relying on off-resonant optical pumping, polariton emission is typically observed between the excitation spot and a point several microns further away [38,40,42,46,51]. While such broad emission pattern is reminiscent of a diffusion process, the match between total distance over which that emission is detected on the one hand and the product of the maximum LP group velocity and cavity lifetime on the other hand, suggest ballistic propagation. The results of our simulations are thus in qualitative agreement with such observations as also our results suggest that, while polariton propagation appears diffusive under off-resonant excitation conditions, the front of the wave packet propagates close to the maximum group velocity of the LP branch.
Based on the analysis of our MD trajectories we propose that on the experimentally accessible timescales, polariton propagation appears diffusive due to reversible population transfers between stationary dark states and propagating bright states. For lossy cavities, radiative decay of the cavity modes further slows down polariton transport such that the excitation reaches a maximum distance before decaying completely. Because a large fraction of the population resides in the nondecaying dark states, the lifetime of the moleculecavity system is extended [71], and polariton propagation can be observed on timescales far beyond the cavity lifetime, in line with experiment [40].
Note that in our simulations we couple excitons only to the modes of the Fabry-Pérot cavity, whereas in experiments with micro-cavities constituted by metal mirrors, excitons can in principle also couple to the surface plasmon polaritons (SPP) below the light line that are supported by these metal surfaces. While their role will depend on the details of the set-up (e.g., the materials used, energy of the relevant molecular excitations, etc.), we cannot rule out that reversible population transfer between the dark states and SPP-exciton polaritons also contributes to the effective diffusion constant observed in those experiments [40,48]. However, because the SPP decays exponentially away from the metal surface, and SPP-exciton polaritons also have group velocity, the qualitative behavior is not expected to change.

Size dependence
Due to limitations on hard-and software, the number of molecules that we can include in our simulations is much smaller than in experiments [74-76]. We therefore investigated how the number of molecules, N , coupled to the cavity affects the propagation by performing simulations for different N . To keep the Rabi splitting (∼ 325 meV) constant, and hence polariton dispersion the same, we scaled the cavity mode volume with the number of molecules N (see SI for details).
While the transport mechanism is not strongly affected by N ( Figures S4-S12), the total population that resides in the bright states decreases when the number of molecules, and hence the number of dark states, increases, in particular in the diffusion phase. Such decrease in bright state population is due to the 1/N scaling of the rate at which population transfers between dark and bright states [73]. Because the number of dark states is proportional to N , while the number of bright states is constant for a fixed number of cavity modes, this dependency affects the ratio between the population in the dark and bright states, with the latter rapidly decreasing with increasing N . As the overall propagation velocity is determined by the population in bright states, also the velocity is inversely proportional to N ( Figure S13). Therefore, in experiments, with 10 5 -10 8 molecules inside the mode volume [74-76], the propagation velocity is much lower than in our simulations.
Nevertheless, because of the 1/N scaling, the effective polariton propagation velocity approaches the lower "experimental limit" of 10 5 coupled molecules [74] already around 1000 molecules. We therefore consider the results of the simulations with 1024 Rhodamines sufficiently representative for experiment and for providing qualitative insights into polariton propagation. Indeed, a propagation of 9.6 µmps −1 in the cavity containing 1024 molecules is about an order of magnitude below the maximum group velocity of the LP (68 µmps −1 ) in line with experiments on organic microcavities [40], and cavity-free polaritons [43]

Conclusions
We have investigated exciton transport in Rhodamine cavities by means of atomistic MD simulations that not only include the details of the cavity mode structure [57, 58], but also the chemical details of the material [56]. The results of our simulations suggest that the transport is driven by an interplay between propagating bright polaritonic states and stationary dark states. Reversible population exchanges between these states interrupt ballistic motion in bright states and make the overall propagation process appear diffusive. While for off-resonant excitation of the molecule-cavity system, these exchanges are essential to transfer population from the initially excited molecule into the bright polaritonic branches and start the propagation process, the exchanges limit the duration of the initial ballistic phase for on-resonant excitation. As radiative decay of the cavity modes selectively depletes the population in bright states, ballistic propagation is restricted even further if the cavity is lossy. Because dark states lack in-plane momentum, the reversible population exchange between dark and bright states causes diffusion in all directions. Therefore, under off-resonant excitation conditions, the propagation direction cannot be controlled. In contrast, because bright states carry momentum, the propagation direction in the ballistic phase can be controlled precisely by tuning the incidence angle and excitation wavelength under on-resonant excitation conditions.
The rate at which population transfers between bright and dark states depends on the non-adiabatic coupling vector, whose direction and magnitude are determined by the Huang-Rhys factor in combination with the frequency of the Franck-Condon active vibrations [73], both of which are related to the molecular Stokes shift [77]. In addition, because the non-adiabatic coupling is inversely proportional to the energy gap [73], the Stokes shift in combination with the Rabi splitting, also determines the region on the LP branch into which population transfers after off-resonant excitation of a single molecule [78][79][80][81]. We therefore speculate that the Stokes shift can be an important "control knob" for tuning the coherent propagation of polaritons. Because our Rhodamine model features the key photophysical characteristics of an organic dye molecule, we speculate that the propagation mechanism observed in our simulations is generally valid for exciton transport in strongly-coupled organic micro-cavities, in which the absorption line width of the material exceeds the Rabi splitting and there is a significant overlap between bright and dark states. To confirm this, we have also performed simulations of exciton transport in cavities containing Tetracene and Methylene Blue and observed that the propagation mechanism remains the same (SI).
Future work will be aimed at investigating how the propagation can be controlled by tuning molecular parameters, temperature, Rabi splitting, or cavity Q-factor [82]. Because we include the structural details of both cavity and molecules, our simulations, which are in qualitative agreement with experiments, pave the way to systematically optimize molecule-cavity systems for enhancing exciton energy transfer.
Supplementary information. The Supporting Information contains: (i) details of the multiscale MD simulation model; (ii) details of simulation setups and analysis; (iii) results of simulations with different numbers of molecules, at 0 K, in symmetric cavities, with positional and energetic disorder, with different vacuum fields, as well as with Tetracene and Methylene Blue instead of Rhodamine; and (iv) animations of all wavepackets.

Declarations
Here,σ + j (σ − j ) is the operator that excites (de-excites) molecule j from the electronic ground ; 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 (PESs) 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 PESs are modelled at the hybrid quantum mechanics / molecular mechanics (QM/MM) level of theory. 6,7 Figure S1: One-dimensional (1D) Fabry-Pérot micro-cavity model. Two 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 .
The third term in Equation 1 describes 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 cavity vacuum field (i.e., |E| = ℏω cav (k z )/2ϵ 0 V cav ), here along the y-direction ( Figure S1); ϵ 0 the vacuum permittivity; and V cav the cavity mode volume.

One-dimensional Periodic Cavity
Following Michetti and La Rocca, 4 we impose periodic boundary conditions in the z-direction of our cavity, and thus restrict the wave vectors, k z , to discrete values: k z,p = 2πp/L z with p ∈ Z and L z the length of the 1D cavity. With this approximation the molecular Tavis-Cummings Hamiltonian  in Equation 1 can be represented as an (N + n mode ) by (N + n mode ) matrix with four blocks: 5 We compute the elements of this matrix 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 the Fock states for all n mode cavity modes are empty. The basis state |ϕ 0 ⟩ is the ground state of the molecule-cavity system with no excitations of neither the molecules nor cavity modes: The upper left block, H mol , is an N × N matrix containing 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, 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 curve in Figure 1b, main text): with ℏω 0 the energy at k 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 for 1 ≤ j ≤ N and n min ≤ p ≤ n max .
Diagonalization of the multi-scale Tavis with eigenenergies E m . The β m j and α m p expansion coefficients reflect the contribution of the molecular excitons (|S j 1 (R j )⟩) and the cavity mode excitations (|1 p ⟩) to polariton |ψ m ⟩.

Mean-Field Non-Adiabatic Molecular Dynamics Simulations
where c m (t) are the time-dependent expansion coefficients of the time-independent polaritonic basis functions |ψ m ⟩ defined in Equation 14. A unitary propagator in the local diabatic basis is used to integrate these coefficients, 9 while the nuclear degrees of freedom of the N molecules evolve on the mean-field PES: Population transfers between polaritonic states |ψ m ⟩ and |ψ l ⟩ are driven by the non-adiabatic coupling, or vibronic coupling V NAC m,l . 10-12 V NAC m,l depends on the overlap between the non-adiabatic coupling vector d m.l and the velocities of the atomsṘ: 13 whereṘ is a vector containing the velocities of all atoms in the system, and the non-adiabatic coupling vector d m,l is given by: where E m and E l are the eigenvalues associated with |ψ m ⟩ and |ψ l ⟩, respectively. The direction of d m,l determines the direction of population transfer, while its magnitude determines the amount of population that is transferred.

Cavity Decay
Radiative loss through the cavity mirrors is 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). 14 Assuming that the intrinsic decay rates γ cav are the same for all modes, the total loss rate is calculated as the product of γ cav and the total photonic weight, n mode p |α m p | 2 , of state |ψ m ⟩. Thus, after an MD step ∆t, ρ m (t) = |c m (t)| 2 becomes: Since ρ m = (ℜ[c m ]) 2 +(ℑ[c m ]) 2 , changes in the real and imaginary parts of the (complex) expansion coefficients c m (t), due to spontaneous photonic loss in a low-Q cavity, are: Simultaneously, the population of the zero-excitation subspace, or ground state, ρ 0 (t+∆t), increases as: 2 Simulation details

Rhodamine Model
Rhodamine molecules, one of which is shown schematically in Figure S2a, were modelled with the Amber03 force field, 15 using the parameters provided by Luk et al. 14 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, 16  The level of theory used in this work is a necessary compromise between efficiency and accuracy.
Because no experimental data is available for the Rhodamine model used here, we instead verified the validity of our model by recomputing the electronic energies along an excited-state molecular dynamics trajectory of a single Rhodamine molecule in water at higher levels of theory: • Time-dependent density functional theory (TD-DFT) 25 with the CAM-B3LYP functional, 26,27 in combination with the 6-31G(d) basis set; 28 • Complete Active Space Self Consistent Field (CASSCF) with 12 electrons in 12 π orbitals, 29 averaged over the first two states (SA2) and expanded in a cc-pVDZ basis set; 30 • Extended Multi-Configuration Quasi-Degenerate Perturbation Theory to second order (xMC-QDPT2), 31 based on the SA2-CASSCF(12,12)/cc-pVDZ reference.
The TD-DFT calculations were carried out with the Gaussian09 program, 32  Note that because the xMCQDPT2 method only corrects the energy, but not the wave function, the transition dipole moment is the same as that of the underlying CASSCF reference. Panel c: the energy gaps between the excited and ground state.
As shown in Figure S3, the absolute excitation energy is overestimated at the CIS/3-21G level of theory with respect to the more accurate xMCQDPT2//SA2-CASSCF(12,12)/cc-pVDZ approach.
However, the positions of the minima and maxima coincide at all levels of theory. Therefore, we conclude that the topology of the PES is not very sensitive to the level of theory and that the nuclear dynamics will be sufficiently similar at all levels of theory.
Furthermore, the Stokes shift, which was shown to have an important effect on polariton relaxation, 34 was evaluated by first optimizing the molecule in the ground and excited state and then taking the difference between the vertical excitation energies at these minima. At the CIS/3-21G level of theory the Stokes shift of 0.12 meV is in somewhat better agreement with the xMCQDPT2//SA2-CASSCF(12,12)/cc-pVDZ value of 0.200 meV than SA2-CASSCF(12,12)/cc-pVDZ, which yields 338 meV, or TD-CAMB3LYP/6-31G(d), which yields 88 meV.
Because the polaritonic states are determined by the coupling between the cavity field and the molecular transition dipole moments, we also compared the transition dipole moment between the different levels of theory ( Figure S3b). As with the electronic energies, there are deviations, but the average transition dipole moments are within 0.1 D.
Because we can compensate the overestimated excitation energy by adding an offset to the cavity dispersion in our simulations, we conclude that for the purposes of our study, the CIS/3-21G method provides an acceptable compromise between desired accuracy and computational efficiency for calculating trajectories of ensembles of Rhodamine molecules strongly coupled to the confined light modes of a Fabry-Pérot micro-cavity.

Minimum Energy Conical Intersection
There are three important points on the S 1 PES of the Rhodamine molecule: (i) the Frank-Condon Results of the optimizations are summarized in Table 1, in which the first colomn contains the absolute S 1 energies in Hartrees, while the second colomn lists the relative energies with respect to the FC energy in kJmol −1 . Because the MECI point is more than 120 kJmol −1 higher in energy than FC point, it is extremely unlikely that the system would sample the MECI region and undergo a radiationless decay into the S 0 ground state on the time scale of our simulations. Instead, the most probably deactivation process for bare Rhodamine is fluorescence, in line with experiment.

Tetracene Model
Tetracene ( Figure S2b) was placed at the center of a rectangular box that was filled with 265 cyclohexane molecules. Intermolecular interactions were modelled with the Gromos96-54a7 unitedatom force field. 35 The simulation box, which contained 2,785 atoms, was equilibrated for 100 ns.
During equilibration, the coordinates of the Tetracene atoms were kept fixed. The LINCS algorithm was used to constrain bond lengths in cyclohexane, 18

Methylene Blue Model
Methylene Blue ( Figure S2c) was placed at the center of a rectangular box that was filled with 2,031 TIP3P water molecules. 16 Non-covalent interactions between Methylene Blue and water were modelled with the Amber03 force field. 15  After the equilibration at the force field (MM) level, the system was further equilibrated for 10 ps at the QM/MM level. The time step was reduced to 1 fs. The Methylene Blue molecule was modelled with density functional theory, using the B97 functional, 38 in combination with the 3-21G basis set. 28 The water solvent was modelled with the TIP3P force field. 16 The QM system experienced the Coulomb field of all MM atoms within a 1.0 nm cut-off sphere and Lennard-Jones interactions between MM and QM atoms were added. The singlet electronic excited state (S 1 ) was modeled with time-dependent density functional theory, 25 using the B97 functional in combination with the 3-21G basis set for the QM region (i.e., TD-B97/3-21G), 28

Molecular Dynamics of Rhodamine-Cavity Systems
After equilibration at the QM/MM level, the Rhodamine molecules in solution were placed with equal inter-molecular distances on the z-axis of a 1D, 4 50 µm long, optical Fabry-Pérot microcavity ( Figure S1). The dispersion of this cavity was modelled with 160 discrete modes (i.e., k z,p = 2πp/L z with 0 ≤ p ≤ 159 and L z = 50 µm). The microcavity was red-detuned with respect to the Rhodamine absorption maximum, which is 4.18 eV at the CIS/3-21G//Amber03 level of theory used in this work, such that the energy of the cavity photon at zero incidence angle (k 0 = 0) was ℏω 0 = 3.81 eV, corresponding to a distance of L x = 0.163 µm between the mirrors ( Figure S1).
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. Both lossless (ℏγ cav = 0 eV or γ cav = 0 ps −1 , Equation 19) and lossy ( ℏγ cav = 0.04 eV or γ cav = 66.7 ps −1 ) mirrors were considered such that the lifetime, τ cav , of the lossy cavity is comparable to the 2-14 fs lifetimes of metallic Fabry-Pérot cavities used in experiments of strong coupling with organic molecules. [39][40][41] With a lifetime of 15 fs and a resonance at 3.81 eV, the Quality-factor, defined as Q = ω cav τ cav , for the lossy cavity would be 87. We note, however, that this Q-factor is artificially For systems without disorder, bright states and dark states are easily distinguishable because dark states lack contributions from cavity mode excitations. In contrast, in the presence of disorder due to molecules adopting different conformations in the simulations (and in reality), the cavity modes are smeared out over many states. Nevertheless, because the cavity modes are not distributed evenly over all states, the upper (UP) and lower (LP) polaritonic branches remain visible, even if the number of eigenstates that make up these branches exceeds the number of cavity modes. Therefore, to have sufficient dark states, loosely defined as states with a negligible cavity mode contribution, the ratio between the number of dark (N dark ) and polaritonic states (N pol ) must be larger than zero. While this is the case when there are 1024 (N dark /N pol = 2.7 ), 512 (N dark /N pol = 1.1), or 256 molecules (N dark /N pol = 0.3) in the cavity, the smallest cavity systems with 160 molecules (N dark /N bright = 0), has no real dark state manifold, which, as we show below, affects the propagation mechanism.
We used a leap-frog Verlet algortihm with a 0.1 fs timestep (0.5 fs, for simulations at diffent Rabi splittings) to perform the integration of the nuclear dynamics, while the unitary propagator in the local diabatic basis was used to propagate the polaritonic degrees of freedom. 9 To avoid selfinterference of the wavepacket, the simulations were terminated far before a wavepacket could reach the periodic boundary of the 1D cavity at z N +1 = z 1 . All results were obtained as averages of at least two simulations. For all cavity simulations we used Gromacs 4.5.3,22 in which the multi-mode Tavis-Cummings QM/MM model was implemented, 5 in combination with Gaussian16. 42

Molecular Dynamics of Tetracene-Cavity Systems
After equilibration at the QM/MM level 1024 Tetracene molecules in cyclohexane were placed with equal inter-molecular distances on the z-axis of a 1D, 4 50 µm long, optical Fabry-Pérot micro-cavity ( Figure S1). The dispersion of this cavity was modelled with 160 discrete modes (i.e., k z,p = 2πp/L z with 0 ≤ p ≤ 159 and L z = 50 µm). The microcavity was red-detuned with respect to the Tetracene absorption maximum, which is 3.87 eV at the CIS/3-21G//Gromos96-54a7 level of theory used in this work, such that the energy of the cavity photon at zero incidence angle (k 0 = 0) was ℏω 0 = 3.55 eV. To maximize the collective light-matter coupling strength, the transition dipole moments of Tetracene molecules were aligned to the vacuum field at the start of the simulation. A leap-frog Verlet algortihm with a 0.1 fs timestep was used to integrate the nuclear dynamics, while the unitary propagator in the local diabatic basis was used to propagate the polaritonic degrees of freedom. 9 To avoid self-interference of the wavepacket, the simulation was stopped before a wavepacket could reach the periodic boundary of the 1D cavity at z N +1 = z 1 . Results were obtained as averages of three independent simulations that were performed with Gromacs 4.5. 3,22 in which the multi-mode Tavis-Cummings QM/MM model was implemented, 5

in combination with
Gaussian16. 42

Molecular Dynamics of Methylene-Blue Cavity System
After equilibration at the QM/MM level 1024 Methylene Blue molecules in solution were placed with equal inter-molecular distances on the z-axis of a 1D, 4 50 µm long, optical Fabry-Pérot micro-cavity ( Figure S1). The dispersion of this cavity was modelled with 160 discrete modes (i.e., k z,p = 2πp/L z with 0 ≤ p ≤ 159 and L z = 50 µm). The microcavity was red-detuned with respect to the Methylene Blue absorption maximum, which is 2.5 eV at the TD-B97/3-21G//Amber03 level of theory used in this work, such that the energy of the cavity photon at zero incidence angle (k 0 = 0) was ℏω 0 = 1.9 eV. To maximize the collective light-matter coupling strength, the transition dipole moments of the Methylene Blue molecules were aligned to the vacuum field at the start of the simulation. With a vacuum field strength of E y = 0.000025 a.u. (0.128 MVcm −1 ) this configuration yields a Rabi splitting of ℏΩ Rabi = 175 meV. Because the purpose of this additional simulation was to confirm the general transport mechanism, rather than explore the effect of molecular or cavity parameters, we only considered an ideal lossless cavity (ℏγ cav = 0 eV or γ cav = 0 ps −1 , Equation

19).
A leap-frog Verlet algortihm with a 0.5 fs timestep was used to integrate the nuclear dynamics, while the unitary propagator in the local diabatic basis was used to propagate the polaritonic degrees of freedom. 9 To avoid self-interference of the wave packet, the simulation was ended before a wavepacket could reach the periodic boundary of the 1D cavity at z N +1 = z 1 . The simulation was where β = 10 −12 m 2 is a coefficient characterising the shape of the wavepacket and k z,m is the in-plane momentum of polariton |ψ m ⟩, given by: with k z,p = 2πp/L z the discrete wavevector in a periodic 1D cavity of length L z .

Off-resonant Excitation
In experiments off-resonant excitation is often realized by pumping a higher-energy electronic excited state of a single molecule. According to Kasha's rule, this molecule then rapidly relaxes into the lowest energy electronic excited state (S 1 ), which is resonant with the cavity. To avoid the computationally highly demanding description of this ultra-fast initial relaxation process, we started our simulations in the S 1 electronic state of one of the molecules. Because we performed our simulations in the adiabatic representation, in which the polaritonic wave functions are delocalized over all molecules, we took a linear combination of adiabatic states that localizes the excitation onto that molecule as initial conditions. To find the expansion coefficients for such localisation, we used the completeness relation to transform the adiabatic eigenstates of the Tavis-Cummings Hamiltonian (|ψ α ⟩, Equation 14) back into the basis of molecular excitations and cavity mode Fock states (i.e., |ϕ j ⟩ =σ + j |ϕ 0 ⟩ for j ≤ N and |ϕ j ⟩ =â † j |ϕ 0 ⟩ for i > N ): where H TC is the matrix representation of the Tavis 3 Simulation Analysis

Contributions of Molecular and Cavity Mode Excitations to the Wavepacket
The total contribution of the molecular excitations (|Ψ exc | 2 ) to the wavepacket (black curve in the fourth column of Figures S8, S9, S5 and S4 and panels d and h of Figures 2 and 4 of the main article), was obtained by summing the projections of the N excitonic basis states (i.e., |ϕ j ⟩ =σ + j |ϕ 0 ⟩ for 1 ≤ j ≤ N , Equation 7) onto |Ψ(t)⟩ (Equation 15) at each time-step of the simulation: (25) where N + n mode is the number of polaritonic states.

Monitoring Wavepacket Dynamics
To monitor the propagation of the wavepacket, 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 (first column of Figures S8, S9, S5 and S4 and panels a and e of Figures 2 and 4 in the main text). We thus represent the density as a discrete distribution at grid points that correspond to the molecular positions, rather than a continuous distribution. In addition to the total probability density, |Ψ(t)| 2 , we also plotted the probability densities of the molecular |Ψ mol (t)| 2 and photonic |Ψ pho (t)| 2 contributions separately (second and third columns of Figures S8, S9, S5 and S4 and panels b, f and c, g, respectively, of Figures 2 and 4 in the main text).
The amplitude of |Ψ exc (t)⟩ at position z j in the 1D cavity (with z j = (j−1)L z /N for 1 ≤ j ≤ N ) is obtained by projecting the excitonic basis state in which molecule j at position z j is excited (|ϕ j ⟩ =σ + j |ϕ 0 ⟩) to the total wave function (Equation 15): where the β m j are the expansion coefficients of the excitonic basis states in polaritonic state |ψ m ⟩ (Equation 14) and c m (t) the time-dependent expansion coefficients of the total wavefunction |Ψ(t)⟩ (Equation 15).
The cavity mode excitations are described as plane waves that are delocalized in real space. We therefore obtain the amplitude of the cavity modes in polaritonic eigenstate |ψ m ⟩ at position z j by Fourier transforming the projection of the cavity mode excitation basis states, in which cavity mode where the α m p are the expansion coefficients of the cavity mode excitations in polaritonic state |ψ m ⟩ (Equation 14) and we normalize 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 is then obtained as the weighted sum over the Fourier transforms: with c m (t) the time-dependent expansion coefficients of the adiabatic polaritonic states |ψ m ⟩ in the total wavefunction |Ψ(t)⟩ (Equation 15).

Average Velocity
To determine the average velocity, v av , of the total wave packet, we first computed the expectation value of the position as a function of time: Then, v av was obtained from the slope of a linear fit to ⟨z(t)⟩.

Photo-Absorption Spectra
Following Lidzey and coworkers, 43 we define the "visibility", I m , of polaritonic state ψ m as the total photonic contribution to that state (i.e., I m ∝ nmax p |α m p | 2 ). Thus, the angle-resolved, or wave vector-dependent, one-photon absorption spectra of the Rhodamine cavity systems were computed from the QM/MM trajectory of the uncoupled Rhodamine as follows: For each frame of this trajectory, the polaritonic states were computed and the energy gaps of these states with respect to the overall ground state (i.e., E 0 , with all molecules in S 0 , no photon in the cavity, Equation 7), were extracted for all cavity modes, p, and summed up into a superposition of Gaussian functions, weighted by |α m p | 2 : Here, I abs (E, k z ) is the absorption intensity as a function of excitation energy E and in-plane

Photo-Luminescence Spectra
As in Berghuis et al., 36 we analysed the photo-luminescence of the propagating wavepacket after non-resonant excitation through a (virtual) mask placed at various distances from the molecule that was excited initially. This analysis thus mimics the Fourier Transform microscopy technique that is used experimentally to determine how different polaritonic states contribute to the propagation. 36 In such imaging experiments, the photoluminescence from a specific area of the sample is collected through a pinhole and converted into momentum space. In the analysis of our simulations we modeled this pinhole, or mask, as a rectangular function of width ∆z = 10 µm centered at z = 10 and 20 µm. At each time-step t, the cavity photoluminescence at the position of the pinhole was collected as follows: |Ψ pho hole (z, t)⟩ = z j ≤z+∆z/2 Second, |Ψ pho hole (z, t)⟩ was Fourier transformed from real space into k z space: such that the total wave function under the mask can be written as |Ψ hole (t)⟩ = |Ψ exc hole (z, t)⟩ ⊗ |Ψ pho hole (k z , t)⟩. Since the polaritonic states {|ψ m ⟩} form a complete basis (Equation 15), the part of the wavepacket under the mask at time t can be expanded as a linear combination of these eigenstates: To obtain the N + n mode expansion coefficients, d m (t), we computed the overlap between the polaritonic eigenstates and the fraction of the wavepacket visible through the pinhhole: Following Lidzey and co-workers, 43 we approximated the photoluminescence detected through the pinhole at time t as the "visibility" of the polaritonic states weighted by their expansion coefficients under the mask, d m (t): 5 with I hole (E, p, z, t) the photoluminescence intensity at time t as a function of excitation energy E, in-plane momentum p (k z,p = 2πp/L z ) and position of the mask z; ∆E m the excitation energy of polaritonic state |ψ m ⟩ at that time frame (∆E m = E m − E 0 , with E 0 the energy of the ground state, |ϕ 0 ⟩, Equation 7); and σ = 0.05 eV the convolution width. The visibility through the pinhole, I hole (E, p, z), was accumulated over 100fs: where s is the number of time steps (1000 used for the pinhole analysis of simulations under off-resonant excitation conditions with no photon losses in this work).

Polariton Transport for Different Numbers of Molecules N
Because the number of molecules we can include in our simulations is necessarily much smaller than in experiments, we investigated how that number affects energy transport by performing simulations with 160, 256, 368, 512, 768 and 1024 molecules. To keep the Rabi splitting (∼ 325 meV) constant, and hence polariton dispersion the same, we scaled the cavity mode volume by the number of molecules N (see Section 2 for details).

On-resonant Excitation
In Figure S4 we show the time-evolution of the probability density of the polaritonic wave function,  (Figures S6 and S7).
When photon loss is included ( Figure S5), population builds up in the ground state (blue curves) Figure S4: Polariton propagation in ensembles of 256, 512 and 1024 molecules (first, second and third row, respectively) strongly coupled to a lossless cavity (i.e., γ cav = 0 ps −1 ), after on-resonant excitation of a wavepacket of LP states centered at the k z -vector corresponding to the maximum LP group velocity. First, second and third columns: total probability density, |Ψ(t)| 2 , probability density of the molecular excitons, |Ψ exc (t)| 2 , and of the cavity mode excitations, |Ψ pho (t)| 2 , respectively, as a function of distance (horizonal axis) and time (vertical axis). The magenta dashed line indicates propagation at the maximum group velocity of the LP (68 µmps −1 ). Fourth column: Contributions of the molecular excitons (black) and cavity mode excitations (red) to |Ψ(t)| 2 as a function of time. Without cavity decay, the population of the ground state (GS, blue) remains zero.
at the expense of states with cavity mode excitations (red). Because the number of dark states increases with N , and the dark states "protect" against decay, the total decay decreases with N , albeit that for the ensemble sizes simulated here, the effect is rather small. Figure S5: Polariton propagation in an ensemble of 256, 512 and 1024 molecules (first, second and third row, respectively) strongly coupled to a lossy cavity (i.e., γ cav = 66.7 ps −1 ), after on-resonant excitation of a wavepacket of LP states centered at the k z -vector corresponding to the maximum LP group velocity. First, second and third columns: total probability density, |Ψ(t)| 2 , probability density of the molecular excitons, |Ψ exc (t)| 2 , and of the cavity mode excitations, |Ψ pho (t)| 2 , respectively, as a function of distance (horizontal axis) and time (vertical axis). The magenta dashed line indicates propagation at the maximum group velocity of the LP (68 µmps −1 ). Fourth column: Contributions of molecular excitons (black) and cavity mode excitations (red) to |Ψ(t)| 2 and population of the ground state (GS, blue) as a function of time. Figure S6: Expectation value of the position (⟨z⟩) (top panels) and mean square displacement (MSD) (bottom panels) of the total time-dependent wavefunction |Ψ(t)⟩ of 256 Rhodamine chromophores strongly coupled to an ideal cavity (i.e., γ cav = 0 ps −1 , panels a and c) and a lossy cavity with γ cav = 66.7 ps −1 (panels b and d), under on-resonant excitation conditions. In the top panels, the black lines represents ⟨z⟩ = ⟨Ψ(t)|ẑ|Ψ(t)⟩/⟨Ψ(t)|Ψ(t)⟩ while the shaded area indicates the root mean squared deviation ⟨(z(t) − ⟨z(t)⟩) 2 ⟩. Figure S7: Expectation value of the position (⟨z⟩) (top panels) and mean square displacement (MSD) (bottom panels) of the total time-dependent wavefunction |Ψ(t)⟩ of 512 Rhodamine chromophores strongly coupled to an ideal cavity (i.e., γ cav = 0 ps −1 , panels a and c) and a lossy cavity with γ cav = 66.7 ps −1 (panels b and d), under on-resonant excitation conditions. In the top panels, the black lines represents ⟨z⟩ = ⟨Ψ(t)|ẑ|Ψ(t)⟩/⟨Ψ(t)|Ψ(t)⟩ while the shaded area indicates the root mean squared deviation ⟨(z(t) − ⟨z(t)⟩) 2 ⟩.

Off-resonant Excitation
In Figure S8 we show the propagation of the wavepacket after off-resonant excitation in a lossless cavity (γ = 0 ps −1 , Equation 19) with 160, 256, 512 and 1024 Rhodamine molecules. In Figure S9 we show the propagation of the wavepacket after off-resonant excitation in a lossy cavity (γ = 66.7 ps −1 for all cavity modes) with 256, 512 and 1024 Rhodamine molecules.
Because the initial wavefunction that localizes the excitation onto a single molecule at the start of the simulation, is a superposition of all states, including the highest-energy, highest-momentum UP states, the fastest propagation in the smaller ensembles (i.e., N < 512) is dominated by wavepackets of these states. Since UP states posses much higher group velocities compared to LP states, these wavepackets cross the periodic boundary around 200 fs and start self-interfering in these ensembles.
We therefore only plot ⟨z⟩ and MSD until such crossings occur in Figures S10, S11 and S12.
In the smallest ensemble, with 160 molecules coupled to the 160 modes of the 1D Fabry-Pérot cavity, there are no dark states and all states are hybrids of molecular and cavity mode excitations.
Therefore, in contrast to the larger ensembles with dark states, all states have group velocity due to their cavity mode excitations, and all of the population propagates ballistically, as evidenced by the quadratic dependence of the mean-square displacement (MSD) on time ( Figure S10). Because in the largest ensembles with many dark states, the propagation appears diffusive, the 160 molecule simulation without dark states provides additional support for the conclusion that transient trapping of the population in stationary dark states is the origin for such diffusion.
Comparing the relative contributions to the wavepacket of the molecular excitons (black) and cavity mode excitations (red) in the population plots for 160, 256, 512 and 1024 Rhodamine molecules in the right column of Figure S8 suggests that as the number of molecules, and hence the number of dark states, increases, the total population that resides in bright light-matter states decreases. As explained in Tichauer et al. 44 , such decrease in bright state population is due to  the latter rapidly decreasing with increasing N .
With a larger fraction of the population residing in the stationary dark states, the wavepacket propagates more slowly, as evidenced by comparing the expectation value of ⟨z⟩ in Figure S10 for 160 molecules (without dark states), Figure S11 for 256 molecules, Figure S12 for 512 molecules and Figure 5 (main text) for 1024 molecules. Average propagation velocities (υ av ) obtained by fitting a linear function to ⟨z⟩, are plotted as a function of the number of molecules in Figure S13. In   to the largest ensemble with 1024 Rhodamine molecules, for which a more linear dependency is observed (Figure 3c, main text). Because a linear dependency of the MSD on time is associated with diffusion, these trends further confirm that stationary dark states are responsible for diffusive propagation.
When cavity losses are added, population transfer into the ground state competes with polariton propagation ( Figure S9). Because only states with a significant contribution from cavity mode excitations can emit, the overall population decays faster for the smaller ensembles, suggesting that dark states can "protect" the excitation from decay under off-resonance excitation conditions and hence increase the lifetime of the excitation. 21

Comparison between Uni-and Bi-directional Cavity Models
In MD simulations of 1024 Rhodamine molecules strongly coupled to a symmetric cavity (i.e.,

Simulations at 0K
We could previously show that population transfer between bright and dark states is mediated by nuclear dynamics. 44 Therefore, in classical MD simulations, a finite temperature is required for such non-adiabatic transitions. To understand the effect of the thermally activated nuclear dynamics on the transport, we repeated the simulations of 1024 Rhodamine molecules strongly coupled to the confined light modes of a cavity (ω 0 = 3.81 eV, 160 modes) with perfect mirrors (γ cav = 0 ps −1 ) at 0 K. This hypothetical temperature was modelled by freezing all nuclear degrees of freedom.

On-resonant Excitation
In Figure S15 we plot the probability density of the wavepacket as a function of time and position after on-resonant excitation at 0K. The wavepacket propagates ballistically at the central group velocity and retains its smooth Gaussian shape, in contrast to simulations at 300 K, in which peaks appear at the locations of molecules due to the fluctuations in their geometries and hence energies i.e. diagonal energy disorder. 8 Importantly, at 0K, no transition from ballistic to diffusive motion is observed ( Figure S16d), in line with experiments on inorganic materials at low temperatures. 46 Figure S15: Polariton propagation after on-resonant excitation of a wavepacket of LP states centered at the k z -vector corresponding to the maximum LP group velocity in a lossless cavity (i.e., γ cav = 0 ps −1 ) containing 1024 Rhodamine molecules at 0 K. Panels a, b and c: total probability density |Ψ(t)| 2 , probability density of the molecular excitons |Ψ exc (t)| 2 , and of the cavity mode excitations |Ψ pho (t)| 2 , respectively, as a function of distance (horizontal axis) and time (vertical axis). The purple dashed line indicates propagation at the maximum group velocity of the LP (68 µm ps −1 ). Figure S16: Expectation value of the position (⟨z⟩, top panels) and mean square displacement (MSD, bottom panels) of the total |Ψ(t)⟩ (left panels), excitonic |Ψ exc (t)⟩ (middle panels), and photonic |Ψ pho (t)⟩ (right panels) wave functions of 1024 Rhodamine molecules strongly coupled to a lossless cavity (γ cav = 0 ps −1 ) after on-resonant excitation at 0 K. As a reference, the curves corresponding to the same system at 300 K are plotted as grey dashed lines.

Off-resonant Excitation
In Figure S17 we plot the probability density of the wavepacket as a function of time and position after off-resonant excitation at 0K. Compared to the time-evolution of the wavepacket at 300 K, the propagation is reduced at 0K, and most of the excitation remains localized at the molecule that was initially excited ( Figure S17b). The reduced mobility of the total wavepacket is visible also in the evolution of the expectation value of ⟨z⟩ in Figure S18a. In contrast to the total wavepacket and the contribution of the molecular excitons, the part of the wavepacket that consists of the cavity mode excitations propagates ballistically ( Figure S18c and f). Thus, while most of the wavepacket remains localized on the molecule that was excited initially, a small fraction of bright states propagate ballistically. Because there are no nuclear dynamics at 0K, this propagation is solely driven by constructive and destructive interferences between bright states that were occupied in the initial superposition, and that are due to differences in their phases (i.e., Rabi oscillations: Figure S17: Polariton propagation after off-resonant excitation of a single molecule located at z = 5 µm into S 1 at 0 K. Panels a, b and c: total probability density |Ψ(t)| 2 , probability density of the molecular excitons |Ψ exc (t)| 2 , and of the cavity mode excitations |Ψ pho (t)| 2 , respectively, as a function of distance (horizontal axis) and time (vertical axis), in a cavity with perfect mirrors (i.e., γ cav = 0 ps −1 ). The purple and yellow dashed lines indicate propagation at the maximum group velocity of the LP (68 µm ps −1 ) and UP (212 µmps −1 ), respectively.
Figure S18: Expectation value of the position (⟨z⟩, top panels) and mean square displacement (MSD, bottom panels) of the total, |Ψ(t)⟩ (left panels), excitonic, |Ψ exc (t)⟩ (middle panels), and photonic, |Ψ pho (t)⟩ (right panels), wave functions of 1024 Rhodamine molecules strongly coupled to an lossless cavity (γ cav = 0 ps −1 ) after off-resonant excitation at 0 K. As a reference, the curves corresponding to the same system at 300 K are plotted as grey dashed lines. Note that we plot both ⟨z⟩ and MSD only until t = 200 fs, since at later times the part of the wave packet composed of upper polariton states has crossed the periodic boundary at z = 50 µm.

Non-Adiabatic Transitions
As shown in the previous section, polariton propagation at 300 K differs significantly from motion at 0K. To understand if these differences solely arise from dynamically induced disorder or also involve non-adiabatic transitions, we plot the time-evolution of the expansion coefficients of the adiabatic LP, UP and dark states (i.e., m |c m (t)| 2 with m ∈ LP (160 low-energy states), UP (160 higher-energy states), and DS (the remaining states), respectively) in Figure S19 for 256 molecules in an ideal lossless cavity after on-resonant and off-resonant excitation. The variation in the magnitude of these coefficients with time suggest that there is significant non-adiabatic population transfer between bright and dark states in our simulations.

Effect of Positional Disorder on Polariton Propagation
In our simulations the solvated Rhodamine molecules were located equidistantly along the z-axis of the 1D cavity. To verify if such ordering affects the propagation process, we repeated a simulation under off-resonant excitation conditions with 512 Rhodamines randomly distributed along the z-axis of the cavity. This was achieved by adding a random displacement to the initial equidistant positions. We only considered a perfect cavity with no losses (γ cav = 0 ps −1 ) and with a vacuum field strength of 0.0000707107 au (0.364 MVcm −1 ), yielding a Rabi splitting of ∼ 325 meV. Direct visual comparison between the wavepacket dynamics in this system, shown in Figure S20, to that in the system with an equal spacing between the molecules (third row in Figure S8), suggest that positional disorder does not influence polariton propagation, nor the time evolution of the wave function composition ( Figure S20d and Figure S8h). Figure S20: Polariton propagation after an off-resonant excitation into the S 1 electronic state of a single molecule located at z = 5 µm, in an ensemble of 512 Rhodamine molecules randomly distributed along the z-axis of an ideal cavity (i.e., γ cav = 0 ps −1 ). Panels a, b and c: total probability density, |Ψ(t)| 2 , probability density of the molecular excitons, |Ψ exc (t)| 2 , and of the cavity mode excitations, |Ψ pho (t)| 2 , respectively, as a function of distance (horizontal axis) and time (vertical axis). Panel d: Contributions of molecular excitons (black) and cavity mode excitations (red) to |Ψ(t)| 2 and population of the ground state (GS, blue) as a function of time.

Effect of Initial Energy Disorder on Ballistic Transport
Starting simulations with the same molecular geometries, and thus identical excitation energies and transition dipole moments, facilitates the identification of bright and dark states. However, in experiments the molecular geometries span a distribution weighted by their Boltzmann factors. To assess the effect of such initial (thermal) disorder on polariton propagation after resonant excitation, we repeated a simulation of 256 Rhodamines inside a lossless cavity (γ cav = 0 ps −1 ) with a vacuum field strength of E y = 0.0001 a.u. (0.514 MVcm −1 ), with different initial geometries and velocities.
These configurations were selected from the equilibrium QM/MM trajectory described in section 2. Figure S21: Panel a: Total probability density, |Ψ| 2 , of 256 Rhodamine molecules strongly coupled to an ideal cavity with (green line at t = 0 fs and blue line at t = 100 fs) and without (black line at t = 0 fs and red line at t = 100 fs) initial energy disorder after on-resonant excitation. Panel b: Mean square displacement (MSD) of the total wave function, |Ψ(t)⟩, with (black) and without (gray) initial geometric disorder among the molecules. The green and red dashed lines are quadratic fits to the MSD in the ballistic phase.
As shown in panel a of Figure S21, with energetic disorder, the initial wavepacket is no longer a smooth Gaussian curve (green line), but its propagation is very similar to that of a wavepacket starting without such disorder (black line): Initially, there is ballistic motion at a velocity close to the maximum group velocity of the LP, followed by diffusion (see Figure S21b). However, because with the initial energy disorder non-adiabatic population transfers into lower energy dark states occur more readily from the start of the simulation, the transition between ballistic and diffusive transport happens earlier.

Effect of Rabi Splitting on Polariton Propagation
The results of our MD simulations suggest that reversible population transfers between bright polaritonic states with group velocity and dark states without group velocity, play an important role in the propagation of a polaritonic wavepacket. Because the non-adiabatic coupling vectors that control these population transfers are inversely proportional to the energy gap between the states, 44 altering the overlap between dark and bright states by changing the Rabi splitting could impact polariton propagation. To facilitate direct comparison, the same random velocities were used for each Rabi splitting. The latter was realized by using the index of the molecule as a random seed for selecting the velocities in the first series of simulations, and adding 1024 to that index in the second series.

On-resonant Excitation
In Figure  Based on previous results, which suggest that increasing the vacuum field strength decreases the overlap between the bright polaritonic branches and dark state manifold, 21 we speculate that the decrease in group velocity upon increasing the Rabi splitting, may be compensated by a reduction of the population transfer into the dark states. To minimize that population transfer and extend the for typical organic materials would require entering the ultra-strong coupling regime, where the Rabi splitting is more than 10% of the excitation energy. 1 However, in this regime, the Tavis-Cummings models is no longer strictly valid as the RWA breaks down, and the counter-rotating terms (σ + jâ † kz andσ jâkz ), as well as the dipole self energy (Ĥ DSE = 1 2ϵ 0 Vcav u y · N iμ i 2 ) would need to be included, 47 which is computationally intractable for the system sizes used in this work.
Nevertheless, even if our results do not reveal a clear trend, there may be a trade-off between the absorption line width of the material and the Rabi splitting, for optimal exciton transport. However, because the focus of the current work is on the mechanism of polariton propagation, we consider finding that trade-off beyond the scope of our work, and will address the effect of the Rabi splitting on polariton transport in future work.

Off-resonant Excitation
In Figure S22b we show the Mean Square Displacement of the total wavepacket for the various Rabi splittings after off-resonant excitation of a single molecule located at z = 5 µm. As for on-resonance excitation, we do not observe a clear correlation between the Rabi splitting on the one hand and the Mean Square Displacement on the other hand. Very many additional simulations are probably needed to extract a trend. However, because the focus of the current work is on the mechanism of polariton propagation, we consider a thorough investigation of the effect of the Rabi splitting beyond the scope of the current work. Yet, since we anticipate that the Rabi splitting may have an effect that could possibly be exploited and can be directly controlled via the molecular concentration, we will address the effect of the concentration on polariton transport in future work.

Contribution of High-Energy Polariton States to the Wavepacket Propagation
Following off-resonant excitation into the S 1 electronic state of a single Rhodamine molecule, we observe that the photonic part of the wavepacket has a small front that propagates faster than the rest of the wavepacket. In Figure S23 we show the probability density of the cavity mode excitations in the 512 molecule cavity system at 50 and 100 fs. The wavepacket is composed of a slower and faster component, with the front of the faster component reaching distances corresponding to the maximum group velocity υ max UP = 212 µmps −1 of the UP branch. To confirm that the fastest components in the wavepacket are due to propagation of population in the UP states that get rapidly populated due to Rabi oscillations (Equation 24), we analyzed the photoluminesce through 10 µm wide pinholes located at z = 10 µm and at z = 20 µm (see section 3.5 for details of these calculations). 36 The photoluminescence spectra collected through these pinholes are plotted as a function of the in-plane momentum (k z ) in Figures S23b and c. While at z = 10 µm the photoluminescence is due to emission from both LP and UP states, emission through the pinhole at z = 20 µm originates exclusively from UP states. These observations thus confirm that the fastest moving part of the wave packet is composed of UP states. Figure S23: Panel a: Probability densities of Ψ pho 2 at t = 50 and 100 fs after off-resonant excitation of a cavity including 512 Rhodamine molecules. The initial excitation is localised on molecule placed at z = 5 µm. The dashed-dotted lines indicate the product between the maximum group velocity υ max UP = 212 µmps −1 and the simulation time, and thus corresponds to where something moving at that group velocity would be located. Panels b and c: Photoluminescence spectra collected through pinholes located at z = 10 µm and z = 20 µm, respectively.

Polariton Propagation in Tetracene and Methylene Blue cavities
With an absorption linewidth of 100 meV (FWHM), a Stokes shift of 120 meV, and a minimum energy conical intersection sufficiently above the Franck Condon region to prevent internal conversion (1.3 eV), our Rhodamine model captures the key photophysical features of typical organic dye molecules. Therefore, we consider the propagation mechanism observed in our simulations, representative for organic polaritons in general. Nevertheless, to demonstrate that the mechanism is robust with respect to changes in the optical properties of the molecules, we repeated the sim-

Tetracene
In Figure S24 we show polariton propagation in a cavity filled with Tetracene in cyclohexane following an on-resonant excitation of a wavepacket of LP states for both a lossless (γ cav = 0 ps −1 ) and a lossy cavity (γ cav = 66.7 ps −1 ). After a very short ballistic phase, which is only discernible as the steep initial rise of ⟨z⟩ in the top panels of Figure S25, propagation becomes diffusive, as evidenced by a linear time-dependence of the Mean Square Displacement ( Figure S25c and d for the lossless and lossy cavity, respectively). The absence of a clear coherent transport regime is due to the large absorption linewidth of Tetracene. As we could show previously, 21 the rate at which population transfers from the LP into the dark states depends on the energetic overlap between the LP and the molecular density of states, with the latter matching the absorption spectrum). Because for Tetracene the absorption spectrum is much broader than the Rabi splitting, this overlap is very large and hence the population transfer is very fast. As a consequence, the excitonic contributions start dominating the wavepacket immediately after the excitation, ending the ballistic propagation phase earlier than in cavities filled with Rhodamine, for which the overlap between the LP and the molecular dark states is much smaller. Because these transfers are reversible, the propagation continues in a diffusive manner. Adding losses does not change this picture, but because the coherent phase is much shorter than in cavities filled with Rhodamine, also the total radiative decay is smaller (bottom rows in Figures 2 and S24). Figure S24: Polariton propagation in cavities filled with Tetracene molecules after on-resonant excitation of a wavepacket of LP states, centered at z = 5 µm. Panels a, b and c: total probability density |Ψ(t)| 2 , probability density of the molecular excitons |Ψ exc (t)| 2 and of the cavity mode excitations |Ψ pho (t)| 2 , respectively, as a function of distance (horizontal axis) and time (vertical axis), in a cavity with perfect mirrors (i.e., γ cav = 0 ps −1 ). The magenta dashed line indicates propagation at the maximum group velocity of the LP (68 µmps −1 ). Panel d: Contributions of molecular excitons (black) and cavity mode excitations (red) to |Ψ(t)| 2 as a function of time in the perfect cavity. With no photon losses, no ground state population (blue) can build up. Panels e, f, and g: |Ψ(t)| 2 , |Ψ exc (t)| 2 and |Ψ pho (t)| 2 , respectively, as a function of distance (horizontal axis) and time (vertical axis), in a lossy cavity (γ cav = 66.7 ps −1 ). Panel h: Contributions of the molecular excitons (black) and cavity mode excitations (red) to |Ψ(t)| 2 as a function of time. The population in the ground state, created by radiative decay through the imperfect mirrors, is plotted in blue.
5 µm. As in cavities filled with Rhodamine ( Figure 3 in the main text), propagation is diffusive, driven by reversible exchange of population between bright states with group velocity and dark states without group velocity. Because for the on-resonant excitation, the initial coherent ballistic phases is very short due to the broad absorption of Tetracene strongly overlapping with the LP branch, the Mean Square Displacements are similar for both excitation schemes, in particular for the lossless cavity ( Figures S25 and Figures S27). This is in contrast to cavities filled with Rhodamine for which the longer lasting coherent phase with ballistic propagation, leads to a faster initial rise of the Mean Square Displacement, until the turn-over into the diffusion regime occurs ( Figure 2 in the main text). Figure S26: Polariton propagation in cavities filled with Tetracene molecules after off-resonant excitation of a molecule located at z = 5 µm. Panels a, b and c: total probability density |Ψ(t)| 2 , probability density of the molecular excitons |Ψ exc (t)| 2 and of the cavity mode excitations |Ψ pho (t)| 2 , respectively, as a function of distance (horizontal axis) and time (vertical axis), in a cavity with perfect mirrors (i.e., γ cav = 0 ps −1 ). The magenta and yellow dashed lines indicate propagation at the maximum group velocity of the LP (68 µmps −1 ) and UP (220 µmps −1 ), respectively. Panel d: Contributions of molecular excitons (black) and cavity mode excitations (red) to |Ψ(t)| 2 as a function of time in the perfect cavity. With no photon losses, no ground state population (blue) can build up. Panels e, f, and g: |Ψ(t)| 2 , |Ψ exc (t)| 2 and |Ψ pho (t)| 2 , respectively, as a function of distance (horizontal axis) and time (vertical axis), in a lossy cavity (γ cav = 66.7 ps −1 ). Panel h: Contributions of the molecular excitons (black) and cavity mode excitations (red) to |Ψ(t)| 2 as a function of time. The population in the ground state, created by radiative decay through the imperfect mirrors, is plotted in blue. Figure S27: Top panels: Expectation value of the position of the total-time dependent wavefunction ⟨Ψ(t)|ẑ|Ψ(t)⟩/⟨Ψ(t)|Ψ(t)⟩ after off-resonant excitation of lossless (a) and lossy cavities (b) filled with Tetracene molecules. The black lines represent ⟨z⟩ while the shaded area around the lines represents the root mean squared deviation (RMSD, i.e., ⟨(z(t) − ⟨z(t)⟩) 2 ⟩). Bottom panels: Mean Square Displacement (MSD, i.e., ⟨(z(t) − ⟨z(0)⟩) 2 ⟩) in the ideal lossless cavity (c) and the lossy cavity (d).

Methylene Blue
Because the broad absorption line-width of Tetracene reduces the duration of the ballistic propagation regime in Tetracene to a few fs, we also performed a simulation of Methylene Blue, which has an absorption line-width that is narrower than that of both Tetracene and Rhodamine. In Figure S28, we show the wavepacket dynamics after exciting a wavepacket of LP states. Indeed, with a reduced overlap between the LP states and the molecular dark states, ballistic propagation is observed. As in Rhodamine (Figure 2 in main text), the ballistic regime ends when the excitonic contributions to the total wave function become dominant, which happens around 75 fs for the cavity filled with Methylene Blue ( Figure S28d). In line with results for Rhodamine, the Mean Square Displacement then changes from a quadratic to linear time dependence, as indicated by the quadratic (magenta) and linear (cyan) fits in Figure S29b. Thus, after on-resonant excitation of the cavity filled with Methylene Blue, propagation proceeds as in cavities containing Rhodamine or Tetracene, with a short ballistic phase when bright states are predominantly populated, followed by a diffusion phase when the population mostly resides inside the dark state manifold.
The additional simulations thus support the generality of a polariton transport mechanism based on reversible exchange of population between propagating polaritonic states and stationary dark states. Our results furthermore suggest that QM/MM simulations open up the possibility to systematically investigate the role of various molecular and cavity parameters on exciton transport.
While further simulations are needed to disentangle these effects and formulate design rules for enhancing exciton transfer with optical resonators, the results we have obtained so far already suggest that this process can be controlled by tuning the overlap between the absorption line-width of the molecules and Rabi splitting of the strongly coupled light-matter systems. Figure S28: Polariton propagation in cavities filled with Methylene Blue molecules after onresonant excitation of a wavepacket of LP states centered at z = 5 µm. Panels a, b and c: total probability density |Ψ(t)| 2 , probability density of the molecular excitons |Ψ exc (t)| 2 and of the cavity mode excitations |Ψ pho (t)| 2 , respectively, as a function of distance (horizontal axis) and time (vertical axis), in a cavity with perfect mirrors (i.e., γ cav = 0 ps −1 ). The magenta dashed line indicates propagation at the maximum group velocity of the LP (145 µmps −1 ). Panel d: Contributions of molecular excitons (black) and cavity mode excitations (red) to |Ψ(t)| 2 as a function of time in the perfect cavity. With no photon losses, no ground state population (blue) can build up.