Identifying Vibrations that Control Non-adiabatic Relaxation of Polaritons in Strongly Coupled Molecule-Cavity Systems.

The strong light-matter coupling regime, in which excitations of materials hybridize with excitations of confined light modes into polaritons, holds great promise in various areas of science and technology. A key aspect for all applications of polaritonic chemistry is the relaxation into the lower polaritonic states. Polariton relaxation is speculated to involve two separate processes: vibrationally assisted scattering (VAS) and radiative pumping (RP), but the driving forces underlying these two mechanisms are not fully understood. To provide mechanistic insights, we performed multiscale molecular dynamics simulations of tetracene molecules strongly coupled to the confined light modes of an optical cavity. The results suggest that both mechanisms are driven by the same molecular vibrations that induce relaxation through nonadiabatic coupling between dark states and polaritonic states. Identifying these vibrational modes provides a rationale for enhanced relaxation into the lower polariton when the cavity detuning is resonant with specific vibrational transitions.

Recent observations of low-threshold lasing,[1−4] of enhanced exciton transport,[5−9] and of modified reactivity under strong light–matter coupling[10−14] suggest that the strong coupling regime may have great potential in material science.[15,16] This regime is reached when the rate of energy exchange between excitations in the material and confined light modes, as exist for example in Fabry–Pérot cavities (Figure a) or near plasmonic surfaces, exceeds the intrinsic decay rates of both the material and photonic modes.[17] Under these conditions, excitations in the material hybridize with the confined light modes to form new light–matter states, called polaritons.[18]

Figure 1

Panel a shows a schematic representation of a Fabry–Pérot cavity filled with tetracene (Tc) molecules (not to scale). Panel b shows the (normalized) absorption as a function of wave vector (k = 2πn/Lcav) of a one-dimensional cavity of length Lcav = 15 μm with a vacuum field strength of 0.000176777 au (0.909024 MVcm–1) and a cavity resonance at k = 0 μm–1 of 2.69 eV, filled with 512 Tc molecules. The dispersion of the empty cavity is shown as a dashed white line, as is the excitation energy maximum of Tc at 300 K (3.22 eV, evaluated at the TDA-CAMB3LYP/6-31G(d) level of theory). The absorption (violet) and emission spectra (turquoise) of bare Tc at 300 K are plotted on the right y axis. The upper (UP) and lower polariton (LP) branches are indicated.

Within the rotating wave approximation (RWA), valid under weak driving conditions typically employed in strong coupling experiments, the multimode Tavis–Cummings model provides conceptual insight into polariton formation when N optically active molecules strongly couple to the nmode confined light modes of an optical cavity:[19,20]Here, σ̂+ (σ̂–) is the operator that excites (de-excites) the jth molecule at position z with excitation energy hν (horizontal dashed white line in Figure b) from the electronic ground (excited) state |S0⟩ (|S1⟩) to the electronic excited (ground) state |S1⟩ (|S0⟩); â (â†) is the photonic annihilation (creation) operator for an excitation of a cavity mode with wave vector k and energy ℏωcav(k) (curved dashed line in Figure b). The strong light–matter interaction is modeled within the dipolar or long-wavelength approximation through g:[21]where ucav is a unit vector in the direction of the electric field of the cavity vacuum field, ϵ0 is the vacuum permittivity, Vcav is the cavity mode volume, and μTDM is the transition dipole moment of molecule j.

The eigenstates of the multimode Tavis–Cummings Hamiltonian are coherent superpositions of excitations in the molecules and of excitations in the cavity modes:with eigenenergies E. The β and α expansion coefficients reflect the contribution of the molecular excitations (|S1⟩) and of the cavity modes (|1⟩) to polariton |ψ⟩. Due to thermal fluctuations, the excitation energies and transition dipole moments of the molecules span a distribution rather than a single value. Therefore, at thermal equilibrium, the expansion coefficients are all different.

For large, realistic molecule–cavity systems, only few states have a significant contribution from the cavity modes.[22] These bright states can absorb light and because they inherit the dispersion of the cavity modes, they form the wave-vector dependent lower (LP) and upper polariton (UP) branches, as shown in Figure b. These branches are separated by the Rabi splitting, defined as the energy gap between the UP and LP at the wave vector where the molecular excitation energy matches the cavity mode energy. In contrast, the large majority of states lack such photonic contribution and are therefore dark. Without contributions from the cavity modes, dark states have energies that do not depend on the wave vector and thus lack dispersion.

