Simulating molecular polaritons in the collective regime using few-molecule models

Significance Polariton chemistry holds the exciting promise of manipulating matter properties by strong light–matter coupling with optical microcavities. A crucial step toward achieving this goal is to have a tractable theoretical and computational framework that can describe both, the myriad of molecules in the ensemble, and complex internal vibrational structure characteristic of each individual molecule. Our CUT-E formalism demonstrates that a system of only one or two effective molecules strongly coupled to the cavity mode can describe all relevant phenomena, making the system easily solvable using existing quantum dynamics methods. Besides being a powerful computational tool, our formalism provides intuition on how collective effects can affect local excited state dynamics, establishing a framework to design new strategies for polariton chemistry.

The study of molecular polaritons beyond simple quantum emitter ensemble models (e.g., Tavis-Cummings) is challenging due to the large dimensionality of these systems and the complex interplay of molecular electronic and nuclear degrees of freedom. This complexity constrains existing models to either coarse-grain the rich physics and chemistry of the molecular degrees of freedom or artificially limit the description to a small number of molecules. In this work, we exploit permutational symmetries to drastically reduce the computational cost of ab initio quantum dynamics simulations for large N . Furthermore, we discover an emergent hierarchy of timescales present in these systems, that justifies the use of an effective single molecule to approximately capture the dynamics of the entire ensemble, an approximation that becomes exact as N → ∞. We also systematically derive finite N corrections to the dynamics and show that addition of k extra effective molecules is enough to account for phenomena whose rates scale as O(N −k ). Based on this result, we discuss how to seamlessly modify existing single-molecule strong coupling models to describe the dynamics of the corresponding ensemble. We call this approach collective dynamics using truncated equations (CUT-E), benchmark it against well-known results of polariton relaxation rates, and apply it to describe a universal cavity-assisted energy funneling mechanism between different molecular species. Beyond being a computationally efficient tool, this formalism provides an intuitive picture for understanding the role of bright and dark states in chemical reactivity, necessary to generate robust strategies for polariton chemistry.
polariton chemistry | collective strong-coupling regime | permutational symmetries | quantum dynamics Molecular polaritons are quasiparticles arising when vibrational and/or electronic excitations of an ensemble of molecules are collectively coupled to a confined electromagnetic mode such as those found in optical microcavities. These systems have attracted much interest in the last decade due to their potential applications in altering chemical reactions dynamics (1)(2)(3)(4)(5)(6)(7)(8)(9), energy transfer and energy conversion mechanisms (10)(11)(12)(13)(14)(15)(16)(17)(18)(19), and as a framework to achieve room-temperature exciton-polariton condensation (20,21). Theoretical work aimed at explaining experimental results or predicting new phenomena emerging in polaritonic architectures face the formidable challenge of properly modeling the molecular (local) degrees of freedom of each molecule while describing the superradiant interaction of the molecular ensemble with the field (collective). The dynamics arising from the complex interplay of vibrational and electronic degrees of freedom in molecules renders simple quantum optics models such as the original Tavis-Cummings Hamiltonian (22), limited in their applicability to molecular polaritons. Thus, molecular polaritons face unique challenges and opportunities that are not encountered in more traditional polariton systems (23), such as atomic or artificial qubit ensembles (24), or cryogenic inorganic semiconductors (25). Most of the reported simulations of molecular polaritons can only deal with one of two aforementioned challenges at a time. On the one hand, theoretical studies that acknowledge the collective nature of the lightmatter coupling are typically limited to a few dozen molecules at a time and involve sophisticated numerical treatments (26), simplifications such as single vibrational mode descriptions (27), or semiclassical trajectories (8,28). On the other hand, models that implement ab initio treatments are often restricted to a single or few molecules in a cavity (29)(30)(31). Regardless, from a computational standpoint, it seems suspicious that it is necessary to explicitly simulate the dynamics of N molecules, especially if they are identical to each other. Indeed, there are numerous symmetries in the system that should significantly reduce the computational cost of these simulations (32)(33)(34)(35)(36)(37)(38)(39)(40)(41)(42).

Significance
Polariton chemistry holds the exciting promise of manipulating matter properties by strong light-matter coupling with optical microcavities. A crucial step toward achieving this goal is to have a tractable theoretical and computational framework that can describe both, the myriad of molecules in the ensemble, and complex internal vibrational structure characteristic of each individual molecule. Our CUT-E formalism demonstrates that a system of only one or two effective molecules strongly coupled to the cavity mode can describe all relevant phenomena, making the system easily solvable using existing quantum dynamics methods. Besides being a powerful computational tool, our formalism provides intuition on how collective effects can affect local excited state dynamics, establishing a framework to design new strategies for polariton chemistry. Fig. 1. (A) Molecular polaritons in the collective strong-coupling regime. Molecules are not well described by structureless two-level systems, and the interplay between their internal (e.g., vibrational) degrees of freedom and the collective interaction of their optical (e.g., electronic) transitions with the optical mode cannot be described using simple models such as the standard Tavis-Cummings Hamiltonian. (B) Pictorial representation of the collective dynamics using truncated equations (CUT-E) method. Processes with O(N −k ) rates can be described by a model of k + 1 effective molecules coupled to a cavity mode.

