Thermalization rate of polaritons in strongly-coupled molecular systems

Polariton thermalization is a key process in achieving light-matter Bose--Einstein condensation, spanning from solid-state semiconductor microcavities at cryogenic temperatures to surface plasmon nanocavities with molecules at room temperature. Originated from the matter component of polariton states, the microscopic mechanisms of thermalization are closely tied to specific material properties. In this work, we investigate polariton thermalization in strongly-coupled molecular systems. We develop a microscopic theory addressing polariton thermalization through electron-phonon interactions (known as exciton-vibration coupling) with low-energy molecular vibrations. This theory presents a simple analytical method to calculate the temperature-dependent polariton thermalization rate, utilizing experimentally accessible spectral properties of bare molecules, such as the Stokes shift and temperature-dependent linewidth of photoluminescence, in conjunction with well-known parameters of optical cavities. Our findings demonstrate qualitative agreement with recent experimental reports of nonequilibrium polariton condensation in both ground and excited states, and explain the thermalization bottleneck effect observed at low temperatures. This study showcases the significance of vibrational degrees of freedom in polariton condensation and offers practical guidance for future experiments, including the selection of suitable material systems and cavity designs.


I. INTRODUCTION
Electronic and vibrational states hold a central place in the molecular systems that drive photo-induced processes in nature [1] and underlie many of the technologies we interact with on a daily basis [2].Despite the large energetic disparity between electronic and vibrational states, they can undergo substantial electronphonon type interactions, known as exciton-vibration coupling [3,4].This coupling defines absorption and emission spectra of molecular systems [5][6][7] and shapes the relaxation dynamics at the microscopic level [8][9][10].When placed inside optical cavities, molecules can engage in strong light-matter interactions, leading to an effective tripartite interaction [11][12][13].This results in the formation of new eigenstates, known as vibronic-type excitonpolaritons [13,14], also referred as polaron polaritons in exciton-plasmon systems [15,16].Due to the typically strong exciton-vibration coupling, these polariton states are prevalent in molecular systems.Being hybrid lightmatter states, polaritons inherit properties from both molecules and the cavity electromagnetic field.From the latter, they acquire a small effective mass and a typically short lifetime ∼ 100 fs.
As bosons, polaritons can form Bose-Einstein condensate (BEC), where their small effective mass and low density of states favor condensation at elevatedtemperatures [17].However, their short lifetime means that polariton BEC in molecular systems is significantly out of thermal equilibrium, requiring constant re-population of the reservoir to compensate for polariton losses [18].Despite the intrinsic nonequilibrium nature of polariton condensates, they can still achieve some form of effective local equilibrium with the environment.This equilibrium is defined by the complex interplay of gain, loss, and thermalization within the polariton system.Recent theoretical reports have rigorously demonstrated this regime through an exact solution for the density matrix in the fast thermalization limit [19,20].The formation of po-lariton BEC, in general, requires two conditions: 1 -the effective rate of polariton thermalization overcomes the energy dissipation, and 2 -the total number of lower polaritons surpasses the critical number [19,20].While the effective thermalization rate overcomes the dissipation rate for the lower polaritons above the condensation threshold, it is not always true at the onset of polariton BEC.Polariton distributions below and just at the condensation threshold typically demonstrate higher temperature than the lattice [15,[21][22][23] which is due to insufficient thermalization rate with respect to fast polariton decay.Therefore, understanding polariton thermalization is pivotal for the physics of light-matter condensates in molecular systems.The actual microscopic mechanisms behind polariton thermalization remain an open question in the field.The current status is mostly relied on multimode mean-field theories, where thermalization is introduced through the Lindblad master equation or similar using effective thermalization constants [12,24].Thermalization rates are typically defined by comparison with experimental data as the best-fit parameters.Recently, Satapathy et al. have shed light on the nature of thermalization in organic microcavities, experimentally demonstrating the central role of the emission Stokes shift of molecules in polariton thermalization towards the BEC [23].However, the microscopic mechanisms of polariton thermalization in molecular systems remain largely unexplored.
The material properties of molecular systems significantly influence thermalization behavior of polaritons [15,21,23,25].In exploring the microscopic origins of thermalization, one must consider the typical energy scale between nearest neighbor polariton states ∼ 0.1 meV.Indeed, the finite size of a realistic system leads to the quantization of polariton states separating them in energy by ∆ω min ∼ S −1 [17] (see Section 3 for details).This energy scale is inconsistent with photon emission from excited states (∼ 1000 meV) and with strong high-energy molecular vibrations (∼ 100 meV).Although coupling to high-energy molecular vibrations enables an efficient energy relaxation mechanism in organic polariton systems [9] driving them towards polariton BEC [12,26], it is unlikely to be the thermalization mechanism within the lower polariton branch.However, in addition to high-energy vibrations, molecules, especially in densely packed molecular layers, such as polymer films, exhibit a wide range of low-energy vibrational modes [27,28].The same type of exciton-vibration coupling between electronic states and low-energy molecular vibrations can bridge nearest neighbor polariton states, matching this energy scale very efficiently.Similar mechanism of electron-acoustic phonon interaction is known to be important for polariton thermalization in crystalline semiconductor microcavities bearing Wannier-Mott type excitons [17,29].
Direct observation of low-frequency vibrational modes below 10 meV (≲ 100 cm −1 ) in most spectroscopic experiments is challenging.Nonetheless, these modes are omnipresent and influence the spectral properties of molecules implicitly as well as the dynamics of excited states [30,31].This influence is evident in phenomena such as the Stokes shift of the 0 − 0 vibronic emission peak relative to absorption [31] or temperaturedependent broadening of the emission spectrum [30].Being identified as torsional and librational degrees of freedom of conjugated rings at a molecular backbone [8,[31][32][33] as well as longitudinal acoustic modes of the backbone [5,31,34], low-energy vibrations constitute one of the first relaxation processes that take place within approximately 10 fs after photoexcitation to excited electronic states, known as structural (or geometric) relaxation [33,35,36].This relaxation, akin to high-energy molecular vibrations, is driven by exciton-vibration coupling [30,31,33].Considering that any subsequent relaxation or deactivation of electronic excitations typically follows this geometric relaxation [36] and given the short lifetime of polaritons in molecular systems, we propose that low-energy vibrations are one of the primary candidates for polariton thermalization.
In this work, we developed a microscopic theory for polariton thermalization via low-energy molecular vibrations coupled to electronic degrees of freedom and derived a simple analytical expression based on the Stokes shift and linewidth of the 0-0 vibronic peak in emission spectra.We calculated the polariton thermalization rate in a practical microcavity structure and provided a recipe for cavity design and the choice of molecular system to achieve the desired thermalization rate.Last but not least, we revealed important temperature dependence that provides a quantitative understanding of the observed thermalization bottleneck effect at low temperature [37].

II. HAMILTONIAN OF A MOLECULAR SYSTEM WITH STRONG LIGHT-MATTER AND VIBRONIC INTERACTIONS
In this Section, we develop a microscopic model to describe molecules that exhibit strong exciton-vibration (vibronic) interactions and are coupled to an optical cavity, e.g. the class of organic microcavities.We start our description with electronic and vibrational degrees of freedom of a molecular system itself, such as thin molecular films.
We consider our molecular system consisting of N mol molecules.Each of them hosts Frenkel exciton and N vib vibrational modes.Strong localization of Frenkel excitons [38], enable us to define total Hamiltonian of the organic molecular film Ĥmol as a sum of the individual terms Ĥ(m) The Hamiltonian of mth molecule is [6,7,11,39] where ω ) is the annihilation (creation) operator of the jth vibrational mode, Ω (m) is the interaction constant between the molecule and the incident light with the frequency ω Ω .The parameters for each molecule is slightly different due to inherent disorder of organic systems [40][41][42].
Given the strong exciton-vibration interaction in organic molecules that may exceed vibrational eigenfrequencies [6,11,31,40,43], we transition to the dressed exciton and vibrational operators.
where operators Ŝ(m) and B(m) j are the annihilation operators of the dressed excitons and dressed molecular vibrations.This operators fulfill the following commutation relations where we introduce the displacement operator and the energy of the vibrationally dressed excitons In the next step, we incorporate the modes of the electromagnetic field inside the cavity into the system.We also replace the energy of the dressed excitons, ω (m) 0 , dressed vibrations, ω (m) vj , and their interaction constant, Λ (m) j , of each individuall molecule by the mean values, ω 0 , ω vj , and Λ j .
The cavity modes possess different transverse (inplane) momenta, represented by ℏk ∥ (hereinafter referred as ℏk).Consequently, the full Hamiltonian of the systems now reads: where â † cavk (â cavk ) is the creation (annihilation) operator of a photon in the cavity with the wave vector k and frequency ω cavk , which obey bosonic commutation relation âcavk , â † cavk ′ = δ k,k ′ .Here, we assume that electric field of the kth mode is distributed in a plane parallel to the mirrors according to e ikr [12,26,37,[44][45][46][47].Vector r (m) points to the position of the mth molecule.The single molecule light-matter interaction energy is [48], where d (m) is the transition dipole moment of the molecule, and E (m) k represents electric field amplitude for "one photon" in the cavity with the in-plane momentum ℏk at the mth molecule position.
Expanding upon operators for the dressed excitons given in Eq. (3) and the vibrational displacement in Eq. ( 6) using the approximation e −Λj ( B(m) † ) we bring the total Hamiltonian of the system, as shown in Eq. ( 9) to the following form We restrict ourselves to the negative exciton-photon detuning; the case of positive detuning can be considered similarly.
It is convenient to introduce the collective operators effectively describing phase coherent excitonic and vibrational states and the Rabi frequency The operators (11) obey bosonic commutation relations ĉlk , ĉ † In the limit of large N mol and small enough excitation density, the operators for excitons (10) exhibit bosonic properties, i.e. ĉexck , ĉ † . This approximation is consistent with the most experiments [12,26] where the number of molecules in the illuminated area is ∼ 10 8 and the occupation of the excitons does not surpass ∼ 0.1.The Hamiltonian of the system Eq.( 9), in terms of the collective operators, reads: Considering the strong light-matter interaction with the collective vibrationally dressed excitonic states, we introduce operators for both the lower exciton-polaritons, denoted as ŝlowk , and the upper exciton-polariton, denoted as ŝupk , These polariton operators describe delocalized lightmatter states characterized by the in-plane momentum ℏk.Note that, despite disorder in molecular systems, the value k is a good quantum number, at least in strongly coupled organic microcavities.Indeed, the bottleneck effect in polariton thermalization, observed experimentally at low temperatures [37], underscores the crucial role of the density of states change in the vicinity of the inflection point in the momentum space [50].Angle-resolved experiments on non-ground state condensation further demonstrate the in-plane momentum conservation [51].The importance of a proper quantization of states in the momentum space is also evident in recent single-photon stimulated experiments [12], where the momentum distribution of a seed beam has to be matched to the size of the system in order to provide nonlinearity at the single photon level.Now, we can transform our Hamiltonian in Eq. ( 13) into the basis of the exciton-polariton states where we omit the non-resonant terms.We define the polariton eigenfrequencies as follows: The ratio between the excitonic (material) and photonic parts of polariton states are defined by the Hopfield coefficients sin φ k and cos φ k [29], where φ k is We note that the lower polaritons inherit from cavity photons the quadratic dispersion in the vicinity of the ground state, ω lowk ≈ ω low0 + α pol k 2 , where As evident from the Hamiltonian in its final form given by Eq. ( 16), polaritons are the hybrid light-matter states that include matter components featuring both electronic and vibrational degrees of freedom.In the next section, we show that the nonlinear interplay between them gives rise to polariton thermalization through low-energy molecular vibrations.

III. THERMALIZATION RATE OF POLARITONS
Due to the complexity of the vibrational landscape with its various degrees of freedom we treat low-frequency vibrations through introducing a reservoir and exclude them using Born-Markov approximation.To justify this assumption, we compare the degrees of freedom of the system with those of the reservoir.The Born-Markov approximation is applicable for the reservoir having significantly more degrees of freedom than the system itself [52,53].This scenario is precisely what occurs in our system.Indeed, the number of states for lower polaritons roughly matches the number of states in the corresponding cavity.For instance, a cavity with its fundamental mode characterized by a wavelength λ has a number of states N states ≈ πk 2 max S/(2π) 2 , where S is the area of interest [17].In room temperature BEC experiments [12,26,37] these parameters are typically the following: λ ≈ 500 nm and S ≈ 500 µm 2 , which result in N states ≈ 10 4 number of states.Furthermore, the corresponding number of states within the energy range k B T at the ground polariton state is about 10 3 .The number of molecules / molecular chains is good estimate for the total number of molecular vibrations, which is around 10 8 within the region of interest.In addition, we would like to highlight some empirical evidences in support to the Born-Markov approximation.For example, the thermal distribution of polaritons above the BEC threshold [37] with the temperature close to the environment is a signature of at least partial thermal equilibrium.This observation aligns with the expectations of the Born-Markov approximation.
The Hamiltonian (16) allows us to derive the contribution of the molecular vibrations to the thermalization rate of polaritons.We separate all the vibrational modes of the molecules into two groups: high-frequency vibrations and low-frequency vibrations.We order all the vibrational modes such that ω vn < ω vn+1 and denote M as the number of the vibrational mode for which ω vM < Γ and ω vM +1 > Γ, where Γ is the standard deviation of dressed excitons transition frequencies.All the vibrational modes with the natural frequency below and equal ω vM we call low-frequency vibrations, the rest vibrational modes we call high-frequency vibrations.We take into account the low-frequency vibrations effectively, as a reservoir, while we consider high-frequency vibrations explicitly.Thus, from Hamiltonian (16), we obtain the thermalization rate between two arbitrary lower polariton states with the wave vectors k and k ′ (ω lowk > ω lowk ′ ) where we use the continuous limit of the distribution of the low-frequency vibrations.ν(ω) is the density of the states of the molecular vibrations at frequency ω, Λ(ω) is the square of Huang-Rhys factor at the frequency ω, The frequency difference ∆ω kk ′ is limited from below due to the finite size of the system.In general, the finite size of the system leads to the discrete spectrum [54].For example, in 2D organic microcavities, the minimal value of ∆ω kk ′ is reversely proportional to the area of the system [17], S, and equal as long as we can approximate the real dispersion of the lower polaritons with the quadratic one [20].We note, that this procedure is independent of the correlation length of the polaritons.Indeed, the correlation radius is determined by the density matrix of the polaritons [20,55].Moreover, even for a very large system, the correlation length can be small compared to the size of the system below the condensation threshold and infinite above the condensation [20,55,56].The thermalization rates ( 21)-( 22) obey the Kubo-Martin-Schwinger relation [57] The corresponding Lindblad operator [12,58] for the density matrix of polaritons ρ describes the energy flow from the lower polaritons having wave vector k 2 towards ones with wave vector k 1 .Equations ( 21) and (22) show that the thermalization rates, denoted as γ therm ∼ Λ 2 = S H , are directly proportional to the Huang-Rhys parameter (S H ).This parameter determines the strength of the exciton-vibration interaction and establish connection to spectroscopic properties of bare molecules.Interestingly, a Hamiltonian of the same type is also used to describe Raman scattering [39,59], as well as the transition process from bright excitons to lower polaritons, which occurs through highfrequency vibrations [12,13].Indeed, the cross-section for Raman scattering is proportional to S H . Therefore, polariton thermalization via low-frequency vibrations predominantly involves Raman-allowed transitions.Whether these vibrations are also IR-allowed is not particularly important for this mechanism.
Given the ultrafast timescale of geometric relaxation in the electronically excited states of highly conjugated molecular systems [33,35,36], the low-frequency vibrations emerge as the primary candidates for polariton thermalization in organic cavities.This thermalization process conserves the total number of polaritons, but does not preserve the total energy [19,20].As a polariton transitions to a state with lower in-plane momentum, it loses some energy.This released energy is resonantly absorbed by molecules through the excitation of low-frequency molecular vibrations exhibiting the highest exciton-vibration coupling, which, in turn, dictates the net polariton thermalization rate.
The thermalization rates ( 21) and ( 22) depend on the properties of the molecular vibrations with particular frequency, namely Λ 2 (ω) and ν(ω), the angles φ k and φ k ′ (see Eq. ( 19)), and the frequency difference between them ∆ω kk ′ = |ω lowk − ω lowk ′ |.To estimate the thermalization rates (21) and (22) we should consider the properties of both molecular vibrations and polaritons.For polariton states with wave vectors k and k ′ for which ∆ω kk ′ , ∆ω k0 ≪ |ω exc − ω low0 |, we can expand sin 2 (φ k ′ − φ k ) to the following approximate equation Although local properties of low-frequency molecular vibrations are hard to access experimentally [60,61], we can reliably estimate ωΛ 2 (ω)ν(ω) and ω 2 Λ 2 (ω)ν(ω) in the low-frequency range, by their average values A 1 /ω M and A 2 /ω M , respectively, where It is very useful to consider Eq. ( 21)-( 22) in high and low temperature limits separately.Inspired by the recent experimental progress in room temperature BEC and related phenomena we focus on high temperature limit as the most practical.
In the high temperature limit k B T ≫ ℏω M we have ℏ∆ω kk ′ ≪ k B T and therefore can approximate n v (∆ω kk ′ ) ≈ k B T /(ℏ∆ω kk ′ ) obtaining the following thermalization rates: for any relation between ω lowk and ω lowk ′ .Then we estimate the average thermalization rate over frequencies ∆ω kk ′ from 0 to ω M using Eq. ( 29) For 2D organic microcavities, the minimal value of ∆ω kk ′ is reversely proportional to the area of the system, S, and equal ∆ω min = 4πα pol /S.The major contribution to the average thermalization rate γ kk ′ therm and γ k ′ k therm according to Eq. ( 32) corresponds to the transitions between the nearest polariton states.
Thus, we obtain thermalization rate at the high temperature limit In the low-temperature limit, where k B T ≲ ℏω M , we cannot neglect unity in the factor (1+n v (∆ω kk ′ )).Therefore, we must proceed with Eq. ( 21)- (22) in its complete form.Here we use average A 2 /ω M to estimate ω 2 Λ 2 (ω)ν(ω) and obtain the rates at the low temperature limit.
where we set ω lowk > ω lowk ′ .In the next Section, we show how the features of the emission and absorption spectra of an organic film is linked to the integral properties of low-frequency molecular vibrations ( 29) and (30).

IV. THE ROLE OF LOW-ENERGY VIBRATIONS IN ABSORPTION AND EMISSION SPECTRA
Determining the individual parameters for each of the low-frequency vibrational modes in densely-packed molecular systems even without optical cavities stands as an extremely challenging experimental problem.Luckily, ensemble averaged properties of vibrational modes and their coupling to electronic states can be obtained from emission and absorption spectra, for example by analysis of the widths and relative spectral positions of the maxima as function of temperature.In this section we provide theoretical analysis for the spectral properties of bare molecular systems (such as thin organic films) without a cavity and establish connection to A 1 and A 2 parameters that define the polariton thermalization rate.
In general, the theoretical calculation of the emission and absorption spectra of the molecular system requires the knowledge of Hermitian evolution of the organic film and the relaxation processes.The Hermitian evolution is given by the Hamiltonian (1), while the relaxation processes can be described within Lindblad approach by introducing the density matrix ρ of excitons and vibrations hosted by the molecules in the organic film.The nonequilibrium dynamics of ρ is governed by the Lindblad master equation [48,62] where L Ân (ρ) is Lindblad operator, Ân is the relaxation operator We consider energy dissipation of both: electronic and vibrational subsystems of molecules.The incoherent pumping of dressed excitonic states is also introduced via Lindblad operator.Therefore, all relaxation processes and the incoherent pumping correspond to the following operators γ vj .Considering that the relative differences in the parameters, stemming from disordered nature of molecular systems, are all of comparable magnitude, the most significant absolute difference lies in ω (m) 0 .Therefore, we can state that the molecules within the system are distinguishable solely by the energy of the dressed excitons.Given this, we assume that the dephasing of excitons in the system originates primarily from two factors: (i) inhomogeneous broadening due to the difference in the energies of the dressed excitons, and (ii) the impact of lowenergy vibrations on the dynamics of the excitons.These factors are explicitly considered in Hamiltonian (1).And thus we do not need to account for them separately in the Lindblad operators given by Eq. (37).
In light of the abovementioned, we will disregard the variations in all other molecular parameters, therefore, the index (m) is omitted hereinafter as identified in Eq. (39).

A. Emission and absorption spectra
The stationary emission spectrum of the molecular system itself can be calculated from two-time correlator σ † (t + τ )σ(t) [62].The presence of the weak incoherent pumping implies the assumption: γ pump ̸ = 0, γ pump ≪ γ diss , Ω = 0 in Eq. ( 37), which corresponds the low probability of any of the molecules to be in the excited state [62].
The absorption spectrum is determined by the imaginary part of the permittivity that characterizes linear response on the weak external field [63].Therefore, absorption spectrum is proportional to −⟨σ(t)e iωΩt ⟩/Ω at the condition γ pump = 0, Ω ̸ = 0.
Given the typical molecular density of organic thin films utilized in polariton research is 10 20 −10 21 cm −3 we deal with the continuous distribution of exciton energies.Moreover, taking into account the disorder mentioned above it is reasonable to assume normally distributed energies with the standard deviation Γ.We also assume that Γ ≫ γ diss /2 + γ vΣ .Thus, the inhomogeneous broadening leads to the Gaussian lineshape [2,[64][65][66][67][68] of the spectral peaks.
Using the master equation ( 37) we find emission and absorption spectra of the bare molecular system without a cavity (see Supplementary Information)  40) and ( 41), the spectra (c) and (d) take into account low-frequency vibrational modes effectively according to Eq. ( 43) and ( 44), the spectra (e) and (f) utilize the exact model while neglecting low-frequency modes.
where the sum for k 1 , k 2 , ..., k M runs over all nonnegative integers and The expression for the absorption spectrum (Eq.( 41)) is distinguished from that of the luminescence spectrum (Eq.( 40)) by the substitution of the sign preceding N vib j=1 ω vj k j .Consequently, the non-linear interaction between electronic states and low-frequency vibrations of molecular nuclei is a contributing factor to the Stokes shift [31].
We can effectively take into account low-frequency vibrational modes and obtain from Eq. ( 40) and ( 41) the reduced expression for emission and absorption spectra (see Supplementary Information) where we assume n vj ≈ 0 for high-frequency vibrations and we introduce In the limit of the continuous distribution of the lowfrequency vibrations Eq. ( 45) and ( 46) become here we make use of the same notations as in Eq. ( 21)-Eq.( 22), in particular, n v (ω) is determined by Eq. ( 23).

B. Stokes shift
From Eq. ( 43)- (46), it is evident that, at a given temperature, the line shape of the emission spectrum does not enable us to differentiate between inhomogeneous broadening and the effects of low-frequency vibrations.However, a combined analysis of both the emission and absorption spectra allows for this distinction.Unlike low-frequency vibrations, inhomogeneous broadening does not cause the Stokes shift between 0-0 vibronic peaks in the absorption ω abs and emission ω em .Therefore, the Stokes shift, ω abs − ω em , provides insights into the low-frequency vibrational dynamics.
Indeed, from Eq. ( 45) and ( 47) one can obtain the Stokes shift: The Stokes shift is nearly independent of temperature.Therefore, analysis of the emission and absorption spectra of the molecular system itself (without a cavity) enables us to infer the net properties of low-frequency vibrations -essential for the polariton thermalization rate as delineated in Eq. (29).

C. Temperature-dependent broadening
The second integral property of the low-frequency vibrations can be determined from the temperature dependence of the linewidth of the 0 − 0 vibronic peak in the emission spectrum.At high temperatures, k B T ≳ ℏω M , expression for the linewidth Eq. ( 48) takes the following form: The obtained equation above not only outlines the temperature dependence of the linewidth for the 0-0 vibronic peak with respect to the Stokes shift but also emphasizes the self-consistency of the theory presented.Furthermore, Eq. ( 50) provides the means to calculate the degree of inhomogeneous broadening denoted by Γ 2 .At the low-temperature limit, the linewidth of the 0−0 emission peak (47) approaches a constant value.
Given the value of Γ 2 is extracted from the analysis at high temperatures limit discussed above, using Eq. ( 51) we can now obtain the second net property of lowfrequency vibrations (30) that define polariton thermalization rate.

V. CALCULATION OF THE THERMALIZATION RATE IN AN ORGANIC MICROCAVITY
Here we focused on organic Fabry-Pérot microcavities as the prototype system for these calculations; however, the findings are applicable to other cavity types, including plasmonic [15,16,22,69], photonic crystal [70], microdisk/microring [71], defect and gap-based cavities [72][73][74][75][76][77], among others.The decision to choose a particular material system is motivated by two main reasons: (i) the system must exhibit strong light-matter interaction and exciton-vibration coupling, and even more importantly, it should be capable of supporting polariton Bose-Einstein condensation (BEC); (ii) the chosen molecular system itself should be extensively studied to ensure that our results can be benchmarked against established findings.Organic microcavities based on methyl-substituted conjugated ladder-type polymer (MeLPPP) emerge as the ideal candidates for this analysis.In addition, the temperature-dependent polariton thermalization investigated experimentally [37] offers invaluable insights for this study shedding light onto the role of low-energy molecular vibrations.
To extract polariton thermalization rate we proceed with the following steps.First, we calculate the emission and absorption spectra of the bare MeLPPP layer at different temperatures and extract the net properties FIG. 3. Polariton thermalization rate calculated from the spectral analysis according to Eq. ( 35) as the function of lightmatter interaction strength (Rabi frequency ΩR) and the energy of the ground polariton state (ℏω lowk=0 ).
of the low-frequency vibrations according to Eq. ( 29)- (30).Second, we use the extracted net values of A 1 and A 2 to estimate polariton thermalization rate Eq. ( 35) and (36) as the function of light-matter interaction (12) and the energy of the ground state of polaritons.
Our comparison of the emission and absorption spectra starts by employing three separate approaches: (i) the exact microscopic model as defined by Eq. ( 40)-( 41); (ii) the reduced model accounting for the net effect of low-energy vibrational modes Eq. ( 43)- (44); and (iii) the exact model again, but this time excluding all low-energy modes.Figure 1 presents the results of these calculations at two different temperatures, 6 K and 300 K. Our reduced method agrees well on a quantitative level with the exact microscopic model that includes all mentioned vibrational modes.Furthermore, it is clear that neglecting low-energy modes fails to accurately reproduce the Stokes shift and linewidth in the spectra.
The comparison between our theoretical predictions and the experimental absorption and emission spectra reported in the literature demonstrates a quantitative match, including the observed slight asymmetry in the spectral lines [42,61,81].It is important to note that our theory primarily accounts for the effects of low-energy vibrational (geometric) relaxation and does not take into account factors such as intermolecular interactions or spin multiplicity change.As a result, this theoretical approach is most effective for molecular systems that ex- hibit a pronounced mirroring of vibronic replicas in their spectra.
The Stokes shift and the linewidth of the 0-0 emission peak can be extracted from Figure 1(a,b).We have expanded this dataset with additional calculations using the exact model across a range of temperatures from 6 K to 400 K. Figure (2) shows the Stokes shift and the linewidth as the function of temperature.Importantly the results obtained from the reduced model using analytical expressions Eq. ( 45)-( 46) are in good agreement with the exact method (Fig. 2) and experimental data [42].Although the Stokes shift remains constant with temperature changes (Fig. 2a), the linewidth stays constant only at low temperatures and increases linearly with temperature, as illustrated in Fig. 2b.Finally, the joint analysis of the emission and absorption spectra at different temperatures allows us to determine the net values A 1 ≈ 18 meV and A 2 ≈ 200 meV 2 of the lowfrequency vibrations, defining polariton thermalization rates Eq. ( 35)- (36).
Next, we demonstrate how the net values A 1 and A 2 , obtained from the spectral analysis of the bare MeLPPP layer, can be applied to calculate the polariton thermalization rate in practical microcavities [12,26,37].We make use of the equation (36) to extract thermalization rates between neighboring states in the momentum space, for example ℏk ′ = 0, and ℏk = ℏk ′ + ℏδk, with the corresponding energy being ℏω lowk = ℏω lowk ′ + 4πα pol /S, where α pol is defined by Eq. (20).As the cutoff frequency for low-energy vibrations we use ω M ≈ 200 cm −1 .
Figure 3 presents the extracted values of the polariton thermalization rate as the function of light-matter interaction (Rabi frequency Ω R ) and the energy of the ground polariton state (ω lowk=0 ).It should be noted that the polariton thermalization rate has not been directly measured to date and has always been treated as a variable parameter in the fitting of experimental data by microscopic and/or mean-field models [12,22,[82][83][84].Here, we introduce an independent method that provides direct access to the thermalization rate in molecular polariton systems without relying on any adjustable parameters.
Even though thermalization processes reflect the matter behaviour in polariton states, the design of the cavity is expected to have a notable impact, especially regarding aspects tied to the Hopfield coefficients of light-matter states, such as exciton-cavity detuning and mode volume.The general rule of a thumb is that an increase in the material component of the polaritons results in a higher thermalization rate (Fig. 3).As depicted in Figure 3 the polariton thermalization rate monotonically increases with the Rabi frequency.Despite the explicit dependence on the number of molecules in Eq. ( 21)- (22) and Eq. ( 35)-( 36), the thermalization rate does not decrease with the increase in N mol , because Ω R ∝ √ N mol .Moreover, the thermalization rate grows with the concentration of the molecules due to the increase in material component of the polaritons, which is described by the Hopfield coefficients in Eq. ( 21)-( 22).This observation is in agreement with polariton plasmon systems [22].
Temperature is another critical parameter in polariton thermalization.The thermalization rate given by Eq. ( 35)-( 36) exhibits very strong dependence on temperature due to the factor n v (∆ω kk ′ ). Figure 4 demonstrates the peculiar thermalization behavior.First, the temperature rise leads to an almost linear increase in the thermalization rate to the nearest neighbor polariton state.Secondly, at higher temperatures, the range of polariton states capable of effective thermalization to a specific state expands, or alternatively, it increases the thermalization length, defined as the number of states to which a polariton can efficiently thermalize.These results explain the bottleneck effect in polariton thermalization observed in organic microcavities at low temperatures [37].It underlines the distinct role of thermalization processes in nonequilibrium Bose-Einstein condensation of molecular exciton-polaritons.
Our theory enables direct access to the polariton thermalization rate from standard spectroscopic measurements of bare molecules.Here we provide joint theoretical analysis for the ground and excited state polariton condensation recently achieved in experiments [51].We use numerical model developed in Ref. [12] with the value of thermalization rate obtained here and experimental FIG. 5.
Numerical simulation of the ground (left) and excited state (right) polariton condensation.Here we use microscopic model Ref. [12] with experimental parameters from Ref. [51] and effective thermalization rate of γ therm = 5 • 10 −10 eV extracted from the analytical theory.parameters from Ref. [51] to simulate E, k−distributions of polariton occupation.Figure 5 shows results of numerical simulations for the ground state (left) and excited state k = 2.5 µm −1 (right) polariton condensation.It is worth mentioning that the condensation in excited states with high in-plane momenta is quite sensitive to the thermalization rate.Stable condensation in this regime requires a delicate balance between the gain stimulated through the seed beam and the losses mainly coming from finite cavity photon lifetime and thermalization processes favoring downward relaxation to the ground state.Therefore excited state condensation appears to be an excellent testbed to benchmark thermalization theories.The effective polariton thermalization rate extracted from our analytical theory γ therm = 5 • 10 −10 eV demonstrates quantitative agreement with the experiment [51] throughout the simulations.Note, the model does not use any free parameters, see Supplementary Information for further details.

VI. CONCLUSION
In this work, we investigate the microscopic origins of polariton thermalization in organic systems.We have developed a theoretical framework that encompasses both strong light-matter and exciton-vibration couplings.The main focus of our study is the role of low-energy vibrations in shaping the spectral properties of molecular systems and their interplay with polariton thermalization when strongly coupled to an optical cavity.Analytical expressions for emission and absorption spectra have been derived to effectively incorporate low-frequency vibrational coupling.We introduced net parameters accounting for the ensemble averaged properties of lowenergy modes, which can be extracted directly from the Stokes shift (A 1 ) and the temperature-dependent linewidth of the 0-0 vibronic peak in the emission spec-trum (A 2 ) of the molecular system without a cavity.In the next step, we established the correspondence between the spectral properties of the molecular system and the polariton thermalization within the cavity.
We have derived a simple analytical expression for the polariton thermalization rate that is proportional to A 1 and A 2 at high temperatures, meaning a high thermalization rate requires a large Stokes shift and offset linear dependence on the spectral broadening with temperature.Moreover, the polariton thermalization rate depends strongly on temperature and the cavity properties.Finally, we applied the developed formalism to a practical microcavity structure based on the MeLPPP conjugated polymer and calculated the thermalization rate without use of any free parameters.Our results showcase remarkable agreement with recent experimental reports of nonequilibrium polariton condensation in the ground [12,26,37] and excited [51] states, and explain the thermalization bottleneck effect at low temperatures [37].
Our research lays the groundwork for understanding nonequilibrium polariton condensation in molecular systems, its interplay with vibrational degrees of freedom, and serves as a guideline for future experimental studies, providing recipes to choose the proper material system and cavity design, as well as conditions to control the thermalization behavior in strongly coupled molecular systems.
To access spectral properties we apply quantum regression theorem [12] to generate equations for two-time correlator functions in a form: ⟨ŝ † k (t+τ )ŝ k (t)⟩, where amplitudes of polariton states within region of interest defined as ŝk = ⟨ŝ k ⟩ (mean-field approximation).On the next step, we solve them numerically and apply Fourier transformation to the two-time correlators at each mode followed by integration over time.
We numerically solve the differential equations for N = 31 modes at the lower polariton branch including the ground state.The parameters for the model adopted from experimental data Ref.[51], the main ones are the following: • Light-matter interaction: ω exc = 2.72 eV, ω cav = 2.64 eV, Ω R = 85 meV, α cav = 2.2 meV • µm 2 (C3) where we set A 1 ≈ 18 meV, A 2 ≈ 200 meV 2 and ω M = 200 cm −1 (see main text Section V).The numerical simulations ignore some of the details in the thermalization process, accounting for them only effectively.For example, here we assume that thermalization rate is constant within the region of interest, i.e. it does not depend on the in-plane momentum ℏk: γ k1→k2 therm = γ therm when ω 1 > ω 2 and is equal to γ therm exp (−ℏ(ω 2 − ω 1 )/k B T ) otherwise.Therefore, we cannot use our estimation Eq. (C3) directly.Nevertheless, by averaging over many thermalization steps of Eq. (C3), our analytical theory can provide the effective thermalization rate per step of γ therm ≃ 5 • 10 −10 eV that we use in the numerical simulations.
of the exciton of the mth molecule, ω (m) vj is the eigenfrequency of the jth vibrational mode of the mth molecule, Λ (m) j is the square of Huang-Rhys factor of the jth vibrational mode of the mth molecule, σ(m) (σ (m) † ) is the annihilation (creation) operator of the exciton, b(m) transition to the dressed operators (3)-(4) diagonalizes the part of Hamiltonian (2) corresponding to excitonvibration interaction and bring the Hamiltonian to the following form Ĥ(m) mol = ℏω (m) 0 energy k B T can exceed low-energy molecular vibrations ℏω (m)

FIG. 1 .
FIG. 1.The emission (blue solid line) and absorption (red dashed line) spectra of a bare MeLPPP layer at different temperatures (a), (c), (e) T = 6 K; (b), (d), (f) T = 300 K.The spectra in (a) and (b) explicitly take into account all the vibrational modes according to exact expressions Eq. (40) and (41), the spectra (c) and (d) take into account low-frequency vibrational modes effectively according to Eq. (43) and (44), the spectra (e) and (f) utilize the exact model while neglecting low-frequency modes.

FIG. 2 .
FIG. 2. Temperature dependence of the Stokes shift -(a) and linewidth -(b) of 0 − 0 emission peak.The dashed line shows the asymptotic behaviour of the linewidth at high temperatures.Red double arrows mark the net values of the low-frequency vibration, namely A1 and A2 parameters.

FIG. 4 .
FIG. 4. Temperature dependence of the polariton thermalization rate between the ground polariton state (k ′ = 0) and the state with given wavevector k (hor.axis).Here we set ℏω cavk=0 = 2.64 eV and ℏΩR = 85 meV.For the sake of consistency all the other parameters are taken the same as for Fig 3.For these parameters and temperature T = 300 K Eq. (34) gives the estimation 3.4 • 10 −8 eV.