While entropy favors population of the much more numerous dark states,[23] altered chemistry, enhanced transport, and lasing are attributed to the LP. Therefore, it is essential to understand how to maximize transfer into LP states. Optical pumping provides a convenient way to populate those states directly via resonant excitation into the LP branch[24] or indirectly either via resonant excitation into the UP branch[25,26] or via nonresonant excitation into an uncoupled electronic excited state of a molecule.[27−31] Alternatively, molecules can also be promoted into an (uncoupled) excited state via electrical injection.[32] In contrast to direct excitation into the LP, the indirect routes require relaxation from the pumped state into the LP.[33] Although polariton relaxation has been investigated to a great extent in both experiments[26−31,34−39] and theory,[22,40−54] it remains unclear what drives this process or how to reverse-engineer the insights from these works for rational design of new molecule–cavity systems. The latter requires relaxation to be predicted from first-principles with quantum chemistry methods, which is the purpose of the present work.

On the basis of previous theoretical models,[40,41,43,47,50,55] as well as results from both photoluminescence[29,31,56−58] and polariton condensation experiments,[1−3] two main mechanisms have been proposed for driving relaxation from the dark state manifold into the LP branch: vibrationally assisted scattering (VAS) and radiative pumping (RP). In VAS, transitions between dark and bright states are accompanied by discrete changes in the eigenstates of specific molecular vibrations.[29,31,40,41] Therefore, VAS is manifested by an increased photoluminescence intensity from points on the LP branch with energy gaps to the dark state manifold that match the frequency of these vibrational modes. Although energy gaps at which such enhancements are observed often coincide with molecular Raman frequencies,[26,29,56] not all Raman-active modes contribute to VAS.[31] Consequently, it remains difficult to decide in advance which vibrational modes to target through the cavity design in order to control relaxation into the LP. The second mechanism, radiative pumping (RP), is speculated to involve direct exchange of photons between uncoupled molecules and polaritonic states, but the details of this mechanism are even less understood. On the basis of a reasonable agreement between the molecular fluorescence and polaritonic photoluminescence lifetimes,[57,58] it was proposed that a photon, emitted by an uncoupled molecule is reabsorbed by a polaritonic state within the cavity and from there emitted to the outside.

To acquire atomistic insights into these mechanisms and understand their driving forces, we performed multiscale molecular dynamics (MD) simulations,[54,59] in which we tracked the relaxation from the dark state manifold into the LP. As model system, we used tetracene (Tc), a member of the acene family of polycyclic aromatic compounds that have been extensively used in strong coupling experiments.[1,60,61] The electronic ground (S0) and excited (S1) states of Tc molecules were modeled with density functional theory (DFT),[62] and time-dependent density functional theory (TDDFT),[63] within the Tamm–Dancoff approximation (TDA),[64] respectively, using the CAM-B3LYP functional[65,66] in combination with the 6-31G(d) basis set.[67] At this level of theory, the vertical excitation energy of Tc is 3.22 eV in the S0 minimum geometry, while the energy gap to the ground state is 2.85 eV in the S1 minimum. Using Ehrenfest,[68] or mean-field, dynamics, we computed classical MD trajectories of N Tc molecules with their geometric centers evenly distributed along the z axis of the cavity (Figure S1) and with their electronic transition dipole moments aligned to the vacuum field of both single-mode and multimode optical cavities. These trajectories thus evolved on the mean-field potential energy surface provided by the total time-dependent polaritonic wave function, |Ψ(t)⟩, which was expanded in the basis of the time-independent eigenfunctions of the Tavis–Cummings Hamiltonian (eq ) that depend only on the positions of the atoms in the molecules:[52,54]Here, the c(t) are the time-dependent expansion coefficients that are integrated along with the classical trajectories, using the unitary propagator in the local diabatic basis, which is inherently stable also in the case of (trivial) crossings between potential energy surfaces.[69] A more detailed description of the methods employed in this work is provided as Supporting Information (SI).

To understand the nature of the dark states, defined here as states with a total contribution of cavity mode excitations below a numerical threshold (i.e., ∑ α < 0.05, eq ), and, in particular, identify the lowest-energy dark state (DS0), we simulated the dynamics after nonresonant excitation into the S1 electronic state of one of the N molecules (molecule 1). These simulations were performed with N = 1, 2, 4, 16, 32, 64, and 128 Tc molecules strongly coupled to a single confined light mode with energy ℏωcav = 3.22 eV. Here, we use a single rather than a multimode cavity to control the energy gap between the dark states and the single LP state via the Rabi splitting, defined asThe Rabi splitting was kept constant at ∼429 meV for all systems by scaling the mode volume of the cavities, Vcav, with the number of molecules, N, at the start of the simulation (Table S1).