A B
In this work, we outline a wavefunction-based formalism that makes use of such symmetries to significantly reduce the complexity of quantum dynamics simulations of the single-excitation manifold of molecular polaritons. Moreover, this formalism naturally provides a hierarchy of approximations to further simplify the problem in a way that, in the N → ∞ limit, polaritonic properties can be calculated using a modified effective single molecule coupled to a cavity with the collective coupling. This provides grounds for some single-molecule strong coupling phenomena to appear in the collective regime, consistent with previous work where linear optical properties can be calculated from effective single-molecule models in the thermodynamic limit (43,44). Moreover, we show that a system of two effective molecules can describe all the effects with rates that scale as 1/N . In general, we show that processes with O(N −k ) rates are described by k + 1 effective molecules strongly coupled to the cavity. This implies that, for a large ensemble of molecules, it is enough to consider only a few effective molecules to solve for the dynamics of the original polariton system (Fig. 1). The model can be applied to study disordered ensembles (e.g., a mixture with two chemical species) without a significant increase in the computational cost.
The article is organized as follows: In Section 1, we present the Hamiltonian and the multiconfigurational representation of the total wavefunction of the system, in which permutational symmetries become evident. Then, we uncover a convenient mathematical structure of the equations of motion (EoM) where approximate symmetries emerge, and which become useful for large N , while keeping the collective light-matter coupling √ N g finite, with g being the single-molecule coupling (this is the physical condition of interest in experiments and which concerns us hereafter). This structure allows us to derive the simple effective Hamiltonians involving only a few molecules, that solve for the dynamics of the entire ensemble. In Section 2, we make use of the effective single-molecule model to demonstrate how both optical and material properties in the original system can be computed using the effective singlemolecule simulation. In Section 3, we benchmark our formalism against a well-known result: the nonradiative relaxation of polariton and dark states. In Section 4, we present a pedagogical and intriguing application that reveals the power of our formalism, describing how to exploit polariton dynamics to obtain nonstatistical outcomes in photoproducts. Finally, we summarize the work in Section 5.

