Molecule-photon interactions in phononic environments

Molecules constitute compact hybrid quantum optical systems that can interface photons, electronic degrees of freedom, localized mechanical vibrations and phonons. In particular, the strong vibronic interaction between electrons and nuclear motion in a molecule resembles the optomechanical radiation pressure Hamiltonian. While molecular vibrations are often in the ground state even at elevated temperatures, one still needs to get a handle on decoherence channels associated with phonons before an efficient quantum optical network based on opto-vibrational interactions in solid-state molecular systems could be realized. As a step towards a better understanding of decoherence in phononic environments, we take here an open quantum system approach to the non-equilibrium dynamics of guest molecules embedded in a crystal, identifying regimes of Markovian versus non-Markovian vibrational relaxation. A stochastic treatment based on quantum Langevin equations predicts collective vibron-vibron dynamics that resembles processes of sub- and superradiance for radiative transitions. This in turn leads to the possibility of decoupling intramolecular vibrations from the phononic bath, allowing for enhanced coherence times of collective vibrations. For molecular polaritonics in strongly confined geometries, we also show that the imprint of opto-vibrational couplings onto the emerging output field results in effective polariton cross-talk rates for finite bath occupancies.


I. INTRODUCTION
Molecules are natural quantum mechanical platforms where several atoms are interlinked via electronic bonds. The inherent coupling between the electronic transitions at optical frequencies and the mechanical nuclear motions (vibrons) at terahertz frequencies renders molecular systems ideal for the realization of quantum optomechanical effects. This is however different from the radiation pressure coupling mechanism in macroscopic systems, as optomechanical interactions in molecules intrinsically occur in a hybrid fashion involving a two-step process of photon-electron and vibron-electron (vibronic) interactions [1][2][3][4]. The vibronic coupling resembles the radiation pressure Hamiltonian (via a boson-spin replacement) which can be in the strong coupling regime since the strength of the coherent coupling can be comparable to the vibrational frequency. At cryogenic temperatures (e.g., at T ∼ 4 K), molecular vibrations are in their quantum ground state thus circumventing usual complications arising from additional optical cooling requirements [5]. Moreover, naturally occurring or engineered differences in the curvatures of the ground and excited state potential surfaces of the molecular electronic orbitals can lead to the direct generation of non-classical squeezed vibrational wavepackets [6]. These aspects suggest that molecular systems offer natural platforms, where one can exploit the inherent opto-vibrational coupling as a quantum resource.
When molecules couple to their condensed-matter environment, e.g. in the solid state, the mechanical modes of localized intramolecular vibrations (vibrons) are augmented by collective delocalized vibrational excitations of the host material (phonons), which allow for electron-phonon (polaron) couplings. In practice, coupling to a large number of phonon modes makes the study of molecular vibrations in the solid state notoriously challenging. Some of the challenges can be tamed under cryogenic conditions where experiments manage to reduce phonon coupling on the so-called zero-phonon line (ZPL) of the transition between |g, n ν = 0⟩ and |e, n ν = 0⟩ sufficiently to reach its natural linewidth limit. This can be verified in ensemble measurements, e.g. via hole burning, or in single-molecule spectroscopy [7]. A good example of an experimental platform is provided by dibenzoterrylene (DBT) molecules embedded in anthracene crystals [see Fig. 1(a)], exhibiting a lifetime-limited linewidth and near-unity radiative yields at cryogenic temperatures [8][9][10][11][12]. However, even if vibrational spectroscopy at the single-molecule level is readily accessible in the laboratory [13,14], a quantitative understanding of the couplings between the molecular vibrational modes and their internal and external degrees of freedom is still largely missing. In particular, a detailed study of decoherence sources is necessary. An open quantum system approach, such as employed in our treatment, can shed light onto a few aspects of coherent and incoherent vibrational dynamics and onto the light-matter interactions in the presence of vibrons, phonons and cavity-localized photon modes. Our formalism makes use of quantum Langevin equations which allows us to follow the evolution of system operators such as the electronic coherence and vibrational quadratures and to derive analytical results for the time dynamics of both expectation values and two-time correlations (needed for the computation of emission and absorption spectra). We find that closely spaced molecules can experience col-lective vibrational relaxation, an effect similar to the suband superradiance of quantum emitters in the electromagnetic vacuum. This can be exploited to decouple collective two-molecule vibrational states from the decohering phononic environment leading to the possibility of coherently mapping motion onto light and vice versa. In addition, at the level of the pure light-matter interface, coupling to confined optical cavity modes can increase the oscillator strength of the molecule by effectively reducing vibronic couplings [12].
Our formalism also allows us to treat problems relevant to experiments in cavity quantum electrodynamics with molecules, where standard concepts such as strong coupling or the Purcell effect can suffer important modifications once couplings between electronic transitions and vibrations are taken into account. To this end, we make use of analytical tools based on quantum Langevin equations [15] to account for an arbitrary number of vibron and phonon modes. Earlier theoretical works have either traced out the typically fast vibrational degrees of freedom [16,17], used limited numerical simulations, or focused mostly on aspects such as vibrational relaxation in solids [18][19][20][21], electron-phonon and electron-vibron couplings [22][23][24], temperature dependence of the zerophonon linewidth [25,26] and anharmonic effects [27,28].
However, it should also be borne in mind that the relevance of our treatment is not restricted to the physical system considered here as very similar effects also occur in related solid-state emitters such as quantum dots or vacancy centers in diamond. The coupling of such systems to photonic nanostructures has been studied quite extensively over the last years [29][30][31][32][33][34][35][36]. There is, furthermore, a general current interest in impurities interacting with a quantum many-body environment, such as molecular rotors immersed in liquid solvents [37,38], Rydberg impurities in quantum gases [39] or magnetic polarons in the Fermi-Hubbard model [40]. Our treatment can then be understood as a general model for the coupled dynamics of spin systems to many, possibly interconnected, bosonic degrees of freedom as illustrated in Fig. 1(c).

A. General considerations
We develop here a complex model where all interactions between light, electronic transitions, vibrons and phonons are taken into account for finite temperatures. We derive general expressions for the light scattered by a molecular system (of one or more molecules) embedded in a solidstate environment outside or inside an optical cavity [see Fig. 1(a)]. As schematically illustrated in Fig. 1(c) the light (mode a) couples to electronic transitions (Pauli operator σ) via a Tavis-Cummings Hamiltonian. These are in turn affected by the vibronic coupling to one or more molecular vibrations which leads to the red-shifted Stokes lines in emission [cf. Fig. 1(b)]. We focus here on a single mode with relative motion coordinate Q for the sake of simplicity. The solid-state matrix supports a multitude of bosonic phonon modes with displacements q k (k from 1 to N ) which directly modify the electronic transition leading to the occurrence of phonon wings in the emission and absorption spectra. In addition, molecular vibrons can deposit energy into phonons as a displacement-displacement interaction, leading to an irreversible process of vibrational relaxation. We will start with the description of the vibrational relaxation process in Sec. III since all subsequent effects will depend on this mechanism. We show that linear phonon-vibron couplings can already result in irreversible vibrational relaxation involving both single-and multi-phonon processes. Moreover, such dynamics can be either Markovian or non-Markovian, depending on the relation between the vibrational frequency and the maximum phonon frequency. For closely spaced molecules, the same formalism allows for the derivation of collective relaxation dynamics exhibiting effects similar to super/subradiance in dense quantum emitter systems. Classical light driving is included in Sec. IV by calculating absorption spectra for coherently driven molecules under the influence of vibronic and phononic couplings as well as thermal effects. We show that interestingly, the vibronic and electron-phonon couplings do not cause any dephasing dynamics even at high temperatures, i.e. the zero-phonon line is mainly lifetime-limited in the linear coupling model. Following a quantum Langevin equations approach, we derive absorption spectra for coherently driven molecules under the influence of vibronic and phononic couplings as well as thermal and finite-size effects. Finally, for molecular polaritonic systems in a cavity setting, we derive transmission functions of the cavity field (see Sec. V), showing the reduction of the vacuum Rabi splitting with increasing vibronic and phononic coupling, as well as phononic signatures in the Purcell regime. The effect of temperature on the asymmetry of cavity polaritons is quantified by deriving effective rate equations for the polariton cross-talk dynamics.