The results of these simulations, discussed in the SI, suggest that after nonresonant excitation into the S1 electronic state of the first molecule, this molecule relaxes into the minimum on its S1 potential energy surface, while the other molecules remain in their S0 minimum geometries (Figure S3, SI). This finding is in line with results from previous simulations, which also suggest that in the strong coupling regime molecules can still access the S1 minimum.[70−73] We next inspected the eigenstates of the Tavis–Cummings Hamiltonian in this combination of molecular configurations. The contribution of the molecule in the S1 minimum energy geometry (molecule 1) to this lowest energy dark state (i.e., |β1DS|2, eq ), listed in Table S1, rapidly converges to unity upon increasing the number of molecules. Since in experiment the number of molecules is much larger than we can include in our simulations, these results suggest that in reality one-photon excitations are almost completely localized onto single molecules in the dark state manifold, confirming an earlier assumption by Agranovich and co-workers about the nature of the lowest energy dark state.[40,41]

Before addressing the molecular mechanism of relaxation from the lowest energy dark state (DS0) into the LP, we first investigated how the relaxation rate depends on the number of molecules, N, that are strongly coupled to the cavity. This information is important, because in our simulations we typically employ a much smaller number of molecules than in experiments. We thus performed simulations with N = 16, 32, 64, and 128 Tc molecules strongly coupled to a single confined light mode with the same energy as before (ℏωcav = 3.22 eV). The simulations were started in the lowest energy dark state DS0 (i.e., |cDS(0)|2 = 1 in eq ), identified above, in which the excitation is localized onto a single Tc molecule that has relaxed into the S1 minimum, far below the cavity frequency and hence no longer that strongly coupled.[40,41] Initial velocities were randomly selected from a Maxwell–Boltzmann distribution at 300 K. For quantitative comparisons between the different ensemble sizes, we assigned the same initial velocities to the separate molecules in each of the ensembles by setting the seed of the random number generator that determines what velocities are selected for each atom in the molecule, equal to the index of that molecule. These initial conditions can cause inhomogeneous broadening, but only on longer time scales.

In Figure a, we plot the time evolution of the population in the LP state for the four ensemble sizes. The total population transferred from the lowest energy dark state DS0 into the LP during the first 10 fs of the simulation decreases with increasing numbers of molecules. To understand this trend, we invoke the results from first-order time-dependent perturbation theory applied to a two-level system. Because at the start of the simulation the DS0 and LP are separated from all other states by at least ∼395 meV and ∼210 meV, respectively, we consider this a sufficiently valid approximation for the early time evolution of this system. Within this approximation, the population of the LP is given by the Rabi formula:where cLP is the expansion coefficient of the LP state in the total polaritonic wave function (eq ); ℏωDS = ELP – EDS, the energy gap between DS0 and the LP; and VDS, the perturbation term, which is the nonadiabatic coupling between these states. The nonadiabatic coupling, or vibronic coupling, determines the mixing of polaritonic states as a result of molecular vibrations.[74−76] The nonadiabatic coupling VDS depends on the overlap between the nonadiabatic coupling vector dDS and the velocities of the atoms Ṙ:[77]where Ṙ is a vector containing the velocities of all atoms in the system, and the nonadiabatic coupling vector dDS is given byThe direction of this vector determines in what direction population is transferred between the adiabatic states, while its magnitude determines how much population is transferred.

Figure 2

(a) Population of the LP after excitation into the lowest energy dark state (DS0) for N = 16 (black), 32 (red), 64 (green) and 128 (blue) Tc molecules strongly coupled to a single confined light mode. The vacuum field of the cavity is adjusted to achieve a Rabi splitting of 429 meV. (b) Zoom-in on the population during the first femtosecond. Open circles are the data points from the trajectories, continuous lines fits to the Rabi formula (eq ). (c) Absolute value of the perturbation term |VDS|, obtained from the fits in panel b, for the various ensemble sizes of strongly coupled molecules. Open circles are data points; the continuous line is a fit to a/√N + b, with a = 0.10 fs–1 and b = 4 × 10–3 fs–1.