Theory
A. Hamiltonian and Multiconfigurational Wavefunction. Consider a system of N molecules collectively coupled to a single cavity mode. The Tavis-Cummings Hamiltonian, extended to include vibrational degrees of freedom missing from original models, can be written as (hereafterh = 1)Ĥ whereĤ (i) are the Hamiltonians for the ith molecule, the cavity mode, and the interaction between them. Here, µ is the reduced mass of the nuclei, |g i and |e i are the molecular ground and excited electronic states, V g/e (q i ) are the ground and excited Potential Energy Surfaces (PES),â is the photon annihilation operator, and q i is the vector of all molecular vibrations of molecule i (Fig. 1). When the PESs are harmonic, the model above reduces to the Holstein-Tavis Cummings model, which has been subject of recent studies (27,43). For the time being, we assume there is no disorder, neglect intermolecular interactions, and use the rotating wave approximation by considering the single-molecule coupling strength g to be much smaller than the bare photon frequency ω c . Finally, we also ignore (nonradiative) couplings between the molecular ground and excited states whose PESs can form conical intersections (45)(46)(47). Under these hypotheses, the excitation number [sum of electronic excitations (Frenkel excitons) and photon number] is conserved.
In this article, we shall focus on the so-called first-excitation manifold, which describes processes in which the cavity triggers excited state dynamics of one molecule at a time. The position-representation ansatz takes the form with the electron-photon states |e i = |g 1 , g 2 , ..., e i , ..., g N , 0 ph and |1 = |g 1 , g 2 , ..., g N , 1 ph , and the vibrational wavefunctions In this expansion, multiconfigurational vibrational wavefunctions ψ are built with sets of m single particle orthonormal functions ϕ j i (q i ) and φ j i (q i ), respectively, that are equal for identical molecules. In this work, we choose such functions to be the eigenstates of the molecular HamiltonianĤ Additionally, there are permutations between coefficients of different excitons, given the interaction of the molecules with the cavity is assumed identical, This means that every coefficient of the electronic state i can be obtained directly from those in the electronic state i. In other words, Eqs. 4 and 5 yield the crucial observation that it is enough to calculate the dynamics of a single excitonic state to know the evolution of all of them. Moreover, the vibrational states of the ensemble of molecules can be completely characterized by specifying the number of molecules in each vibrational state; therefore, the number of vibrational degrees of freedom also reduces drastically. The ground and excited state coefficients can be written in terms of permutationally symmetric states, where N k is the number of ground-state molecules in the vibrational state k. This new notation removes the information about the state of each specific molecule (A configurations. This corresponds to reducing the complexity from an exponential scaling to a polynomic one, as has been shown in previous work that also exploit permutational symmetries for studying ensembles of identical systems (42,48). We will make use of these simplifications to redefine the EoM and subsequently the Hamiltonians. Such analysis is general and can also be done for higher excitation manifolds. In that case, the number of coefficients increases with the number of excitations due to permutations between electronically excited molecules; however, after the number of excitations is half of the number of molecules, this trend reverses.
The EoM for the coefficients can be obtained by using the Dirac-Frenkel variational principle or simply by inserting the ansatz wavefuntion into the time-dependent Schrödinger equation. Rewriting the ansatz as Since Eq. 1 ignores couplings between the molecules in the absence of the photon mode, With these considerations, Eq. 7 becomes The last equation can be simplified using Eq. 6 to write the dynamics in terms of permutationally symmetric states: [10] The first equation represents the absorption of the photon that takes a molecule from the state ϕ k |g into the state φ l |e . The second equation describes the conjugate process.
C. Structure of the Wavefunction. Even though Eq. 10 are exact and represent a significant improvement over Eq. 9, they also allow us to systematically introduce approximations by virtue of the factors N k , which represent the number of ground state molecules in the vibrational state k. For initial states in which one of the values N k is exceptionally large, the dynamics is such that N k is almost conserved, as discussed below. As an example, assume we start with all the molecules in their ground vibrational state and 1 photon in the cavity mode, i.e., A N 00···0 (0) = 1. Thus, we can simplify even more our notation by only reporting the number of such vibrational excitations. By renormalizing the coefficients with N -dependent factors that indicate the number of states that correspond to the same amplitude, we can recover the basis originally introduced by Spano for harmonic modes but which we now use for arbitrary PESs 33: With this final notation, the EoM are written in Eq. 12, where we have removed a constant NE g,1 and defined ω eg,l = E e,l − E g, 1 and In Fig. 2, we provide a pictorial representation of the states associated with these coefficients.

[12]
The first line of Eq. 12 reveals that the initial photonic stateÃ   k (t) upon emission. The first of these two processes depends on the collective light-matter coupling √ N g, while the second one depends on the single-molecule light-matter coupling g. This structure is repeated throughout the system of equations: coupling between states conserving the number of molecules with phonons is collective, while processes that increase the number of such molecules are proportional to the single-molecule coupling. It is interesting to note that this phenomenon is well known in the literature of molecular aggregates. In particular, for J-aggregates, Spano and Yamagata (35) have noted that the ratio of the photoluminescence into the electronic ground state with no phonons versus that into the electronic ground state with one phonon is proportional to the coherence length N of the aggregate. While this phenomenon is routinely used as a spectroscopic probe for N , we use it in our case to drastically simplify the simulations of molecular polaritons, as explained below. Ng couplings). These fast dynamics are linked by bottlenecks (due to single-molecule g couplings) which slowly change the number of molecules featuring such ground-state phonons. Zeroth-order approximation in g corresponds to restricting the dynamics to states with a fixed number of ground-state molecules with phonons, while adding the first-order correction allows the dynamics to create (or annihilate) phonons in (from) 1 additional molecule.
Assuming that we have a large number of molecules, we can treat the mixing of states with differing numbers of ground-state molecules with phonons perturbatively, with the single-molecule light-matter coupling g as the perturbation. This is schematically represented in Fig. 3.
In the limit where g → 0 or N → ∞ (while keeping √ N g constant), which is the regime that interests us, the number of ground-state molecules that support vibrational excitations is conserved during the dynamics. For our initial state, the zeroth-order approximation implies that the wavefunction is described only by the basis states of the left most block of Fig. 3. Adding the first-order correction allows states in the immediate next block (featuring only one ground-state molecule with phonons) to contribute to the wavefunction. The timescale at which the first-order terms contribute is much longer than the ultrafast vibrational dynamics on each excited molecule. The exact wavefunction is recovered as one spans the entire Hilbert space from Left to Right in Fig. 3, but as can be appreciated, there are only a few molecules with phonons in the electronic ground state for large N , even for long times of interest, justifying the convenience of the shorthand notation in Eq. 11.
Although the fact that there is only one molecular specie in the system that is central for the permutational symmetries of Eqs. 4 and 5 to hold, addition of different molecular species can be done without dramatically increasing the computational cost since permutational symmetries still apply for each type. The renormalization of the coefficients associated with each species as in Eq. 11 will depend on their concentration. Similarly, disorder due to inhomogeneous broadening or spacial variation of the coupling to the photon mode in multimode cavities, which has been shown to play a central role in molecular polaritonics systems (50)(51)(52), can be included by adding new molecular species for each value that is sampled according to the corresponding distribution of excitonic frequencies or interaction strength. The number of molecular species (disorder bins N bins ) required to simulate a disordered ensemble with frequency width W can be roughly estimated as N bins ≈ W /(2π /T ), where T is the total simulation time that allows to capture the physical phenomena of interest. We believe that this estimate is an upper-bound, as it does not take into account efficient ways to simulate disorder (for instance, ref. 53), but future work will aim to study these subtleties. D. Zeroth-Order Approximation. Let us assume that temperature T = 0, meaning our initial state in the original notation is given by A (0) 111...1 (0) = 1. The time-dependent wavefunction (Eq. 2) at the zeroth-order approximation is given by the following vibrational wavefunctions [13] By using the renormalized permutationally symmetric coefficients, the EoM in Eq. 12 transform into those of an effective singlemolecule strongly coupled to the cavity mode (Eq. 14).
Eq. 14 is consistent with previous results where an impurity (in this case the optical mode) coupled to a large environment can be simplified into an impurity interacting with an effective harmonic bath that includes the relevant frequencies of the environment (54), in this case, the optical transitions of an individual molecule. Alternatively, it is also consistent with the classical optics treatments of polaritons arising as the result of a photonic oscillator coupling to a set of effective oscillators representing the molecular transitions (e.g., transfer matrix methods) (55,56). This has also been justified by Keeling and coworkers within a quantum mechanical framework (43,44). We can rewrite Eq. 14 in a form that is better suited for implementation in quantum dynamics packages such as the multiconfiguration time-dependent Hartree (MCTDH) method (49,57,58), , [15] with the projector over the ground vibrational state P = |ϕ 1 ϕ 1 | and the photonic and excitonic wavefunctions Here, we readily identify the zeroth-order Hamiltonian in Eq. 15 aŝ withĤ g/e = 1 2µp 2 + V g/e (q). Eq. 15 can be used to recover the original time-dependent many-body wavefunctions in Eq. 13, As far as we are aware, the separation between collective and single-molecule emission processes was first pointed out by Spano when calculating vibrationally dressed lower polariton states (34), although crucial ideas were introduced by the same author much earlier (35), as we shall discuss later. In a more recently article, the same author combined such ideas with the permutationally symmetric basis to compute photoluminescent spectra starting from approximated vibro-polaritonic eigenstates (33). Our work formalizes Spano's observations into the hierarchy summarized in Fig. 3 and capitalizes it to compute dynamics in complex molecular polaritonic systems featuring arbitrary PESs. Three important comments are in order concerning Eq. 17. First, the artificial effective molecule in the ground state is only allowed to be in its ground vibrational state (no phonons). The physical intuition for this is simple: Imagine the molecules have been collectively excited at their Franck-Condon (FC) configurations by the cavity field; each of them can in principle reemit such energy into the cavity creating any vibrational state of their ground electronic state. However, only when they go back to their vibrational ground state, the number of ground state molecules with phonons is conserved and the process becomes collective. This seemingly unremarkable observation has an important implication: it states that single-molecule polariton simulations are not applicable to the collective regime unless the light-matter coupling is restricted to the FC region, where the vibrational ground state ϕ 1 (q) has a significant amplitude (Fig. 4). However, as Eq. 15 reveals, such simulations can be simply adjusted for the collective regime by the appropriate inclusion of the projector P. This projector can be readily implemented in various ways (e.g., by restricting the vibrational basis in the ground state to |ϕ 1 or by projecting the light-matter coupling, among several options). Notice that Eq. 17 does not formally lead to the definition of polaritonic PESs that govern the nuclear dynamics in the presence of polaritons as in ref. 9. However, the key message in the latter work is still consistent with ours: Nuclei feel different forces due to collective effects only near the Franck-Condon configuration. Second, since our formalism explicitly considers the vibrational degrees of freedom, the simplification to a single effective molecule does not imply the elimination of the dark-state manifold from the Hilbert space, as we will explain in detail in Section 2B.1. Third, the results of this section are exact in the thermodynamic (N → ∞) limit because they imply g → 0. We next add the first-order correction, where the wavefunction now has states with amplitude of the order of g. [19] The general case can be written in a compact form using projection operators, where Q i = 1 vib,i − P i (Eq. 20). The first 2 × 2 block corresponds to the zeroth-order approximation and describes the aforementioned effective molecule (say, molecule 1) coupled to the N − 1g. This effective light-matter interaction nonlinearity, which conditions the light-matter interaction of the second molecule on the presence of phonons in the first one, might look odd at first sight but must be endorsed as the emergent physics resulting from the hierarchy. Eq. 20 is in a suitable form for implementation in wavepacket dynamics methods.
Notice that ignoring the matrix elements proportional to g (boxed) amounts to block-diagonalizing the Hamiltonian according to the number of electronic ground-state molecules with phonons, reiterating the approximate symmetry featured in Fig. 3. Eq. 20 can be interpreted as two effective molecules coupled to the cavity. Finally, the original many-body wavefunction in Eq. 3 can be rewritten in terms of the 2-molecule wavepackets as with i as an arbitrary molecule of the ensemble and lk (t)φ l (q 2 )ϕ k (q 1 ).

