Nonlinear opto-vibronics in molecular systems

We analytically tackle opto-vibronic interactions in molecular systems driven by either classical or quantum light fields. In particular, we examine a simple model of molecules with two relevant electronic levels, characterized by potential landscapes with different positions of minima along the internuclear coordinate and of varying curvatures. Such systems exhibit an electron-vibron interaction, which can be comprised of linear and quadratic terms in the vibrational displacement. By employing a combination of conditional displacement and squeezing operators, we present analytical expressions based on a quantum Langevin equations approach, to describe the emission and absorption spectra of such nonlinear molecular systems. Furthermore, we examine the imprint of the quadratic interactions onto the transmission properties of a cavity-molecule system within the collective strong coupling regime of cavity quantum electrodynamics.


I. INTRODUCTION
Opto-vibrational interactions in molecular systems occur in an indirect fashion as light couples to electronic transitions, which are in turn coupled to the vibrations of nuclei [1][2][3][4].A standard description of electron-vibron interactions, under the Born-Oppenheimer approximation, is given by the Holstein Hamiltonian [5,6] which is a spinboson model linear in the vibrational displacement.Some analytical treatments based on quantum Langevin equations (QLEs) [7][8][9][10][11] have been shown to provide approximate analytical results for this model for a large number of vibrational modes and in the presence of fast vibrational (a) Standard scenario where the excited state potential landscape is a copy of the ground state landscape slightly shifted.Electronic excitation is accompanied by the action of a conditional displacement operator D σ † σ , where σ is the ladder operator from the excited to the ground electronic state.(b) Scenario with unshifted potentials with different curvatures.Electronic excitation is accompanied by a conditional squeezing operation S σ † σ .(c) Combined model where electronic excitation leads to a displacing and squeezing operation.
relaxation typically occurring in both bulk [8,12] and solvent environments [13].Similar methods have been used in cavity optomechanics [14,15], where cavity-confined quantum light modes are coupled to macroscopic oscillators via the radiation pressure Hamiltonian, to study the strong photon-phonon coupling regime [16,17].
Such theoretical treatments are based on a polaron transformation which allows for the diagonalization of the bare Holstein Hamiltonian [18].This can be understood as a conditional displacement operation, where the electronic state dictates whether or not a displacement in the vibrational subspace should be performed.In consequence, when a photon excites an electronic transition between two copies of the same harmonic potential landscape slightly shifted (see Fig. 1(a), the vibrational state is excited to a coherent state.The underlying assumption here is however that the potential landscapes are identical.In reality it can happen that the curvatures of the two potential energy surface are different, as illustrated in Fig. 1(b): an electronic transition will then be accompanied by a squeezing of the vibrational wave-packet.In such a case the polaron transformation is modified by an operation involving a conditional squeezing operator.Most generally, one can imagine the situation depicted in Fig. 1(c) where the proper diagonalizing transformation involves a conditional displacement followed by squeezing.In optomechanics, this corresponds to a quadratic photon-phonon interaction [19].
We provide here an analytical treatment based on a set of QLEs for effective spin operators dressed by vibrations, which can be solved under some approximations to provide information about emission and absorption spectra.Additionally, we investigate the transmission properties of an optical cavity within the strong coupling regimes of cavity quantum electrodynamics.By studying the interaction between the molecular systems and the cavity, we gain insight into the nature of light-matter interactions in these complex environments.
The paper is organized as follows: in Sec.II we introduce the modified Holstein model obtained from first principle derivations of the electron-vibration coupling for a scenario depicted in Fig. 1(c).Our analytical treatment is based on a set of simplified QLEs for vibrations and electronic degrees of freedom as derived in Sec.III.We proceed with solving the QLEs under the approximation of weak excitation of the upper electronic state to obtain absorption and emission spectra under illumination with classical light.Finally, in Sec.IV we add a quantum confined light field coupled to the electronic transition via the Tavis-Cummings Hamiltonian and derive the transmission profile of the cavity in the weak and strong coupling regimes of light-matter interactions.

II. THE MODIFIED HOLSTEIN MODEL
We consider a molecule with two relevant electronic states denoted by |g⟩ and |e⟩ for ground and excited, respectively.Transitions between these two states are characterized by Pauli lowering operators σ = |g⟩ ⟨e| and its corresponding Hermitian conjugate.As illustrated in Fig. 1(c), the ground/excited potential landscapes are assumed to have a parabolic shape, with the minima of these two potential landscapes separated by R eg and with different curvatures, thus having different vibrational frequencies: ν g for the electronic ground state and ν e for the electronic excited state.The Hamiltonian describing the molecular system can be expressed as (ℏ = 1) where V e and V g denote the potential landscapes in the electronic excited and ground state, respectively, defined onto the direction of the nuclear coordinate as respect to the ground state such that it diagonalizes the ground state vibrational problem.The Hamiltonian in Eq. ( 1) can now be written as The linear coupling parameter results from the mismatch in the positions of the minima λ 1 = −µν 2 e R eg R zpm /ν g while the quadratic coupling parameter is proportional to the relative change in vibrational frequencies λ 2 = ν 2 e − ν 2 g / 4ν 2 g .The bare electronic frequency splitting is modified by the vibronic coupling ω 0 = ω e − ω g + λ 2  1 ν 3 g /ν 2 e .However, it is more convenient to use a single basis formulation where only the eigenstates of the harmonic oscillator in the ground state are considered, i.e., the eigenstates of ν g b † b denoted by {|m g ⟩}.To this end, one can take the level-dependent unitary transformation H = U † HU with The definitions of the displacement and squeezing operators are the standard ones employed in quantum optics which employ the following displacement r d and squeezing r s parameters defined as and r s = 1 2 (ln ν e − ln ν g ) .
Finally, the Hamiltonian is expressed in diagonal form where the effective frequency ω 00 = ω e − ω g + (ν e − ν g )/2 relates to the zero-phonon line.This is nothing more than a generalized polaron transformation where the electronic coherence operator σ is dressed by the vibrational modes as σD(r d )S(r s ) via both a displacement and a squeezing operation.This offers a recipe to obtain the intensity of vibronic transitions in the emission and absorption processes.Assuming the molecule initially in the excited state with zero vibrations |e; 0 e ⟩, the probability of ending up in the state |g; m g ⟩ is governed by the overlap between the two vibrational wave functions [see Fig. 2(a)] as where H m (x) are Hermite polynomials, α = tanh r s + 1, and β = (tanh r s ) /2.Similarly, we can find the absorption probability amplitude for the absorption transition |g; 0 g ⟩ → |e; m e ⟩ via the Hermitian adjoint operator with α ′ = tanh r s − 1.
We numerically illustrate the departure from such a statistics with various values of λ 2 in Figs.2(b)-(c).Given the commutator [D(r s ), S(r s )] ̸ = 0, the presence of the product D(r s )S(r s ) renders an asymmetry between the emission event |e; 0 e ⟩ → |g; m g ⟩ and the absorption event |g; 0 g ⟩ → |e; m e ⟩.Also, as a simple check, in the limiting case where λ 2 = 0, i.e., ν e = ν g , both transition strengths follow the same Poissonian distribution e −λ 2 1 λ 2m 1 /m!, as expected, reproducing the mirroring effect of emission and absorption spectra usually exhibited by most molecular transitions.

III. ABSORPTION AND EMISSION SPECTRA
In order to derive spectroscopic quantities, we will assume a continuous wave classical drive coupled to the electronic transition incorporated in the following Hamiltonian with the Rabi frequency η ℓ and laser frequency ω ℓ .Since the molecule is also coupled to the electromagnetic vacuum and additional vibrational relaxation baths, we will make use of open system dynamics methods, first formulated in terms of a master equation.First, we include a spontaneous emission channel with the collapse operator σ at rate γ.In addition, as the electronic transition is modified by the vibrational mode [20][21][22], the influence of the environment onto the dynamics of the vibrational mode can be well described by a collapse operator UbU † at the rate Γ.For numerical investigations, the master equation for the system is given where the standard Lindblad superoperator is written as erator O and a corresponding decay rate γ O .In particular in the polaron transformation ρ = U † ρU, the last term in Eq. ( 11) is going to the familiar form L Γ [b] ρ.The dot stands for the position where the density operator, on which the Lindblad superoperator is applied on, is to be included.
It is convenient, for deriving analytical results, to map the master equation into an equivalent set of QLEs.For any system operator A this can be done as follows [7,23] where O in is the zero-averaged and delta-correlated input noise operator associated with the collapse operator O and γ O is the associated decay rate.
For molecules in solid-state environments or in solvents, the vibrational relaxation rate is usually very large greatly surpassing both γ and η ℓ .Therefore, fluorescence occurs preferentially from the state |e, 0 e ⟩, which lies at the bottom of the excited state manifold: this is generally referred to as Kasha's rule [24].The same mechanism is valid for the absorption process, where absorption occurs from the state |g, 0 g ⟩, the lowest in energy.We will make use of this fast vibrational relaxation to impose a quick timescale for the modification of the bosonic b operators and use their quasi-steady state values in the following.First, however, let us partition the total Hilbert space into two orthogonal subspaces (ground and excited electronic state manifolds) via the following two projection operators P g = σσ † and P e = σ † σ.Let us first pay attention to the dynamical equation in the manifold of P e .For convenience reasons, we introduce a projected bosonic operator b e = UbU † P e acting only in this manifold and more explicitly expressed as In a completely similar fashion, projected operators in the ground electronic state manifold can be defined.Let us introduce the ground state polaron operator via the As above, the new displacement operator is )]P g , and the new squeezing operator is The nonvanishing correlation of the zero-average noise operator is given by ).We are now in the position of reconstructing the full solution of the coherence operator in steady state by summing over the contributions in the ground and excited state manifolds.This can be done by formal integration of Eq. (14b) and Eq.(15b) to obtain a solution for ⟨σ⟩ expressed as Here, we have used the Heaviside step function Θ(t) and the initial value ⟨σ(0)⟩ = 0. Considering that the vibrational mode has a large relaxation rate (i.e., Γ ≫ γ), we then decouple the vibronic and electronic degrees of freedom.The S em m e −m(iνg+Γ)(t−τ ) ⟨P e (τ )⟩ , (17a) Replacing the infinite sums from above back into Eq.( 16) leads to a convolution in time.This can be dealt with by employing a Laplace transformation defined as f (s) = ∞ 0 dt f (t) exp(−st) for a time-dependent function f (t) at t ≥ 0. In such a case, Eq. ( 16) takes a much simpler form with the following functions identified corresponding to emission and absorption events, respectively From these expressions, one can proceed in evaluating analytically the population of the excited state p ss e = lim t→∞ ⟨σ † (t)σ(t)⟩ in steady state (as detailed in Ap-pendix D) The coefficients γ ↑ m (ω) and γ ↓ m (ω) represent the dynamic equilibrium population transfer rates for absorption from the ground state to the excited state |g; 0 g ⟩ → |e; m e ⟩ and emission from the excited to the ground state |e; 0 e ⟩ → |g; m g ⟩ as illustrated in Fig. 3(a).The rates are analytically expressed as Specifically, these rates contribute to the rate equation for the population of the excited state, given by (see Appendix E for detailed derivations): This equation holds true under the condition η ℓ ≪ Γ.
Remarkably, one could also obtain the same expression for the population of the excited state in steady state and compare with full numerical simulations to a very good fit, as illustrated in Fig. 3(b).The parameters are given in the caption and are chosen in close attention to other works [25,26].Additionally, we can employ the pump-probe scenario to analyze the absorption and emission processes.In this scenario, the molecule absorbs a photon at the frequency ω ℓ , transitioning to the excited state |e; m e ⟩ under the resonant condition ω ℓ = ω 00 + mν e .Subsequently, after undergoing fast vibrational relaxation, the molecule emits a photon centered around the frequency ω 00 −m ′ ν g , which can be detected with a modified linewidth γ + m ′ Γ.The absorption and emission profiles are then obtained by summing up the contributions from all possible cases, resulting in Lorentzian profiles represented by γ ↑/↓ m , as shown in Fig. 3(c): Here, the scaling of the vibrational rates has been on purpose exaggerated in order to clearly point out the difference in energies expected for the smaller and higher energy sidebands.The presence of the quadratic electronvibron coupling under realistic conditions, is expected to only slightly break the symmetry between the emission and absorption spectra, as the expected values for λ 2 lie well below in the subunit region.More details on the procedure we have followed for the above derivations is presented in Appendix F and basically follows the quantum regression theorem formalism [23,27].

IV. MOLECULAR POLARITONICS
Let us now ask what is the imprint of the asymmetry between the ground and excited state potential landscapes on the signal of an optical cavity containing such a molecule in the strong coupling regime of cavity quantum electrodynamics.To this end, we consider a single molecule placed within the optical volume of a single mode optical cavity mediating transitions between the ground and excited potential landscapes.Under strong optical confinement conditions, the interaction of light and matter can lead to the production of hybrid quantum states, i.e., polaritons [1,11,[28][29][30][31][32][33][34][35][36] as superpositions of ground or excited electronic states and zero or single photon states.While polaritons are eigenstates solely of the electron-photon interaction Hamiltonian, the intrinsic electron-vibron coupling can provide a mechanism of polariton cross-talk, leading to a unidirectional loss of energy from the higher state to the lower energy state.This has been shown analytically in Ref. [7] for the standard case of identical ground and excited state potential landscapes and found to be most pronounced when the vibrational mode is resonant to the interpolariton frequency splitting.
Let us now consider the case of N molecules inside the spatial extent of a single-mode of a Fabry-Pérot optical resonatoras, illustrated in Fig. 4. The dynamics of a single molecule is governed by the Hamiltonian H from Eq. ( 3).The interaction between the N molecules and the cavity field mode is characterized by the Tavis-Cummings model, consisting of the free cavity field at frequency ω c and with bosonic mode a and the Tavis-Cummings interaction with the unit light-matter coupling strength g and the laser field drive with amplitude η c and frequency ω ℓ .For convenience, we have made the assumption here that all molecules are identical.Let us proceed with a set of effective QLEs for the cavity mode a and the state dependent polaron operators σe,n and σg,n for the nth molecule in the rotating frame at the laser frequency ω ℓ : Here, the total dissipation for the cavity field κ = κ 1 + κ 2 encompasses the losses via both mirrors.The operators A 1,in = η c / √ 2κ 1 + a 1,in describes the input classical field coming through the left mirror η c / √ 2κ 1 and the zeroaverage input noise with the only non-vanishing two-time correlations ⟨a 1,in (t)a † 1,in (t ′ )⟩ = δ(t − t ′ ).Additionally, zero-average input noise comes through the right side mirror with similar correlations ⟨a 2,in (t)a † 2,in (t ′ )⟩ = δ(t − t ′ ) and uncorrelated with the a 1,in (t).
The Markovian limit is achieved under the large relaxation rate condition for vibrational mode, i.e., Γ ≫ κ and Γ ≫ γ.In this case, the approach to treat the vibrations as a local phonon bath is still applicable.By formally integrating the equations for polaron operator, tracing over the cavity mode as well as electronic degrees of freedom and taking the Laplace transformation, we have The coupling between the cavity mode a and the projection operator P e,n leads to non-linear effects.However, we restrict our analysis to the weak excitation regime, i.e., the cavity photon number is much smaller than unity and the population of the excited electronic state |e⟩ is negligible (under the condition that η c ≪ κ).In other words, this approximation allows for the construction of  a linear response theory formalism where the transmitted light gives information on the position and linewidths of the hybrid light-matter eigenstates of the system.In the case of identical conditions, the expectation value of the electronic coherence operator σ n for the n-th molecule will be equivalent to that of the other molecules, i.e., ⟨σ⟩ = ⟨σ n ⟩ = ⟨σ m ⟩ (m ̸ = n).Then, the equations of motion are written in the vector form (in the Laplace transform domain) as Mv + v c = 0, with the drift matrix and the definitions v = (⟨a⟩ , ⟨σ⟩) T and v c = (η c , 0) T .The diagonalization of the drift matrix (under resonance condition ∆ c = 0) yields the frequencies ω ± and linewidths γ ± [33,37] of the two polaritons as These particularities of the polaritons can be explored in a very simple way by performing a scan of the laser frequency around the cavity resonance and noticing the position of the peaks corresponding to the hybrid lightmatter states.This can be done at the analytical level in the weak excitation regime and compared to full exact numerics.We define the complex cavity transmission amplitude as the ratio of the normalized continuous outgoing field versus incoming field amplitudes and illustrate its behavior with respect to the scanning laser frequency in Fig. 5.The quantity ⟨a⟩ ss is the average value of the cavity mode amplitude in steady state in the linear response regime with χ ab = lim s→0 1/G ab .
We illustrate numerical and analytical results in Fig. 5 where the profile of the cavity transmission at ω c = ω 00 is plotted.The presence of the linear electron-vibron coupling scaling with λ 1 induces an interaction between upper and lower polaritons already presented at the theoretical level in a few treatments [1,7,8].Instead, at the level of a single molecule, the quadratic interaction will suppress the polariton cross talk, as illustrated in Fig 5(a).In essence, the squeezing term is responsible with a shift in the molecular resonance which then in turn brings the cavity off-resonance with the electronic transition except the zero-phonon transition process.Increasing the number of molecules while assuming very weak driving conditions presents a different situation.This is shown in Fig. 5(b) as an effective reduction of the upper polariton with increasing particle number.

V. CONCLUSIONS
We have applied the toolbox of open system dynamics and in particular the QLEs formalism to analytically describe spectroscopic properties of solid-state embedded molecules in free space or in optical cavity settings.In particular, we generalized our previous approach introduced in Ref. [7] to a scenario where the potential landscapes of a molecule have unequal curvatures in the ground and excited electronic state.This has seen the introduction of a generalized polaron operator where the electronic degree of freedom is dressed by vibrations via a displacement operation followed by an additional squeezing operation.The first effect is seen in the emergent asymmetry between absorption and emission profiles for molecular spectroscopy.A second effect that emerges from our analytical calculations is the context of cavity quantum electrodynamics where the additional squeezing operation leads to a detuning between the bare molecular resonance and the cavity resonance.Our calculations can be relevant in the direction of optomechanics or optovibronics, owing to the strong electron-vibron couplings, albeit under very lossy conditions.Under the condition Γ ≫ γ, and Γ ≫ η ℓ , one can assume the correlation between noise operators B in e and B in g is negligible as and the nonvanishing correlation functions obeying the fluctuation-dissipation relation Effective Quantum Langevin Equation for the Electronic Transition.
Let us pay attention to the electronic transition by introducing the "dressed" dipole operator, i.e., polaron operator, ) in a rotating frame at driving frequency ω ℓ can be expressed as Notably, Eq. (B8) only conclude the contribution of the vibrations projecting to the manifold of P e .To get the dynamics of the system in the whole Hilbert space, one need also to get the dynamics of the general polaron operator σg = S † g D † g σ for the vibrational mode projecting to the manifold of P g with S g = exp r s (b g − b †2 g )/2 P g , and Repeat the process to deriving the dynamical equation of σe , one can recast the above equation into with σg = exp (iν e − iν g ) b † g b g t σg .Meanwhile, one could also get the dynamical equation for the population of the excited state given by Appendix C: Vibrational Dynamics As discussed in the previous section, the electronic transition is dressed by vibrations.For the further calculation of the electronic transition, we here analyze the properties of vibrations and derive the expression for the two time-correlation terms for the product of squeezing and displacement operators S † g D † g , S † e D † e .
In the previous section, we have derived the effective Langevin equation (B4) describing the dynamics of the vibrational mode projected onto the manifold of P e .One can easily obtain the exact solution of such an equation, which is given by Due to the correlation time for vibrations being much short than the electrons, the quantity for σ † σ(τ ) varies little around σ † σ(t).We can thus proceed via a Markov approximation with taking it out of the integral, which yields With this, we can derive the two-time correlation for t > τ where ⟨•⟩ vib denotes taking the average over the degrees of freedom of the vibrational mode.
Similar to the case for calculating b e , we can obtaint the two-time correlation term for the vibrational operator b g projected onto the manifold of P g as Thus, the two-time correlation function for the vibrational mode in the whole Hilbert space reads Here, the correlations between different manifolds have been properly dropped because of the much smaller value for such terms ⟨b e (τ )b † g (0)⟩ and ⟨b g (τ )b † e (0)⟩ comparing with that in the same manifold, as illustrate in Fig. (6).The squeezing operator can be written in a disentangled form we can calculate the two-time correlation function for displacement-squeezing operator under the vacuum state.For instance, the two-time correlation for the manifold of P e is given by The two-time correlation in the manifold of P g reads and also where S ab m e −m(iν+Γ)t , (C14a) with α = tanh r s − 1 = −2ν g /(ν e + ν g ).

Appendix D: Stability Analysis
It is straightforward to obtain the expression for the Pauli operator in the whole Hilbert space by formally integrating Eqs.(B10) and (B12), with the initial value σ(0) = 0.Under the assumption that the correlation time for the vibrations is much shorter than that for the electronic transition, we can treat the vibrations as a Markovian phonon bath.By taking the average of the vibrational mode and substituting Eq. (C11) as well as Eq.(C13) into the equation above, we then obtain , where Θ(t) is the Heaviside step function.
Tracing over the electronic transition, we finally obtain the simplified formal solution for ⟨σ⟩ as with the population on the electronic excited state P e = ⟨σ † σ⟩.We proceed our calculation via the Laplace transformation (defined as f (s) = ∞ 0 dt f (t) exp(−st) for a timedependent function f (t) at t ≥ 0).Then Eq. (D3) can be written in the Laplace domain as where G em and G ab are the Laplace transform of G ab (ν e , t) and G em (−ν g , t), expressed as Assuming that the molecule is prepared in the electronic ground state |g⟩ followed by taking an average over the electronic transition on both sides of Eq. (B13), and applying the Laplace transformation, we finally obtain Plugging Eq. (D6) into Eq.(D3), we can get where I[•] and R[•] denote taking the imaginary and real part, respectively.According to the final value theorem, we get the steady values with χ Q = lim s→0 G Q .In the limit of weak driving η ℓ ≪ γ, the equation above can be simplified to Appendix E: Rate Equation For large vibrational relaxation rates Γ ≫ η ℓ , the electronic transition is usually going from the lowest vibrational state in both electronic states |g⟩ and |e⟩.Thus, the motion of the population p e m on the state |e; m e ⟩ as well as the population p g m on the state |g; m g ⟩ are given phenomenology by where γ m describe the incoherent spontaneous emission progress |e; 0⟩ → |g; m⟩ satisfying  We now assume the molecule is initially prepared in the excited state |e; 0⟩ to compute the spectrum of emission in the transient regime.By setting η ℓ = 0, one could obtain the expression of the two-time correlation function ⟨σ † (0)σ(τ )⟩ through Eq. (F3) ⟨σ † (0)σ(t)⟩ = F em (ν g , τ )e −(iω00+γ)t . (F4) Taking the Fourier transformation gives the expression of the emission spectrum S em m (γ + mΓ) (γ + mΓ) 2 + (ω 00 − ω − mν g ) 2 . (F5) Absorption Spectra in the Stability Regime.

(F6)
For simplicity, let us assume the amplitude for the driving field is very weak, so that only few molecules are occupy their excited state, i.e., P e (t) ≪ 1 and P g (t) ≈ 1.The solution will become (F11)

Figure 1 .
Figure 1.(a)Standard scenario where the excited state potential landscape is a copy of the ground state landscape slightly shifted.Electronic excitation is accompanied by the action of a conditional displacement operator D σ † σ , where σ is the ladder operator from the excited to the ground electronic state.(b) Scenario with unshifted potentials with different curvatures.Electronic excitation is accompanied by a conditional squeezing operation S σ † σ .(c) Combined model where electronic excitation leads to a displacing and squeezing operation.

Figure 2 .
Figure 2. (a) Schematic diagram of a molecular system exhibiting two parabolic electronic potential surfaces, slightly shifted and with different curvatures quantified by the vibrational frequencies νg and νe.Histogram of the vibrational state occupancy in the electronic ground state upon emission from |e, 0e⟩ in (b) and in the electronic excited state upon external drive from the |g; 0g⟩ state in (c) for various values of λ2 at fixed λ1 = 1.
) and obeying the relation b † e b e |e; m e ⟩ = m e |e; m e ⟩.Meanwhile, we define a time-dependent generalized polaron operator [7, 8], by the transformation σe = σS † e D † e exp i(ν e − ν g )b † e b e t .This allows the derivation of a set of effective QLEs in the rotating frame at the driving frequency ω ℓ for the emission process (see Appendix B for details) ḃe ≈ − (iν e + Γ) b e + √ 2ΓB in e P e , (14a) σe ≈ − (i∆ ℓ + γ) σe − η ℓ S † e D † e e i(νe−νg)b † e b e t + 2γσ in S † e D † e e i(νe−νg)b † e b e t , with the detuning ∆ ℓ = ω 00 − ω ℓ , the displacement operator D e = exp r d (b † e − b e ) P e and the squeezing operator S e = exp r s (b 2 e − b †2 e )/2 P e .The input noises B in e and σ in are zero-averaged and have the following two-time correlations ⟨B in e two-time correlation functions on the right side of above equation could be expressed as (see Appendix B for details) ⟨D e (τ )S e (τ )e i(νe−νg)b † e b e τ e −i(νe−νg)b † e b e t S † e (t)D † e (t)⟩ = ∞ m=0

Figure 4 .
Figure 4. Schematics of an ensemble of molecules inside a Fabry-Pérot resonator.Cavity photon loss occurs at rates κ1 and κ2 via the mirrors M1 and M2, respectively.Lightmolecule interactions occur at rate g while spontaneous emission and cavity driving are at rates γ and ηc, respectively.

Figure 5 .
Figure 5. Cavity transmission (|T | 2 ) of the molecule with various linear and quadratic electron-vibron couplings λ1 and λ2 for a strong coupling to a single cavity mode.(a) Cavity transmission with a single molecule.(b) Cavity transmission as the function of the number of molecules in the cavity.The white lines represent the profile at N = 25 and 100, respectively.Parameters: ωc = ω00, g = 3κ, νg = 10κ, γ = 0.01κ, Γ = 20κ, 2κ1 = 2κ2 = κ, and the driving is assumed very weak ηc/κ = 0.001.

s→0 1 /
G ab and ∆ eff = I lim s→0 1/G ab denoting the effective decay rate and additional frequency shift.
Two-time Correlation Function for the Displacement-Squeezing Operator.
Emission Spectra in the Transient Regime.