Because we keep the Rabi splitting constant by scaling the cavity mode volume with the number of molecules, the nonadiabatic coupling vector is inversely proportional to the square root of the number of molecules strongly coupled to the cavity (see SI for a derivation). Indeed, fitting the Rabi formula (eq ) to the population of the LP during the first femtosecond (Figure c) yields estimates for |VDS| that scale approximately as 1/√N. Importantly, because the number of molecules that we can include in our simulations is much smaller than the number of molecules in a typical Fabry–Pérot cavity, we overestimate the rate at which the LP is populated from the dark state manifold as compared to experiment. Conversely, the transfer from the LP into the dark state manifold should be comparable to experiment, as the number of dark states into which the population can transfer (i.e., the density of final states) increases with N, canceling the 1/N dependence of the |VLP,DS|2 perturbation term in eq . In principle, therefore, under equilibrium conditions, we would overestimate the population in the LP, which should be taken into account when comparing results between simulations and experiments. Nevertheless, because the rate decreases with increasing N,[22,51,78] we speculate that in real Fabry–Pérot cavities with approximately 106 molecules inside the mode volume,[51,78] the lifetime of the DS0 state will be sufficiently long to reach thermal equilibrium, in particular if the energy gap to the LP states is large. We therefore assume that before relaxation from the DS0 state into the LP occurs, all molecules are in their vibrational ground state.

In the SI, we also show that the nonadiabatic coupling vector dDS connecting the lowest energy dark state (DS0) and the LP is dominated by gradients of the Tavis–Cummings Hamiltonian (eq ) with respect to atomic displacements of the molecule that is in the S1 minimum geometry. Indeed, in simulations in which the atoms of that molecule are kept fixed, there is no population transfer from DS0 into the LP, whereas in simulations with the other molecules fixed, transfer is observed (Figure S5, SI).

Because the molecular displacements responsible for nonadiabatic transitions between the polaritonic states are linear combinations of the vibrational modes, we computed these modes within the harmonic approximation by diagonalizing the Hessian and plotted their overlaps with the nonadiabatic coupling vector as a function of the vibrational mode energy (i.e., ℏωvib for 1 ≤ n ≤ 84) in Figure a. To understand which of the vibrational modes can induce population transfer between the DS0 and LP, we performed 84 short MD simulations of the 16 molecule cavity system. These simulations were initiated in the lowest energy dark state, DS0, in which the first Tc molecule is in the S1 minimum geometry, while the other 15 molecules are in the S0 geometry. In each simulation, we selectively activated one of the 84 vibrational modes of the molecule in the S1 minimum energy geometry by providing initial atomic velocities in the direction of the displacements of that mode (see SI for details). To quantify the extent of population transfer from the DS0 state into the more photonic LP state during the simulation, we projected the excitation of the single cavity mode (i.e., ⟨1|⟨S01S02...S0S0|, eq ) onto the total time-dependent polaritonic wave function (|Ψ(t)⟩, eq ). The top panel in Figure b shows this photonic weight (i.e., |∑∑cα|2) at 1 fs as a function of the vibrational energy of the mode along which the initial velocities were directed. The observation that transfer predominantly occurs if we activate vibrational modes that overlap with the nonadiabatic coupling vector suggests that the relaxation from the dark state manifold into the LP is selectively mediated by these vibrations.

Figure 3

(a) Overlap between the nonadiabatic coupling vector dDS (eq ) and the vibrational modes of Tc, plotted as a function of their energies. (b) Total photonic weight (∑∑cα|2) at 1 fs in simulations of 16 and 512 Tc molecules coupled to a single-mode cavity and of 512 Tc molecules coupled to a Fabry–Pérot cavity modeled with 60 modes (Figure ). The 84 trajectories of these systems were started in the DS0 state, with initial atomic velocities along one of the vibrational modes of the Tc molecule in the S1 geometry. (c) Calculated infrared (IR) and Raman spectra of Tc. (d) Calculated Franck–Condon fluorescence spectrum, plotted with an inverted x axis to match the normal mode energies. The dashed vertical lines are a guide for the eye.