[22]
We make use of this Hamiltonian in Section 3 and show that adding the first-order correction describes collective processes, in addition to those whose rates scale as 1/N . Similarly, corrections of kth order require using a Hamiltonian of k + 1 effective molecules interacting with a cavity mode and would include information about phenomena with rates that scale up to 1/N k . Such decomposition of the many-molecules problem into a few-molecules problem is what we call collective dynamics using truncated equations (CUT-E) method. Let us summarize the central results of the article. Application of the permutational symmetries in Eqs. 4 and 5 to the time-dependent Schrödinger equation for N identical molecules in a cavity gives rise to the EoM in Eq. 12. These EoM are endowed with a convenient hierarchy of timescale separations, schematically illustrated by the scissors in Fig. 3. For a fixed collective coupling √ N g, processes with N -independent (O(N −1 )) rates are captured by a single (two) effective molecule(s) in a cavity, according to Eq. 15 (Eq. 20). This is schematically represented in Fig. 1B. Although high-order corrections might become relevant when a small number of molecules is considered, here we are interested in the case where a large number of molecules is used to reach collective strong light-matter coupling. Thus, we will defer the exploration of corrections beyond first order for future works.

Observables in the Zeroth-Order Approximation
Observables of the real system must be calculated using the many-body wavefunction | (t) and not directly using the effective single-molecule wavefunction |˜ (t) . However, there are particular observables for which both wavefunctions provide the same answer. This will depend on whether the observable is local or collective. For pedagogical purposes, throughout this section, we assume the zeroth-order approximation.
A. Local Properties. Hereafter, we define local properties to be quantities that only depend on one single molecule (e.g., a chemical reaction) or the photon field alone. A.1. Chemical properties. Consider a local observableˆ (i) that depends only on degrees of freedom of molecule i. As an example, assume we are interested in the nuclear dynamics in the excited stateˆ (i) =q i |e i e i |. Using the zeroth-order wavefunction in Eq. 13, we obtain This is a remarkable result because it demonstrates that the excited-state dynamics of N molecules collectively coupled to a cavity mode (described by | (t) ) is identical to the dynamics of an effective single molecule strongly coupled to a cavity (described by |˜ (t) ), except for a constant 1/N dilution factor. This factor is just a consequence of using a single photon to alter the excited-state dynamics of N molecules. The probability of any molecule being at the configuration q is the same as that of the single effective molecule strongly coupled to the cavity. However, we emphasize that the effective molecule is not allowed to have phonons while it is in the ground state, a restriction which can introduce significant differences compared to standard single-molecule calculations. An important corollary of the above analysis is that effects predicted using single-molecule models might actually occur in the collective regime if they rely on changes in the excited PES at the FC region but not the ones relying on changes beyond. This was first pointed out by Galego et al. (9) by analyzing polaritonic PESs for an ensemble of molecules interacting with a cavity mode. This is also consistent with a recent work by Cui and Nitzan (59), who concluded that excited-state dynamics in polaritonic systems is dominated by states that are reachable from the ground electronic state. Yet, we propose that modifications of chemical dynamics that occur on a timescale longer than the decay of the initially prepared excitations can occur in an N -independent manner if they dramatically depend on the ultrafast polariton-modified dynamics at the FC region. In Section 4, we present such an example. A.2. Optical properties. Another set of properties that are equivalent in the single molecule and N → ∞ cases are those that can be extracted only from the dynamics of the field. For example, the linear transmission, absorption, and reflection spectra can be calculated from the photon autocorrelation function c(t). In the zeroth-order approximation, it can be shown that c(t) = (0)| (t) = ˜ (0)|˜ (t) , where |˜ (0) = ϕ 0 (q)|1 is the photonic state. This is consistent with previous work by the Keeling group (43,44).
B. Nonlocal Properties. Another set of observables consists of operators that are delocalized across several molecules or depend on both molecular and optical degrees of freedom. Some of the most common observables of this kind are the populations of the polariton and dark states. Such observables can be obtained with the effective single-molecule model but not from its reduced electronic-photonic density matrix, as we will show next.