B. Hamiltonian formulation
We consider one molecule (later we extend to more than one) embedded in a bulk medium comprised of N unit cells. Our perturbational assumption is that, since the bulk is large, the guest molecule does not significantly change the overall modes of the bulk. The electronic degrees of freedom of the molecule are denoted by states |g⟩ and |e⟩ with the former at zero energy and the latter at ω 0 (we set ℏ to unity), corresponding to a lowering operator σ = |g⟩ ⟨e|. We assume only a pair of ground and excited potential landscapes with identical curvature along the nuclear coordinate and make the harmonic approximation, where the motion of the nuclei can be described by a harmonic vibration at frequency ν and bosonic operators b and b † , satisfying the usual bosonic commutation relations [b, b † ] = 1.
From the displacement between the minima of the two potential landscapes one obtains a vibronic coupling quantified by a dimensionless factor λ (the square root of the Huang-Rhys parameter) and described by a standard Holstein-Hamiltonian [41], where Q = (b + b † )/ √ 2 is the dimensionless position operator of the vibronic degree of freedom (the momentum quadrature is given by . The Holstein coupling also leads to a shift of the electronic excited state energy ω 0 +λ 2 ν, which is removed by the diagonalizing polaron transformation U el-vib . The polaron transformation U el-vib = e i √ 2λP σ † σ = |g⟩ ⟨g| + B † |e⟩ ⟨e| can be seen as a conditional displacement affecting only the excited state, where B † = e i √ 2λP is the inverse displacement operator for the molecular vibration creating a coherent state when applied to vacuum: B † |0 ν ⟩ = |−λ⟩. The Hamiltonian in Eq. (1) does not consider nonadiabatic vibronic coupling which would lead to off-diagonal coupling terms (proportional to σ x and σ y ) and which could drive electronic transitions. Such nonadiabatic terms become relevant if two potential surfaces come close to each other [42]. In Appendix H we briefly discuss how one could treat such terms in the Langevin equations of motion. One could also consider a difference in curvatures between ground (frequency ν) and excited state (frequencyν) potential surfaces which would result in a quadratic coupling term H quad el-vib = βQ 2 σ † σ with squeezing parameter of the vibrational wavepacket β = (ν 2 − ν 2 )/(2ν). We will assume that the vibron quickly thermalizes with the environment (via the fast mechanism of vibron-phonon coupling described below) at temperature T and achieves a steady state thermal occupancyn The electronic transition is coupled to the quantum electromagnetic vacuum which opens a radiative decay channel with collapse operator σ via spontaneous emission at rate γ. For a general collapse operator O with rate γ O we model the dissipative dynamics via a Lindblad term applied as a superoperator to the density operator ρ of the system. The vibronic coupling leads to the presence of Stokes lines in emission and to a mismatch between the molecular emission and absorption profiles. Following the stochastic quantum evolution of a polaron operatorσ = B † σ (vibrationally dressed Pauli operator for the electronic transition) analytical solutions for the absorption and emission spectra of the molecule can be derived in the presence of vibrons [15].
In addition to the coupling to internal vibrations of its nuclei, the electronic transition is also modified through coupling to the delocalized phonon modes of the crystal. We describe the bulk modes as a bath of independent harmonic oscillators with bosonic operators c k and c † k and frequencies ω k . The electron-phonon coupling (see Appendix A for derivations) can then be cast in the same Holstein form as for the vibron where the displacement operators refer to each individual collective phonon mode q k = (c k + c † k )/ √ 2 (the momentum operator is given by . The coupling factors λ k depend on the specifics of the molecule and the bulk crystal. Similarly to the vibronic case, the electron-phonon interaction can be diagonalized by means of a polaron transformation U el-phon = |g⟩ ⟨g| + D † |e⟩ ⟨e|, 2λ k p k is the product of all phonon mode displacements, signifying a collective transformation for all phonon modes. We will assume that the bulk is kept at a constant temperature and is always in thermal equilibrium with the individual mode thermal average occupancies amounting tō n k = [exp(ω k /(k B T ))−1] −1 . The coupling to the phonons gives rise to a multitude of sidebands in the absorption and emission spectra which coalesce into a phonon wing that becomes especially important at elevated temperatures. We will then follow the temporal dynamics of a collective polaron operatorσ = D † B † σ which includes both vibronic and electron-phonon couplings.
Phonons also affect the dynamics of the vibrational mode. Modifications of the bond length associated with the molecular vibration leads to a force on the surrounding crystal (and vice versa), giving rise to a displacementdisplacement coupling, The coupling coefficients α k are explicitly derived in Appendix A. In the limit of large bulk media, this Hamiltonian can lead to an effective irreversible dynamics, i.e. a vibrational relaxation effect. This is the Caldeira-Leggett model widely treated in the literature as it leads to a non-trivial master equation evolution which cannot be expressed in Lindblad form and is cumbersome to solve analytically [43][44][45]. To circumvent this difficulty, we follow the formalism of Langevin equations under the concrete conditions imposed by the one-dimensional situation considered here. We are then in a position to identify the Markovian versus non-Markovian regimes of vibrational relaxation conditioned on the phonon spectrum, namely on the maximum phonon frequency ω max of the system. We can additionally account for a finite phonon lifetime by including a decay rate γ ph k for each phonon mode. To perform spectroscopy, we add a laser drive modeled as H ℓ = iη ℓ σ † e −iω ℓ t − σe iω ℓ t with amplitude η ℓ . We will assume weak driving such that the assumption of thermal equilibrium is still valid. Furthermore, to treat various aspects of molecular polaritonics, we describe the dynamics of a hybrid light-matter platform by adding the coupling of a confined optical mode at frequency ω c to the electronic transition via a Jaynes-Cummings interaction The bosonic operator a satisfies the commutation relation [a, a † ] = 1 and the coupling is given by g = [d 2 eg ω c /(2ϵ 0 V )] 1/2 , where d eg is the electronic transition dipole moment and V is the quantization volume (ϵ 0 is vacuum permittivity). Spectroscopy of the cavitymolecule system can be done by adding a cavity pump H ℓ = iη c a † e −iω ℓ t − ae iω ℓ t with amplitude η c . The cavity loss is modeled as a Lindblad process with collapse operator a and rate κ. In standard cavity QED, depending on the magnitude of the coherent exchange rate g to the loss terms κ and γ one can progressively advance from a strong cooperativity Purcell regime to a strong coupling regime where polaritons emerge. We will mainly focus on analytical derivations of the effects of electron-vibron and electron-phonon couplings at finite temperatures on the emergence of a spectral splitting in the strong coupling regime as well as the transmission in the Purcell regime.

III. VIBRATIONAL RELAXATION
A decay path for the molecular vibration stems from its coupling to the bath of phonon modes supported by the bulk. While it is generally agreed that nonlinear vibron-phonon couplings contribute to the vibrational relaxation process, especially in the higher temperature regime [46], we restrict our treatment to a coupling in the bilinear form of Eq. (3). To understand the physical picture, we first show that in perturbation theory the bilinear Hamiltonian leads to a competition between fundamental processes that involve the decay of a vibrational quantum into superpositions of either single phonon states or many phonons adding together in energy to the initial vibrational state energy. Afterwards we proceed by writing a set of coupled deterministic equations of motion for the vibrational quadratures of the molecule {Q, P } and the collective normal modes of crystal vibrations {q k , p k }. This allows for the elimination of the phonon degrees of freedom and the derivation of an effective Brownian noise stochastic evolution model for the molecular vibrations. We illustrate regimes of Markovian and non-Markovian dynamics and show that an equivalent approach tailored to two molecules can lead to collective vibrational relaxation strongly dependent on the molecule-molecule separation.

A. Fundamental vibron-phonon processes
Let us consider an initial state containing a single vibrational quantum |1 ν , vac ph ⟩ that evolves according to the vibron-phonon bilinear Hamiltonian of Eq. (3). We aim to follow the fundamental processes leading to the energy of the vibration deposited in superpositions of single or multi-phonon states. We move to the interaction picture by removing the free energy with U = e iH0t with the free Hamiltonian H 0 = νb † b + N k=1 ω k c † k c k . The time-dependent interaction picture Hamiltonian, thus, becomes The formal solution of the Schrödinger equation can then be written as a Dyson series |ϕ(t)⟩ = T e −i t 0 dτH(τ ) |1 ν , vac ph ⟩. We can proceed by evaluating the first term in the series which leads to (see Appendix D for details) resonant scattering (ω k = ν) into singlephonon states |0 ν , 1 k ⟩ at perturbative rate α k t as well as off-resonant scattering (ω k ̸ = ν) into states |0 ν , 1 k ⟩ with probability inversely proportional to the detuning ω k − ν. We note that for ν > ω max , only off-resonant transitions are possible. The next order of perturbation theory, however, leads to multi-phonon processes where resonant transitions to states containing three phonons |0 ν , 1 j1 , 1 j2 , 1 j3 ⟩ become possible. The resonance condition reads ω j1 + ω j2 + ω j3 = ν for j 1 ̸ = j 2 ̸ = j 3 , and its amplitude is a sum over terms α j1 α j2 α j3 t/(ω j2 +ω j3 )(ω j3 −ν). These terms are small with respect to the rates of the resonances starting in the first order for ν ≤ ω max and 0 20 40 60  in total are comparable to the single-phonon scattering off-resonant terms.

B. Effective Brownian noise model
Formal elimination of the phonon modes (see Appendix B for details) leads to an effective Brownian motion equation for the momentum of the vibrational modė while the displacement follows the unmodified equatioṅ Q = νP . The effect of the phonon bath is twofold: (i) it can shift the vibrational frequency toν = ν − ν s and (ii) it leads to a generally non-Markovian decay kernel expressed as a convolution Γ * P = ∞ 0 dt ′ Γ(t − t ′ )P (t ′ ). For the particular case considered in the Appendix, the crystal-induced frequency shift is expressed as where k 0 denotes the spring constant of the host crystal, k M represents the spring constant of the vibron, and ∆k is a measure for the coupling of the molecule's relative motion to the bulk. For a discrete system, the expression for the damping kernel Γ(t) = k α 2 k ν/ω k cos(ω k t)Θ(t) involves a sum over all phonon modes which can be turned into the following expression in the continuum limit (N → ∞) Here, J n (x) denotes the n-th order Bessel function of the first kind, Θ(t) stands for the Heaviside function and Γ m = 2νν s /ω max is the decay rate in the Markovian limit. A similar expression is known from the Rubin model [47], where one considers the damping of a single mass defect in a 1D harmonic crystal. The zero-average Langevin noise term ξ is determined by the initial conditions of the phonon bath and can be expressed in discrete form as . We can treat Eq. (6) more easily in the Fourier space where the convolution becomes a product and the Fourier transform of the non-Markovian decay kernel Γ(ω) generally contains a real and imaginary part Γ(ω) = Γ r (ω) + iΓ i (ω). Figure 2(a) shows a plot of Γ r (ω) and Γ i (ω) where we can interpret the imaginary part as a frequency shift which is largest around ω = ±ω max . Together with the transformed equation for the position quadrature −iωQ(ω) = νP (ω) we then obtain an algebraic set of equations which allows us to calculate any kind of correlations for the molecular vibration, both in time and frequency domains. This will be needed later on for computing the optical response of the molecule.

C. Markovian versus non-Markovian regimes
The Markovian limit is achieved when the vibrational frequency lies well within the phonon spectrum ω max ≫ ν and Γ(ω) becomes flat in frequency space: Γ(ω) = Γ m . In this case the memory kernel tends to a δ-function: with the convention Θ(0) = 1/2. In the continuum limit, the correlations at different times are where the noise spectrum is expressed similarly to the case of a standard thermal spectrum for a harmonic oscillator in thermal equilibrium S th (ω) = Γ r (ω)ω/ν[coth (βω/2) + 1] in terms of the inverse temperature β = (k B T ) −1 . The  difference lies in the frequency-dependence of the real part of the decay rate function where the Heaviside functions provide a natural cutoff of the spectrum at ±ω max . While in the time domain, the noise is only δ-correlated at high temperatures and This property is helpful for analytical estimation of the molecular absorption and emission in the presence of non-Markovian vibrational relaxation. In frequency space, the response of the vibron to the input noise of the phonon bath is characterized by the susceptibility In Fig. 2(b), we plot |χ(ω)| 2 for the two cases ω max ≫ ν (Markovian limit) and ω max ≈ ν (non-Markovian regime) for identical Γ m . While in the Markovian regime, the susceptibility has two approximately Lorentzian sidebands with linewidth Γ m centered around ±ν, the finite frequency cutoff in the non-Markovian case leads to an unconventional lineshape with reduced linewidth and slight frequency shift. In the time domain, we can simulate the microscopic classical equations of motion for a large number of phonon modes and compare the results to the standard Markovian limit obtained from Brownian motion theory. This is illustrated in Figs which also includes the non-Markovian regime. At low temperatures β −1 ≪ ν the sideband at −ν is suppressed and the thermal spectrum can be approximated as S th (ω) = [2Γ r (ω)ω/ν]Θ(ω). This two-time correlation function of the momentum quadrature will be required later in the calculation of molecular spectra in section IV.

D. Collective vibrational effects
A collection of impurity molecules sitting close to each other within the same crystal will see the same phonon bath and can, therefore, undergo a collective vibrational relaxation process. This is similar to the phenomenon of subradiance/superradiance of quantum emitters commonly coupled to an electromagnetic environment, where the rate of photon emission from the whole system can be smaller or larger than that of an individual isolated emitter. In order to elucidate this aspect, we follow the approach sketched above, i.e. we eliminate the phonon modes to obtain a set of coupled Langevin equations for two molecules situated 2j sites apart from each other: The mutually induced (small) energy shift Ω = k α k,1 α k,2 /ω k and the mutual damping kernels Γ 12 = ν 2 k α k,1 α k,2 /ω k cos(ω k t)Θ(t) and Γ 21 = ν 1 k α k,1 α k,2 /ω k cos(ω k t)Θ(t) are strongly dependent on the intermolecular separation 2j (see Appendix C for full expression), whereas the individual decay terms Γ 1 and Γ 2 are given by the expressions derived previously. Importantly, now also the noise terms ξ 1 and ξ 2 are not independent of each other but correlated according to a separation-dependent expression specified in Appendix C. In the continuum limit N → ∞, the collective interaction kernels can be approximated with the aid of higher-order Bessel functions (assuming identical molecules ν 1 = ν 2 and consequently Γ 12 (t) = Γ 21 (t)): In Fig. 3(a), we plot the collective decay kernel as a function of time and intermolecular separations 2j. The collective effects do not occur instantaneously but in a highly time-delayed fashion [cf. Fig 3(b)]. We can interpret the collective interaction as an exchange of phonon wavepackets between the two molecules, where the wavepackets are traveling with the group velocity v g = ∂ω/∂q of the crystal (lattice constant a) at a maximum speed v max g ≈ aω max /2 (the high frequency components towards the band edge are slower). This leads to an approximate time of τ = 4jω −1 max for the wavepacket to propagate from one molecule to the other.
The collective interaction will also lead to a modification of the vibrational lifetimes of the molecules which we want to describe in the following. To this end, one can again proceed with a Fourier analysis of Eqs. (13). The expression for the non-Markovian collective interaction kernel in frequency space (between −ω max and +ω max ) reads where we introduced the Chebychev polynomials of first (T n ) and second kind (U n ). We are interested in the real part of the above expression which will give rise to a collectively-induced modification of the vibrational lifetime while the imaginary part corresponds again to a frequency shift. In Figs. 3(c) and (d) we plot the real and imaginary parts of Γ 12 (ω) for small distances j, respectively.
In the Markovian limit ω max ≫ ν everything becomes flat in frequency space and one can approximate Γ 12 (ω) = Γ 12 (0) = Γ m and consequently (Γ 12 * P i ) ≈ Γ m P i with i = {1, 2}. A diagonalization can be performed by moving into collective quadratures P + = P 1 +P 2 and P − = P 1 −P 2 (and identically for the positions) for which the equations of motion decouplė While one of the collective modes undergoes relaxation at an increased rate 2Γ m , the orthogonal collective mode can be eventually decoupled from the phononic environment. Of course, as the derivation we have performed is restricted to one-dimensional crystals, it would be interesting to explore this effect in three dimensional scenarios where both longitudinal and transverse phonon modes have to be considered with effects stemming from the molecular orientation as well as the influence of anharmonic potentials. A recent theoretical work also discusses phonon-bath mediated interactions between two molecular impurities immersed in nanodroplets with respect to the rotational degrees of freedom of the molecules [48].

IV. FUNDAMENTAL SPECTRAL FEATURES
Let us now consider a molecule driven by a coherent light field. We will make use of and extend the formalism used in Ref. [15] to compare the effect of Markovian versus non-Markovian vibrational relaxation, phonon imprint on spectra and temperature effects. To derive the absorption profile of a laser-driven molecule, one can compute the steady-state excited state population P e = ⟨σ † σ⟩ = η ℓ ⟨σ⟩ + ⟨σ⟩ * /(2γ). The average steadystate dipole moment can formally be written as (note that we are assuming weak driving conditions η ℓ ≪ γ such that the laser drive only probes the linear response of the dipole): The important quantity to be estimated is the correlation function for the displacement operators of the molecular vibration ⟨B(t)B † (t ′ )⟩, which is fully characterized by the Huang-Rhys factor λ 2 and the second-order momentum correlation functions: The stationary correlation ⟨P (t) 2 ⟩ = ⟨P 2 ⟩ = 1/2 +n includes the temperature of the environment and does not depend on the details of the decay process. The two-time correlations ⟨P (t)P (t ′ )⟩ (and consequently the vibrational linewidths of the resulting optical spectrum) are crucially determined by the details of the dissipation model derived in section III. In order to capture the non-Markovian character of the vibrational relaxation, we extend the method used in Ref. [15] by computing correlations in the Fourier domain and then transforming to the time domain.

A. The non-Markovian vibrational relaxation regime
Let us first consider the imprint of the particularities of the vibrational relaxation process onto the absorption and emission spectra when molecule-light interactions are taken into account. For the calculation of the momentum correlation function ⟨P (t)P (t ′ )⟩ one has to evaluate the integral in Eq. (12), where the susceptibility weighted with the thermal spectrum is given by the general expression non-Mark. with Γ i (ω) = Γ m ω/ω max between −ω max and +ω max . As discussed in the previous section, the real part of Γ(ω) determines the decay rate while the imaginary part leads to a frequency shift. Generally, performing the integral over the expression in Eq. (19) is difficult since the line shapes can be very far from simple Lorentzians. However, assuming a good oscillator (Γ m ≪ ν) and consequently a sharply peaked susceptibility that only picks frequencies around the vibrational resonance, we can obtain an effective modified frequency ν ′ and decay rate Γ ′ in the non-Markovian regime (however assuming ω max > ν) with ν ′ = ν 2 + Γ i (ν)ν 1/2 and Γ ′ = Γ r (ν ′ ). By expanding Eq. (19) around the poles of the denominator ω = ±ν ′ + δ and assuming |δ| ≪ |ν ′ |, one can then calculate the temperature-dependent momentum correlation function in the non-Markovian regime: with time delay τ = t − t ′ . This allows for an analytical evaluation of the integral in Eq. (17) (see Appendix G for detailed calculation) and leads to a steady-state excitedstate population of where we introduced L(n) = e −λ 2 (1+2n) λ 2n n! and B(n, l) = n l (n + 1) n−ln l . One can immediately obtain the result for the Markovian limit by replacing ν ′ → ν and Γ ′ → Γ m . Figures 4(a),(b) show a comparison between the momentum correlation function and the resulting steady-state population (normalized by the steady-state population of a resonantly driven two-level system P 0 = η 2 ℓ /γ 2 ) in the Markovian and non-Markovian regimes (for fixed Γ m ). We can see that non-Markovianity leads to modified spectral positions and linewidths of the vibronic sidebands while the ZPL is not affected by the dissipation process. The denominator of Eq. (21) contains a sum over Lorentzians with a series of blue-shifted lines with index n arising from the electron-vibration interaction which are weighted by a Poissonian distribution. Thermal occupation of the vibrational states can counteract this effect however by leading to red-shifted lines in absorption [see Figure 4(c)] with index l and weighted by a binomial distribution. However, as shown in Figure 4(d), for large vibrational relaxation rates Γ m ≫ γ the sidebands will be suppressed and absorption and emission of the molecule will mostly occur on the ZPL transitions |g, m ν ⟩ ↔ |e, m ν ⟩. While in the case of zero temperature the ZPL is solely determined by n = 0, for finite temperatures all terms with n = 2l can contribute to it.
An important quantity is the Franck-Condon factor f FC which measures the reduction of the ZPL intensity due to coupling to internal vibrations. This factor is given by the and does not depend on the vibrational relaxation of the molecule. Using the fact that (n + 1)/n = e βν , one can express Eq. (21) as a sum over just a single index in the limit 2λ 2 n(n + 1) ≪ 1 (see Appendix G for derivation) as with I n (x) denoting the modified Bessel functions of the first kind andN = n(n + 1). This expression is similar to the result known from the standard Huang-Rhys theory for emission and absorption [49], but it now additionally includes vibrational relaxations Γ ′ . The ZPL contribution (n = 0) at resonance is thus simply given by P e = η 2 ℓ f FC /γ 2 . The emission spectrum can be calculated from the Fourier transform of two-time correlations ⟨σ † (τ )σ(0)⟩. Considering the decay of an initially excited molecule, one finds that the emission spectrum is simply given as the mirror image (with respect to the ZPL) of the absorption spectrum which is why we restrict ourselves to the calculation of the absorption profile. and have neglected the effect of electron-phonon coupling. However, this can become a dominant mechanism at larger temperatures where all acoustic and optical phonon modes are thermally activated (> 50 K) and the probability of a ZPL transition is very small. To include electron-phonon coupling, the expression for the steadystate dipole moment [cf. Eq. (17)] has to also account for the displacement of the electronic excited state caused by the phonons Here, the coupling to phonons additionally leads to a renormalization of the electronic transition frequencỹ ω 0 = ω 0 − k λ 2 k ω k (polaron shift). The expression in Eq. (23) now jointly contains all of the effects: electron-phonon coupling, electron-vibron coupling and vibrational relaxation (through the correlation function ⟨B(t)B † (t ′ )⟩). Since we consider phonon modes to be independent of each other, the displacement correlation function of the phonons can be factorized When replacing the sum over k with an integral over ω in the continuum limit, this yields (neglecting phonon decay as it will not influence the spectra in the continuum limit): Here, we have introduced the spectral density of the electron-phonon coupling J(ω) = k |λ k ω k | 2 δ(ω − ω k ) = n(ω)λ(ω) 2 ω 2 where n(ω) denotes the density of states. In the one-dimensional derivation considered here we obtain for the spectral density The electron-phonon coupling constant λ 1D e-ph is derived in Appendix E and depends, among other things, on the displacement of the crystal atoms upon excitation of the molecule as well as on the spring constants between the molecule's atoms and the neighboring crystal atoms. Again, the cutoff at ω max arises naturally from the dispersion of the crystal. In the continuum limit considered here, this spectral density would lead to a divergence of the integral in Eq. (24) due to the high density of low-frequency phonons, which is a well known problem for 1D crystals [50,51]. This issue can be addressed by considering only a finite-sized 1D crystal with a minimum phonon frequency cutoff ω min > 0. However, one can instead also consider a spectral density stemming from a 3D density of states: where the electron-phonon coupling constant λ 3D e-ph now has units [s 2 ]. In Figs. 5(a) and (b) we plot the resulting absorption spectrum of the ZPL for 1D and 3D densities of states, whereby the exact shape of the phonon wing is determined by the spectral density function J(ω). While analytical expressions for the integral in Eq. (24) are difficult to obtain in the continuum case, we can express the absorption spectrum of the ZPL including phonon sideband in terms of discrete lines as Here the sum runs over all {n k } = n 1 , . . . , n N and {l k } = l 1 , . . . , l N . This can be seen as a generalization of the result in Eq. (21) for many modes where the N phonon modes are indexed by k and the function L k (n k ) accounts for the displacement of the excited state while the binomial distribution B k (n k , l k ) accounts for the thermal occupation of each mode. As one can see in Figs. 5(a) and (b), thermal occupation of the phonons leads to red-shifted phonon sidebands in absorption and eventually to a symmetric absorption spectrum around the zero-phonon line in the limit of large temperatures. Note that here we did not explicitly include phonon decay γ ph k as it does not influence the absorption spectra in the continuum limit (the phonon peaks overlap and are not resolved). However, one can easily account for a finite phonon lifetime by including it in the momen- Similarly to the Franck-Condon factor for vibrons, one defines the Debye-Waller factor f DW = ⟨D † ⟩ 2 = exp − ∞ 0 dωJ(ω)ω −2 coth(βω/2) which measures the reduction of the ZPL intensity due to the scattering of light into phonons. In Fig. 5(c) we show the behavior of the f DW in the 3D case for different coupling strengths at low temperatures, revealing a stronger temperature-dependence for larger couplings. The total reduction of the ZPL intensity as compared to the twolevel system case is then given by the product f FC · f DW .

C. Dephasing
Within the model we consider, where all interactions stem from a harmonic treatment of both intramolecular vibrations and crystal motion, the zero-phonon linewidth of the electronic transitions is largely independent of temperature. In reality, to account for higher temperature effects one needs contributions quadratic in the phononic displacements which has been theoretically pointed out and experimentally observed [52,53]. However, even in the linear regime, the fact that vibronic and electron-phonon couplings do not lead to significant dephasing is a non-trivial result. One could e.g. expect that the Holstein-Hamiltonian for electron-phonon coupling H H = ω 0 − k √ 2λ k ω k q k σ † σ + k ω k c † k c k which sees a stochastic shift of the excited electronic level should lead to a dephasing of the ground-excited coherence ⟨σ⟩. One reason is the similarity to the pure dephasing of a two-level transition subjected to a noisy laser undergoing evolution with the Hamiltonian [ω 0 +φ(t)]σ † σ where the frequency is continuously shaken by a white noise stochastic term of zero average and obeying ⟨φ(t)φ(t ′ )⟩ = γ deph δ(t − t ′ ). It is straightfoward to show that the time evolution of the coherence in this case becomes ⟨σ(t)⟩ = ⟨σ(0)⟩ e −iω0t e −γ deph t such that the correlations of the noise indicate the increase in the linewidth of the transition [54]. Similarly, one could expect that the zero-averaged quantum noise stemming from the shaking of the electronic transition in the Holstein-Hamiltonian would lead to the same kind of effect. However, computing the exact time evolution of the coherence in the interaction picture [with HamiltonianH H (t)] one obtains: where the time-ordered integral can be resolved by a second-order Magnus expansion, confirming the result already known from the polaron picture. The correlation (for a single mode ω k ) shows a cosine-term similar to the dephasing but which does not continuously increase in time. For small times t ≪ ω −1 k , the cosine-term can be expanded and the dephasing rate can be approximated by γ deph = λ 2 k (n k + 1/2)ω 2 k t while for larger times the rate goes to zero (the time scale is set by γ −1 ). In the continuum limit, the time-dependent dephasing rate γ deph (t) = −ℜ [φ(t)] expresses as In accordance with Figs. 5(a),(b) we can see that linear Holstein coupling can consequently lead to a temperaturedependent zero-phonon line if there is a large density of low frequency (long wavelength) phonon modes with ω k < γ which is the case in 1D but not in higher dimensions. This peculiarity of dephasing in 1D has already been discussed in the literature [50,51,55]. It is however also well established within the literature that the major contribution of the experimentally observed temperaturedependent broadening of the zero-phonon line is caused by a higher-order electron-phonon interaction of the form [52,[55][56][57] H quad el-phon = with the coupling constant of the quadratic interaction β kk ′ . This form of the interaction can stem either, within the harmonic assumption, from a difference in curvatures between the ground and excited state potential surfaces or from anharmonic potentials.

V. MOLECULAR POLARITONICS
It is currently of great interest to investigate the behavior of hybrid platforms containing organic molecules interacting with confined light modes such as provided by optical cavities [16,58,59] or plasmonic nanostructures [60,61]. Such light-dressed platforms have been studied both at the single-and few-molecule level [12,16,60] as well as in the mesoscopic, collective strong-coupling limit [62][63][64]. In these cases, the strong light-matter coupling leads to the formation of polaritonic hybrid states with both light and matter components. Experimental and theoretical works are currently exploring fascinating enhanced properties such as exciton and charge transport [65][66][67][68][69], superconductive behavior [70,71] and modified chemical reactivity [72][73][74][75][76]. There is also recent interest in the modification of nonadiabatic light-matter dynamics at so-called conical intersections leading to fast nonradiative decay of electronic excited states [42,77]. It has been recently shown that the Purcell regime of cavity QED can result in a strong modification of the branching ratio of a single molecule and suppress undesired Stokes lines [12]. Recent theoretical works account for the vibronic coupling of molecules by solving a Holstein-Tavis-Cummings Hamiltonian which leads to the occurence of polaron-polariton states, i.e. light-matter states where the hybridized states between the bare electronic transition and the light field additionally get dressed by the vibrations of the molecules [15,[78][79][80][81][82][83][84]. Many models rely on numerical simulations and are based on following the evolution of state vectors under simplified assumptions assuming only vibronic interactions and finite temperature effects. We employ here the approach of the last section and add a Jaynes-Cummings interaction of a molecule in the phononic environment with a localized cavity mode. A weak laser drive maps the intracavity molecular polaritonics effects to the cavity transmission profile, identifying polariton cross-talk effects at any temperature. Furthermore, we map the combined effect of vibronic and electron-phonon interactions onto the cavity output field.

A. Cavity transmission
We will consider a cavity mode which is driven with amplitude η c and start with a set of coupled Lagevin equations for the electric field operator a as well as the polaron operatorσ(t) = D † (t)B † (t)σ(t) in a rotating frame at the laser frequency ω ℓ : where we defined the effective cavity input A in = η c / √ 2κ+ a in with zero-average input noise a in but non-vanishing correlation ⟨a in (t)a † in (t ′ )⟩ = δ(t − t ′ ). The electronic transition is also affected by a white noise input σ in with non-zero correlation ⟨σ in (t)σ † in (t ′ )⟩ = δ(t − t ′ ). We can formally integrate Eq. (31b) and substitute it in Eq. (31a): where we took the averages and assume factorizability between optical and vibronic/phononic degrees of freedom which is valid if the timescales of vibrational relaxation and cavity decay are separated, e.g. Γ m ≫ κ. We notice that the second term in Eq. (32) represents a convolution since the correlation functions ⟨D(t)D † (t ′ )⟩ and ⟨B(t)B † (t ′ )⟩ only depend on the time difference ⟩, the normalized cavity transmission amplitude T (ω) = ⟨Aout(ω)⟩ ⟨Ain(ω)⟩ can be derived from input-output relations as where H(ω) is the Fourier transform of H(t) and describes the optical response of the molecule to the light field including electron-phonon, electron-vibron and vibronphonon coupling. If we neglect electron-phonon interactions (λ e-ph = 0) and assume, for the sake of simplicity, Markovian decay for the vibration (this can also be extended to the non-Markovian regime, see Section (IV)), H(ω) acquires the form Again, the above expression indicates a series of sidebands with strength determined by the Huang-Rhys factor λ 2 and dependent on the thermal occupationn. In the case of large vibrational relaxation Γ m ≫ γ (corresponds to typical experimental situation), however, those sidebands are suppressed and the cavity will mostly couple to the ZPL transition. We can then define an effective Rabi coupling for the ZPL which takes into account the reduction of the oscillator strength due to Franck-Condon and Debye-Waller factors. In Figs. 6(a) and (b) we plot the cavity transmission at resonance ω c = ω ℓ for increasing thermal occupation with and without the influence of phonons and find that the splitting of the polariton modes is well described by Eq. (35). This also manifests itself in the transmission signal in the Purcell regime characterized by weak coupling g < |κ − γ|/2, but large cooperativity C = g 2 /(κγ) ≫ 1 which is a more realistic regime in currently available single-molecule experiments [12,16]. In Figure 6(c) we compare the transmission of a pure two-level system (obtained by setting λ = 0, λ e-ph = 0) with a molecule in a solid-state environment.
Here the ZPL appears as a dip in the transmission profile with an increase in width γ = γ(1 + C eff ) proportional to the effective cooperativity C eff = g 2 eff /(κγ). As compared to the two-level system case, the coupling to vibrons and phonons leads to a reduction in both width and depth of the antiresonance. If the cavity bandwidth is comparable to the maximum phonon frequency ω max , the imprint of the phonon wing can also be detected in the transmission signal of the cavity, which is relevant for plasmonic scenarios characterized by large bandwidths [60] [see Fig. 6(c)]. The sidebands of vibrons typically lie at frequencies outside the bandwidth of the cavity ν ≫ κ and are consequently unmodified.

B. Vibrationally mediated polariton cross-talk
As shown in the previous sections, vibronic and electronphonon couplings reduce the oscillator strength of the molecule and lead to decoherence and are consequently considered as detrimental. However such couplings can also lead to interesting new physics: In Ref. [15] it was already shown that vibrations can couple upper and lower polaritonic states in a dissipative fashion, resulting in an effective transfer of population from upper to lower polariton and consequently an asymmetric cavity transmission profile with a suppressed upper polaritonic peak (at zero temperaturen = 0) and dominant emission occuring from the lower polariton (this can also be seen in Figs. 6(a) and (b)). We derive here a more general expression for the population transfer between polaritons showing that for finite thermal occupations of the vibrational moden also a transfer from lower to upper polariton can be activated.
This can be interpreted as an exchange interaction which is mediated by either the destruction or creation of a vibrational quantum. From this one can derive equations of motion for the populations of upper and lower polaritonic states:Ṗ with the hybridized decay rates of upper and lower polaritonic state γ ± = (κ + γ)/2 and one can see that the term ℑ ⟨U † L(b † + b)⟩ is the one responsible for population transfer between the polaritons. In the limit of fast textcolorbluevibrational relaxation Γ m ≫ κ this can be turned into a set of rate equations with an effective excitatation transfer from the upper polariton to the lower polariton κ + and a transfer from the lower polariton to the upper polariton κ − (for detailed calculation see Appendix J): Under the assumption of weak vibronic coupling as compared to the splitting between upper and lower polaritonic state λν ≪ 2g, the rates can be calculated to first order as (again we assume Markovian decay for the vibration for that sake of simplicity): Energy transfer between the polaritons can consequently occur if the Rabi splitting ω + − ω − ≈ 2g is roughly equal to the vibrational frequency. In the case of zero temperature (n = 0) the above equations reduces to the results presented in [15] using a Lindblad decay model for the vibration instead of a Brownian noise model. The ratio κ − /κ + =n/(n + 1) which can be inferred from the polariton heights (for normalized Lorentzians the height and width are connected) and which tends to unity in the limit n ≫ 1 can be seen as a direct measure for temperature as it does not depend on any other parameters. While for single molecules the condition ω + − ω − ≈ ν is difficult to achieve for vibrational modes in the THz range, this can be achieved in the collective strong coupling regime for many molecules where the coupling grows as g √ N or for single molecules with phononic modes in the GHz regime. We also note that, in a similar fashion to the linear coupling, also quadratic electron-phonon and vibronic coupling gives rise to a vibrationally-mediated polariton cross-talk with coupling H quad int = β(U † L + L † U )Q 2 /2 (for a single vibrational mode). To this end, one could again derive effective rate equations for the quadraticallymediated population transfer between the polaritons in a similar fashion as for the linear coupling case.

VI. DISCUSSIONS. CONCLUSIONS
We have provided a new approach based on quantum Langevin equations for the analysis of the fundamental quantum states of molecules and their coupling to their surroundings. These features, which lie at the heart of molecular polaritonics, go well beyond the electronic degrees of freedom and address phenomena such as electronvibron and electron-phonon couplings as well as vibronphonon interactions resulting in the relaxation of molecular vibrations. In particular, we have provided analytical expressions for spectroscopic quantities such as molecular absorption and emission inside and outside optical cavities in the presence of vibrations and phonons at any temperature. Moreover, we have presented a model of vibrational relaxation that takes into account the structure of the surrounding phonon bath and makes a distinction between Markovian and non-Markovian regimes. We have demonstrated that the vibrational relaxation of a molecule is crucially determined by the structure of the bath, especially by the maximum phonon frequency ω max . For two molecules embedded in the same crystalline environment, we have shown that the vibrational modes of the spatially separated molecules can interact with each other, resulting in collective dissipative processes that allow for weaker relaxation of collective vibrations. In the strong coupling regime of cavity QED, we have derived temperature-dependent transfer rates for vibrationally mediated cross-talk between upper and lower polaritonic states, i.e. hybrid light-matter states that are normally uncoupled in cavity QED studies of atomic systems. In this work, we based our model on first-principle derivations of the relevant coupling strengths between a single nuclear coordinate of a molecule embedded in a 1D chain. However, the calculations could be readily extended to 3D scenarios and compared with ab-initio calculations for real materials. We point out that our theory could also be relevant for vacancy centers in diamond where similar interactions between electronic degrees of freedom and both localized and delocalized phonon modes occur [85,86]. In the future we want to address the influence of higher-order interactions such as quadratic electronphonon and vibron-phonon couplings, which are known to play an important role at elevated temperatures. It could also be interesting to consider the cavity-modification of the nonradiative relaxation of molecules at conical intersections [42]. We also plan to investigate the collective radiation states of dense molecular ensembles in confined electromagnetic environments such as e.g. occuring in organic semiconductor microcavities.
Note added. Recently, we became aware of a related study [87].

VII. ACKNOWLEDGMENTS
We acknowledge financial support from the Max Planck Society and from the German Federal Ministry of Education and Research, co-funded by the European Commission (project RouTe), project number 13N14839 within the research program "Photonik Forschung Deutschland" (C. S, V. S. and C. G.). We are mostly interested in the coupling between the phonons and the relative motion of the molecule (which will give rise to vibrational relaxation) as well as the coupling between the electronic excitation and the phonons. In the expanded equations above we can see that the first term couples couples the relative motion of the molecule x rm to the motion of the lattice atoms x ℓ , while the second term couples x ℓ to the electronic excitation σ † σ and the last term is a constant energy shift. For left and right atom together we then arrive at the interaction Hamiltonians where for simplicity we assumed equal masses m L = m R and used that for symmetric molecules ∆x L = −∆x R .

Bulk solution
The contribution of the bulk adds as (just considering nearest-neighbour interactions) with the main assumption that the presence of the molecule does not significantly change the bulk modes. Newton's equations of motion for the displacements of the lattice atoms can be arranged in matrix form as where x = (x 1 , . . . , x 2N +1 ) and where we assume a periodic lattice and neglect the edge modes in our analysis. This is a trigonal Toeplitz matrix which can be easily diagonalized (D = T −1 MT) where the k-th eigenvalue is given by with k ∈ [1, 2N + 1] and the associated normalized eigenvectors read such that in the diagonal basis (u = T −1 x) the equation of motion for the k-th mode is given bÿ We can now express the couplings from Eqs. (A9) in terms of the normal modes of the crystal. Further more we will consider only nearest-neighbour coupling for the molecule, i.e. k Rℓ , k Lℓ = 0 for ℓ ̸ = {N, N + 2}. We obtain with and consequently where we used that ∆x N +2 = −∆x N and defined ∆k = k R(N +2 We now quantize the Hamiltonians by introducing the position operators (and analog momentum operators) where we have introduced the couplings We have further more introduced the dimensionless quadratures [q k , p k ′ ] = iδ kk ′ and the Hamiltonian of the bulk can be expressed as H bulk = k ω k p 2 Figure  8 shows a plot of the dispersion ω k as well as the vibron-phonon coupling α k (the electron-vibron coupling shows a similar behavior).

Quadratic interaction
In the previous sections we always assumed identical curvatures of electronic ground and excited states. Here we show that non-identical curvatures give rise to a quadratic electron-vibron interaction (similar for the electron-phonon interaction). Assuming different vibrational frequencies of electronic ground (ν) and excited state (ν) the molecular Hamiltonian [Eq. (A6)] becomes: Quantizing with x rm = 1/(2µν) b † + b and p rm = i (µν)/2 b † − b , this gives rise to an electron-vibron interaction of the form where the second part accounts for the squeezing of the vibrational wavepacket when going from ground to excited state with coupling strength β = (ν 2 − ν 2 )/(2ν).

Appendix B: Vibrational relaxation
The evolution of the molecule's vibrational mode can be calculated from Hamilton's equations of motioṅ The evolution of the phonon bath degrees of freedom is goverend bẏ We can derive an effective equation of motion for Q and P by eliminating the bath degrees of freedom. We can write Eqs. (B2) in matrix formv = Mv + v inhom with v = (q k , p k ) T and v inhom = (0, α k Q) T . The solution is given by with M = TΛT −1 , Λ = diag(iω k , −iω k ) and With this we can derive the solution for q k (t) and subsequently forṖ (after integration by parts): where we introduced the fluctuating force ξ(t) = k (α k q k (0) cos(ω k t) + α k p k (0) sin(ω k t)). One can easily check that the second term in the sum above sums to zero in the continuum limit while the first term is a renormalization of the vibrational frequency withν = ν − ν s . Using the expressionfor α k derived in the previous section we can calculate (for N ≫ 1). All together the vibrational relaxation can be written aṡ where * denotes the convolution (Γ * P )(t) = ∞ 0 dt ′ Γ(t − t ′ )P (t ′ ) and In the continuum limit this can be approximated by where J n (x) are the Bessel functions of first kind. The random force has zero average but non-vanishing correlations where we denote by ⟨•⟩ = Tr [• ρ th ] the thermal average. The commutator reads In the continuum limit N → ∞ (using the prescription k → L π dq → dωn(ω) where πk 2N +1 = qL 2N +1 = qa) we can rewrite the correlation function as The density of states in 1D is given by Multiplying this with the vibron-phonon coupling in frequency domain gives Using this we can write the correlation function as where we denote by Γ m the Markovian decay rate Γ m = 2ννs ωmax . We note that Γ m actually increases with ω max (in our derivation we assumed the reduced mass of the molecule to be equal to the lattice atoms µ ≈ m 0 ): Γ m = ∆k 2 We can see that we obtain an effective frequency-dependent decay rate Γ r (ω) = Γm √ ω 2 max −ω 2 ωmax Θ(ω max − ω)Θ(ω max + ω) (the real part of the Fourier transform of Γ(t)). In the limit ω max → ∞ this becomes the standard result for a harmonic oscillator in a thermal bath Appendix C: Collective vibrational relaxation We now consider two molecular impurities inside the 1D chain located at positions N + 1 + j and N + 1 − j (distance 2j). Again we assume that the presence of the two molecules does not significantly change the bulk modes. The Hamiltonian reads: From this we can calculate equations of motion for the molecular coordinateṡ as well as for the bath degrees of freedomq (C7) Following the same procedure as in Appendix B we can solve for q k to obtain Plugging this into the equations of motion for the system variables we obtaiṅ where we introduced the input noises ξ 1/2 (t) = k α k,1/2 (cos(ω k t)q k (0) + sin(ω k t)p k (0)). Again integrating by parts we findṖ The two noise terms are correlated according to We now still have to determine the coupling coefficients α k,1 and α k,2 . Assuming that the two molecules are placed at q N +1−j and q N +1+j , the couplings are then determined by the neighboring atoms q N +1−j±1 and q N +1+j±1 : where we denote the neighboring couplings for the first molecule by ±1 and for the second molecule by ±2. With this we can calculate the coupling coefficients (using that sin α − sin β = 2 cos( α+β 2 ) sin( α−β 2 )): With this we obtain for the product α k,1 α k,1 (using that cos(α) cos(β) = [cos(α + β) + cos(α − β)]/2): We can finally write Eqs. (C11) as (neglecting the terms that sum to zero in the continuum limit) where we defined Ω = k Let us for simplicity assume two identical molecules ν 1 = ν 2 . We can write the collective decay as In the continuum limit this can be approximated by Taking the Fourier transform of Eqs. (C17) we obtain a set of algebraic equations (again assuming identical molecules) The Fourier transform of Γ 12 can be obtained in two steps. The Fourier transform of the productf 4j (t)Θ(t) wherẽ f 4j (t) = f 4j (ω max t) = ω max (J 4j (ω max t)/(ω max t)), is given by where we have used that F(f 4j )(ω) = (1/ω max )F(f 4j )(ω/ω max ). Since the Fourier transform of the term containing the Bessel function results in where U n−1 is a Chebyshev polynomial of the second kind, we obtain for the collective decay component where y = ω ′ /ω max and x = ω/ω max . For −1 < x < 1 meaning that −ω max < ω < ω max we have with T n being a Chebyshev polynomial of the first kind [88], we obtain and with it In the case ω max → ∞, we obtain F(Γ 12 * P i )(ω) ≈ Γ m P i (ω) which results in the equations of motioṅ Considering the coordinates Q + = Q 1 + Q 2 , P + = P 1 + P 2 and Q − = Q 1 − Q 2 , P − = P 1 − P 2 , we obtain equations of motion for two independent oscillatorsQ where the first one is protected from vibrational decay.

Appendix D: Fundamental vibron-phonon processes
To investigate the evolution of a single vibron state we have to develop an optimal framework first. Starting with the Hamiltonian we go to the interaction picture following where U = e iH0t and H 0 = νb † b + N k=1 ω k c † k c k . With the initial state in the Schrödinger picture given by |1 ν , vac ph ⟩ we obtain |ϕ⟩ = U |1 ν , vac ph ⟩ = e iνt |1 ν , vac ph ⟩. Following the Schrödinger equation i∂ t |ϕ⟩ =H |ϕ⟩ we acquire the Dyson series Evaluating the first order of the Dyson series we obtain In the case where we have a resonant component ω j = ν with j ∈ {1, . . . , N } and since lim ω→0 (e iωt − 1)/ω = it we obtain which can emerge in the case where ν ∈ {0, . . . , ω max }. No resonance condition can be fulfilled in the case where ν > ω max and the off resonant terms of Eq. (D4) describe the dynamics. Besides the single phonon resonances in the case ν ≤ ω max also weak multi-phonon resonances can be found. For our initial state with one excitation in the vibrational state these terms can be found starting from the third order term of the Dyson expansion. For example for These terms are small with respect to the resonances starting in the first order in the case where ν ≤ ω max and in total are comparable to the off resonant terms.
We denote by J(ω) = n(ω)λ(ω) 2 ω 2 = k |λ k ω k | 2 δ(ω − ω k ) the spectral density of the electron-phonon coupling. From Eq. (A19c) we can derive for the 1D case: which leads to a divergence of the integral at small frequencies due to the 1/ω-term. (in 1D the spectral density will always be proportional to ω for low frequencies). Considering a 3D density of states instead: with group velocity v g = ∂ω/∂q and q(ω) the inverted dispersion relation in 3D (we make the Debye assumption and assume a linear dependence between wavevector and frequency). With this we obtain for the spectral density where we introduced the 3D electron-phonon coupling constant λ 3D e-ph = 16 which has units [s 2 ]. Figure 9 shows the schematic behavior of the spectral densities in 1D and 3D.
Let us first consider the purely vibrational part and ignore the electron-phonon part (corresponds to setting all λ k = 0). We can rewrite The expectation value of the coherence in steady state becomes Using that (n + 1)/n = e βν we can also rewrite the displacement correlation function Eq. (G3) as The generating function of the modified Bessel functions I n (x) is given by e which is similar to the expression known from Huang-Rhys theory [49]. Including the phonon modes, the collective phonon mode displacement correlation for N phonon modes can be expressed as (neglecting phonon decay): λ n k n k ! e −λ 2 k (1+2n k ) n k l k (n k + 1) Together with the the vibrational modes, the total absorption spectrum can then be expressed in terms of discrete lines as P e = η 2 ℓ γ ∞ n=0 n l=0 λ 2n n! e −λ 2 (1+2n) n l (n + 1) n−lnl ∞ n1,...,n N n1,...,n N l1,...,l N N k=1 L k (n k )B k (n k , l k ) γ + n Γm 2 γ + n Γm where we introduced L k (n k ) = λ 2n k n k ! e −λ 2 k (1+2n k ) and B k (n k , l k ) = n k l k (n k + 1) n k −l kn l k k .