In the 16 molecule system used above, the LP is higher in energy than the DS0 (in the cavity of Figure , such a situation would correspond to LP states with k > 8.4 μm–1), whereas in experiments VAS is typically observed when the LP state is lower in energy than the dark state manifold.[29,31] We therefore repeated the simulations for a system in which the energy of the LP is below the dark state manifold. In this system, we strongly coupled 512 Tc molecules to a single-mode cavity, resonant with the Tc excitation at 3.22 eV and with a vacuum field strength of 0.000493 au (2.5 MVcm–1, Vcav = 144.1 nm3). In this system, the lowest energy dark state DS0, in which the first Tc molecule is in the S1 minimum geometry while the other 511 molecules are in the S0 geometry, is 39.1 meV higher in energy than the LP. As before, these simulations were initiated in the DS0 state. In each simulation, we again selectively activated one of the 84 vibrational modes of the first Tc molecule, by providing initial atomic velocities in the direction of the displacement of that vibrational mode (see SI for details). The middle panel in Figure b shows the contribution of the cavity mode at 1 fs as a function of the vibrational energy of the vibrational mode that was activated. As in the 16 molecule ensemble, enhanced transfer was observed for vibrational modes that overlap with the nonadiabatic coupling vector, confirming that also when the energy of the LP state is below that of the DS0 state, relaxation into the LP is mediated by specific vibrations, rather than by all vibrations.

Because Raman spectra have been used to identify vibrational modes that can mediate relaxation,[26,29,31] we also computed the Raman spectrum. In Figure c, we compare the Raman intensities of the vibrational modes to the overlap of these modes with the nonadiabatic coupling vector. Because Tc belongs to the D2h point group, both terms in the nonadiabatic coupling vector are fully symmetric (Ag, eq 23 in the SI).[72,73] Therefore, only fully symmetric vibrations of the same Ag symmetry can induce population transfer. Since fully symmetric vibrational modes are also Raman-active, these vibrations appear in the Raman spectrum of Figure c. However, because in Tc vibrational modes of B1g, B2g, and B3g symmetry are Raman-active as well, the Raman spectrum also contains peaks due to vibrations that cannot induce population transfer (e.g., vibrational mode 15 of B1g symmetry at 63 meV), in line with experiments.[29,31]

Whereas for Tc the modes that drive relaxation can be identified based on symmetry arguments, for molecules that lack internal symmetry (i.e., C1), these modes can only be identified by computing their overlap with the nonadiabatic coupling vector. To illustrate this, we have also performed simulations of a single-mode cavity with 16 rhodamine molecules. The results of these simulations, discussed in the SI and shown in Figure S8, suggest that indeed the Raman spectrum does not predict all of the modes that can drive the relaxation between the DS0 and the LP state, whereas these modes are easily identified from their overlap with the nonadiabatic coupling vector.

The nonadiabatic coupling vector contains two terms (eq 20, SI), the first of which is the gradient difference vector (i.e., ∇VS(R) – ∇VS(R)). Therefore, as in previous works,[72,73] the modes that drive the relaxation process could be identified from their overlap with the gradient difference vector as well. Indeed, because the gradient difference vector at the S0 or S1 minima (eq 23 in the SI) describes the displacement between these minima, all modes that are Franck–Condon active can induce population transfer from the DS0 state into the LP. These modes can either be computed or identified from the vibronic progression in an absorption or emission spectrum (Figure d for Tc, Figure S8d for rhodamine). In contrast to the Raman and vibronic spectra, the IR spectrum, also shown in Figure c, is of limited value for predicting vibrational modes that enhance relaxation into the LP for both Tc and rhodamine.

As shown by Pechukas, during nonadiabatic transitions, a semiclassical trajectory evolves under a force in the direction of the nonadiabatic coupling vector.[79,80] The work performed by this force dissipates the energy gap, ΔEDS = ELP – EDS, between the adiabatic states into a displacement.[81] Whereas in our classical MD simulations, population can transfer at any energy gap due to the absence of energy quantization in the vibrational modes, in reality, transfer will only occur if the work performed by the Pechukas’ force equals the energy of a vibrational transition between eigenstates of the vibrational mode along which that force acts. Therefore, in addition to the overlap with the nonadiabatic coupling vector (eq ), population transfer also requires that the energy spacing between the eigenstates of the vibrational mode matches the energy gap between the polaritonic states.