B.1. Polariton and dark-states populations.
Let us write the expectation value of an arbitrary electronic-photonic operator ˆ = Tr ρˆ in terms of the reduced density matrix of the effective single-molecule systemρ, where we have identified |˜ FC = ϕ 1 (q)|e as the FC wavepacket. The above equation implies that the reduced density matrixρ of the effective single molecule is, in general, not enough to calculate any delocalized molecular observables since the last term of Eq. 24 refers to the projection of the wavefunction onto the FC wavepacket, which corresponds to the interexciton coherences: |e i is the totally symmetric excitonic state) being one of the polariton states. Using Eq. 24, population of this state is calculated to be Similarly, dark states can be written as |D = i c i |e i , with |c 1 | = |c 2 | = · · · = |c N | = 1/ √ N , N i=1 c i = 0, and e i e j = e i |D D|e j = c i c * j i.e., they are chosen orthogonal to the |P states, so here we take them to be the Fourier combinations of excitons that are orthogonal to |B ; see for instance (60). This leads to Notice that this calculation is identical for every dark state in the chosen basis; therefore, the population in the dark-state manifold yields where we have made N −1 N ≈ 1 since the zeroth-order approximation becomes exact for N → ∞. From this result, we can extract the intuitive interpretation that the bright state |B in the N -molecule system corresponds to a FC wavepacket in the excited state of the effective single molecule, while the dark states correspond to the rest of the wavefunction whose nuclear configurations lie outside of the FC region (Fig. 4).
For cases where the chemically relevant excited-state dynamics involves fast changes of the nuclear configurations away from the FC point, excited-state reactivity is essentially relaxation to dark states. The question is whether the relaxation occurs along the reactive coordinate of interest (e.g., particular reactive dark-modes) or whether it occurs along orthogonal modes to it. Thus, we have derived a powerful design principle for polariton chemistry which has so far unjustifiably gathered little attention: The strategy is not to avoid decay into dark states, which seems inexorable in most cases, but to use strong coupling to control which dark states to target. In fact, using the zeroth-order approximation, we will illustrate some mechanisms to manipulate ratios for these relaxation pathways in Section 4.

Polariton Vibrational Relaxation
Previous work by del Pino et al. (61) addressed the relaxation dynamics of molecular polaritons using the Hamiltonian in Eq. 1 for a vibrational harmonic bath and linear vibronic coupling, [29] This model of relaxation is the single-photon mode simplification of a previous model by Litinskaya and Agranovich (62). In the limit where vibronic coupling is much smaller than the collective light-matter coupling and there is an energy-dense set of vibrational modes k, we can use our formalism to analytically derive relaxation rates by including a collection of local vibrational modes and using Fermi's Golden Rule with the vibronic couplings as the perturbation. Although such a result is already well known, it serves as a benchmark and illustration for our formalism. To unclutter the calculations, we consider the case where the excitons and the cavity are in resonance (ω c = ω eg,1 = ω). The calculations are explained in detail in supplementary information, and the results are summarized in Table 1. Apart from missing a factor of N −1 N in D←+ , the zeroth-order approximation correctly predicts the relaxation rate from upper polariton to dark states (which in this weak-vibronic coupling model in the effective single-molecule calculation, corresponds to the bare molecular exciton with a single phonon |e, 1 ) but disregards both the upper to lower polariton and dark state to lower polariton Table 1

Zeroth-order approximation
With first-order correction rates (61). This makes sense since the latter two are known to be proportional to 1/N , and the zeroth-order approximation is exact for N → ∞. Thus, for this simple model, the only relaxation process of relevance at ultrafast timescales is the downhill one into dark states. Addition of the first-order correction recovers the decay rates −←− and −←D , as well as the N −1 N factor missing from D←+ in the zero-order approximation. We suspect the factor N −1 N 2 in −←D will become exactly 1/N if the second-order correction is considered. It is unlikely that we would need to go to higher-order corrections for microcavity polaritons or even for polaritons arising in plasmonic antennas, where N = 100 to 1,000 (63), although if strong coupling can be demonstrated with smaller N values, those corrections could start mattering. On the other hand, it would be of interest to compare the performance of our method with the standard one of explicitly simulating several molecules using an optimized algorithm such as multilayer MCTDH (64) and characterize the values of N after which our method can outcompete the latter.
A. Finite Temperature Effects. In the previous sections, we have only dealt with transitions from higher to lower lying polaritonic or dark states. To calculate rates such as D←− , we need to consider all possible initial states allowed by thermal fluctuations. In general, such states involve breaking of the symmetry that is essential for the reduction of the dimensionality in Eq. 10. Therefore, at the current stage, this formalism cannot be easily generalized to finite temperatures. Future works will focus on developing a density matrix approach in which permutational symmetries can be smoothly applied.