Because the latter aspect of molecular vibrations cannot be easily modeled in classical MD simulations, we performed quantum dynamics simulations of a simple model system instead, in which the LP and DS0 polaritonic surfaces are described by one-dimensional harmonic potentials coupled through the nonadiabatic coupling. The details of this model are presented in the SI. The parameters for the vibrational frequency, energy gap, and nonadiabatic coupling in this model were derived from the simulation with 512 Tc molecules. As shown in Figure S9, if the energy gap, ΔEDS, between the DS0 and LP states is the same as the vibrational energy spacing, ℏωvib, population transfers rapidly from the vibrational ground state in DS0 (|0⟩DS) into the vibrational first excited state in the LP (|1⟩LP). In contrast, if the gap between the LP and DS0 states does not match the energy difference between the vibrational levels, even by as little as a few meV, there is no transfer. On the basis of these results, we speculate that the two requirements for vibrations to mediate population transfer from the dark state manifold into the LP, namely, overlap with the nonadiabatic coupling vector dDS and resonance with ΔEDS, provide a rationale for the luminescence enhancement at specific vibrational frequencies observed experimentally.[26,29,31,56] While previous theoretical analyses had already included a phenomenological decay channel involving vibrational modes,[31,40,41,47,50] the added value of our work is that such modes can actually be predicted from first-principles computations of nonadiabatic coupling vectors and vibrational normal modes.

Finally, we investigated the radiative pumping (RP) mechanism, which has been proposed to dominate the relaxation process if the fluorescence spectrum of the molecule overlaps with the LP branch.[26,57,58] This mechanism is speculated to involve a direct exchange of a photon between an excited molecule and the LP branch. To mimic the initial conditions under which RP is speculated to occur, we modeled a molecule–cavity system with 512 Tc molecules in a periodic one-dimensional cavity (Figure b),[42,54] in which the DS0 state at 2.85 eV is nearly degenerate with several states on the LP branch. The cavity was red-detuned with ℏω0 = 2.69 eV at k = 0 μm–1. To increase the number of LP states with energies near the DS0 state, the dispersion of this cavity, , was modeled with 60 discrete modes, i.e., k = 2πn/Lcav with Lcav = 15 μm the cavity length and 0 ≤ n < 59 the mode index. As before, the 84 short trajectories were initiated with velocities along the displacement vectors of each of the vibrational modes of the Tc molecule in the S1 minimum geometry. To reduce the coupling to the cavity and increase localization of the excitation, this molecule was oriented such that the angle between its transition dipole moment and the cavity vacuum field was 80°.

In contrast to the single-mode cavities, there are many states forming the LP branch in this system. Therefore, to quantify population transfer from the DS0 state into all of these LP states, we summed the projections of all cavity mode excitations (eq ) onto the total polaritonic wave function (eq ), i.e., ∑⟨1|⟨S01S02...S0|Ψ(t)⟩. In the third row of Figure b, we plot the sum of these projections at 1 fs as a function of the vibrational mode energy. As in the single-mode cavity systems, population transfer occurs most efficiently for displacements of vibrational modes that overlap with the nonadiabatic coupling vector dDS (Figure a).

Thus, even though we do not explicitly account for the emission or reabsorption of photons in our simulations, we still observe relaxation from the DS0 state into the LP branch when these states nearly overlap, as required for RP. We also observe that this process is driven by the same vibrational modes as when there is an energy gap between these states. We therefore suggest that the mechanism previously identified as RP[57,58] does not involve a direct exchange of a photon but rather should be considered a form of VAS, which does not require a change in the vibrational eigenstates, because there is no gap between the DS0 and LP states.

In conclusion, the results of our atomistic MD simulations of tetracene molecules strongly coupled to the confined light modes of a Fabry–Pérot cavity suggest that displacements of specific vibrational modes drive nonadiabatic relaxation from the dark state manifold into the lower polariton. These vibrations can be identified by their overlap with the nonadiabatic coupling vector dDS between the lowest-energy dark state DS0, in which the excitation is localized onto a single molecule,[40,41] and the LP, in which the excitation is delocalized over many molecules and cavity modes. The efficiency of this relaxation depends on the energy gap between the lowest energy dark state DS0 and the LP, not only because the magnitude of dDS is inversely proportional to that gap (eq ) but also because the frequencies of the vibrational modes need to match that gap for population transfer to occur. Because it is possible to identify these modes with ab initio calculations, our work paves the way for systematic optimization of relaxation channels in molecule–cavity systems for specific applications. Finally, we could also show that the nonadiabatic coupling vector dDS scales as the inverse of the square root of the number of strongly coupled molecules. Because simulations are typically performed with fewer molecules than in experiments, this scaling allows one to identify finite-size artifacts that must be taken into account when comparing results between simulations and experiments. In the context of our work, the scaling implies that, in the thermodynamic limit with approximately 106 molecules inside the mode volume of a Fabry–Pérot cavity,[51,78] the nonadiabatic coupling vectors would approximately be 50 times smaller than in our simulations with 512 molecules but still be in line with the (sub)picosecond relaxation times measured experimentally.[25,39,82,83]