Nonstatistical Excited-State Dynamics
In this section, we look at the mechanism whereby, given two molecular species strongly coupled to a cavity, excitation energy can be selectively funneled to one of the species. We will show that statistical yield estimates based on linear optical spectroscopy are misleading at predicting these outcomes accurately.
For a system with N A molecules of species A and N B molecules of species B inside a cavity, the Hamiltonian in the zeroth-order approximation can be readily generalized from Eq. 14. The Hamiltonian and parameters of the PESs are reported in SI Appendix. In example 1, the excited PES of species A is a displaced harmonic oscillator and that of species B is a dissociative potential. While the latter is mostly relevant for gas phase reactions, a completely analogous photochemical phenomenon in the condensed phases could be obtained with diabatically coupled harmonic surfaces. The cavity frequency is resonant with the FC transitions of both A and B (Fig. 5, 1A).
Given the one-dimensional nature of the PESs, it is straightforward to explicitly construct the eigenstates and eigenenergies of the vibrational Hamiltonians and perform the short-time unitary dynamics; we use the standard discrete variable representation method by Colbert and Miller (65). This calculation was then used to compute spectra of the molecules outside and in the cavity using the formalism illustrated in refs. 32 and 43. The parameters for calculation of the spectra are chosen such that the frequency resolution and total simulation time satisfy Fourier transform relations. The total simulation time of 2.5 ps determines the molecular and cavity linewidths as 0.002 eV (Fig. 5, 1B, C, and D). The absorption spectra of the bare molecules reveal the energy level structure accessed at the respective FC regions (66). For molecules A, we see the strongest peak corresponding to the FC transition accompanied by the other peaks of the vibronic progression. On the other hand, for molecules B, we see a single broad feature due to the dissociative potential. When strongly coupled to a cavity, these peaks form a rich pattern of peak splittings which can be intuitively understood upon diagonalization of Eq. 14 (for instance, for molecule A, six sharp resonances become seven due to coupling to the cavity). We build time-dependent wavefunctions by constructing the corresponding linear combinations of numerically computed eigenstates of the zeroth-order Hamiltonian. This procedure would obviously be impractical for realistic molecular species with many vibrational modes, in which case, an explicit time-dependent approach such as MCTDH would be preferred. This simulation illustrates a phenomenon where the energy initially given to the cavity is eventually channeled preferentially to one of the two molecules. The results of this simulation are consistent with the phenomenology originally theoretically proposed by Groenhof and Toppari, based on computational simulations with at most 1,000 molecules of one of the species (67), and demonstrate the latter remains valid in the thermodynamic limit. For our simulation, we start with an excitation in the cavity, | (0) = ϕ 1 (q A )ϕ 1 (q B )|1 , and calculate excited-state populations for each of the effective molecules, which provides the excited state reactivity of a single molecule of type A/B in the entire molecular ensemble (because we are working in the first excitation manifold). At short times, the cavity evenly distributes the energy into both types of molecules, exciting them at their respective FC regions. We monitor the excited-state populations as a function of time (Fig. 6A). We observe that at short times, the population is transferred equally to both species. However, as time progresses, the excitation is funneled selectively to species B via the cavity. The mechanism of the phenomenon is the following: at short times of the order of the Rabi oscillations, the cavity excites the bright modes of both species with equal populations at the respective FC regions due to the chosen equal collective light-matter couplings. However, as time progresses, excitations in species B decay irreversibly into their dark states upon evolution away from the FC region (Eq. 28) due to the dissociative character of V e0,B (q B ). Molecular species A does not undergo this fast dephasing due to the bound nature of V e0,A (q A ) and is the one that predominantly remains at its FC region and emits a photon to the cavity at the end of the Rabi cycle. The mechanism restarts with the reexcitation of both molecules by the photon. What results from this simulation is a net energy flow from molecular species A to species B mediated by the cavity and which cannot be explained by the contribution of each molecule to the polariton states defined at the FC region. Notice that if the cavity decay is faster than the molecular dephasing, the excitation can leak out before energy transfer from A to B ensues. This issue can be overcome by having the cavity mode under continuous pumping.
Next, we compare the computed long-time populations accumulated in the excited states of molecules A and B to the infinite-time probabilities of molecules of the cavity exciting the different species (p s j for j ∈ {A, B}), Here, {|n } represents the eigenstates ofĤ (0) . Eq. 30 computes the composition of probabilities of an eigenstate simultaneously having photonic and molecule A(B) contribution. However, the entire set of eigenstates is not obtainable using linear absorption spectroscopy alone, given that the latter is a limited projection of information at the FC region (66). As an extreme example, we consider "short-time-resolution eigenstates" consisting of the upper, middle, and lower polaritons computed at the FC region and compute the analogous estimates p s A (FC) and p s B (FC). The infinite and short-time estimates are presented in Fig. 6A. As expected, the plot suggests that the statistical estimate obtained from the full set of eigenvectors gives the right prediction of the yields at long times, while the FC estimate predicts the incorrect outcome of equal populations of A and B. Example 2. Example 1 highlights the danger of inferring excited-state dynamics relying solely on information obtained from linear optical spectroscopy. This issue becomes even more pertinent owing to some observations of Xiang et al. (16), where the authors observe selective cavity-mediated energy transfer into one of two molecules, despite the bare linear absorption spectral lineshapes being very similar outside the cavity. To highlight some of the described subtleties, we explore a second example where we consider a slightly different shape of the excited state PES for species B such that it resembles that of molecular species A in the FC region, while still being dissociative. Energy funneling from species A into species B is observed in both cases, indicating that this phenomenon is not easily predicted from linear optical spectra alone, since at least for Example 2, the absorption spectra of bare molecules A and B are pretty much identical (Fig. 5, 2B and C).

A B
The computed spectra for the bare species and their corresponding spectra under strong coupling to the cavity mode are shown in Fig. 5, 2B, C, and D. The resemblance of the excited state PESs at the FC region translates into very similar spectral lineshapes despite the dissociative vs. bound nature away from the FC region. Regardless of the similarities in the absorption spectra, starting the dynamics with an excitation in the cavity mode still gives rise to the effective energy transfer from species A to B. Accordingly, the statistical ratio computed from the eigenstates at the FC point predicts equal populations of molecular species, while the full set of vibro-polaritonic eigenstates of the system (p s A and p s B in Fig.6B) predicts the correct yields. Since such states are not easily accessible from linear absorption spectroscopy, we consider those measurements insufficient to predict the outcome of these photoproducts.
It is important to note that by virtue of our reliance on the zeroth-order approximation, the cavity-mediated energy transfer mechanism studied in this section is N -independent, and thus differs from that put forward in refs. 68 and 69. The latter relies on a bottleneck transfer of population from dark to polariton modes, whose rate scales as 1/N .
In the future, it will be of great interest to ascertain whether the O(N 0 ) rate mechanism studied in the present work originally put forward in ref. 67 or the mentioned O(N −1 ) rate mechanism is in order in the various experiments of cavity-mediated energy transfer (10,11,16). At least, for the last reference, it is clear that the latter mechanism cannot be operative, for it would yield rates of energy transfer that are much slower than those observed in the experiment.
Finally, we speculate that the mechanism highlighted in this section, where the cavity mediates energy transfer into the fastest dephasing molecular species, should be quite a generic phenomenon and might be at play in the recent study in ref. 70, which shows polariton-enhanced photoconductivity of an organic film. However, a further consideration of this problem merit a careful inclusion of disorder and temperature effects, which will be the focus of future works.

Summary
The method of collective dynamics using truncated equations (CUT-E) exploits permutational symmetries of an ensemble of identical molecules and an emergent hierarchy of timescales to elucidate the excited-state dynamics under collective strong light-matter interaction. Although previous works have used permutational symmetry arguments and have arrived at conclusions that are consistent with ours (31-33, 43, 59), elements of the present work uniquely provide opportunities for a systematic and intuitive treatment of molecular polariton dynamics.
It is important to note that, via quite a different formalism, another method that maps the dynamics of molecular polaritons to a single effective molecule in a cavity has already recently been reported before our work (32). This method, based on density matrices, does not provide O(1/ √ N ) corrections to dynamics and considers a coherent state of the photon at all times, instead of the single-excitation manifold dynamics we have presented. We believe this method is quite complementary in its scope to ours. However, given the different formalisms, it is at present hard to assess the deeper conceptual connections between the methods; this will be subject of future work.
Let us conclude by highlighting some of the features of our approach. First, it allows for the systematic introduction of corrections to the thermodynamic limit. Second, our time-dependent approach generalizes some of the results found in previous work that also exploit permutational symmetries to compute optical properties (32-35, 39, 43); here, we have generalized these concepts to chemical dynamics. Moreover, our work naturally provides an alternative interpretation of bright and dark states based on permutationally symmetric states, that is different from the one inherited from restrictive quantum optics models. These observations provide muchneeded physical intuition to design principles of polariton chemistry control, where rather than avoiding the decay into dark states, one embraces such phenomenon at one's advantage, such as with the examples provided in Section 4. Finally, the method enjoys numerical simplicity and is written in a language that makes it straightforward to modify existing codes for single-molecule strong light-matter calculations, and more generally, convenient for implementation in existing quantum molecular dynamics algorithms. At present, it is unclear whether the strategies here presented can be generalized to multimode cavities and finite temperature scenarios. These issues are subjects of present investigations. Data, Materials, and Software Availability. All study data are included in the article and/or SI Appendix.