# Nonlinear opto-vibronics in molecular systems

Q. Zhang  Max Planck Institute for the Science of Light, D-91058 Erlangen,
Germany  M. Asjad  Department of Mathematics, Khalifa University, Abu Dhabi
127788, United Arab Emirates  M. Reitz  Department of Chemistry and
Biochemistry, University of California San Diego, La Jolla, California 92093,
USA  C. Sommer  Alpine Quantum Technologies GmbH, 6020 Innsbruck, Austria  B.
Gurlek  Max Planck Institute for the Structure and Dynamics of Matter and
Center for Free-Electron Laser Science, 22761 Hamburg, Germany  C. Genes  Max
Planck Institute for the Science of Light, D-91058 Erlangen, Germany
Department of Physics, Friedrich-Alexander-Universität Erlangen-Nürnberg
(FAU), D-91058 Erlangen, Germany

(September 1, 2024)

###### Abstract

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 Neuman and Aizpurua (2018); Clark and Dines
(1986); Smith and Dent (2005); Riede _et al._ (2018). A standard description
of electron-vibron interactions, under the Born-Oppenheimer approximation, is
given by the Holstein Hamiltonian Holstein (1959); Spano (2006) which is a
spin-boson model linear in the vibrational displacement.

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
$\mathcal{D}^{\sigma^{\dagger}\sigma}$, where $\sigma$ 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 $\mathcal{S}^{\sigma^{\dagger}\sigma}$. (c)
Combined model where electronic excitation leads to a displacing and squeezing
operation.

Some analytical treatments based on quantum Langevin equations (QLEs) Reitz
_et al._ (2019, 2020); Zhang _et al._ (2023); Kansanen _et al._ (2021, 2019)
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
relaxation typically occurring in both bulk Reitz _et al._ (2020); Keeling and
Kéna-Cohen (2020) and solvent environments Gilmore and McKenzie (2005).
Similar methods have been used in cavity optomechanics Aspelmeyer _et al._
(2014a, b), where cavity-confined quantum light modes are coupled to
macroscopic oscillators via the radiation pressure Hamiltonian, to study the
strong photon-phonon coupling regime Rabl (2011); Nunnenkamp _et al._ (2011).

Such theoretical treatments are based on a polaron transformation which allows
for the diagonalization of the bare Holstein Hamiltonian Jang (2022). 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 Liao and Nori (2013).

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
$\ket{g}$ and $\ket{e}$ for ground and excited, respectively. Transitions
between these two states are characterized by Pauli lowering operators
$\sigma=\ket{g}\bra{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: $\nu_{g}$ for the electronic ground state and
$\nu_{e}$ for the electronic excited state. The Hamiltonian describing the
molecular system can be expressed as ($\hbar=1$)

$\displaystyle\mathcal{H}=\mathcal{V}_{e}(\hat{R},\hat{P})\sigma^{\dagger}\sigma+\mathcal{V}_{g}(\hat{R},\hat{P})\sigma\sigma^{\dagger},$
(1)

where $\mathcal{V}_{e}$ and $\mathcal{V}_{g}$ denote the potential landscapes
in the electronic excited and ground state, respectively, defined onto the
direction of the nuclear coordinate as

$\displaystyle\mathcal{V}_{e}(\hat{R},\hat{P})$
$\displaystyle=\omega_{e}+\frac{\hat{P}^{2}}{2\mu}+\frac{1}{2}\mu\nu_{e}^{2}(\hat{R}-R_{\mathrm{eg}})^{2},$
(2a) $\displaystyle\mathcal{V}_{g}(\hat{R},\hat{P})$
$\displaystyle=\omega_{g}+\frac{\hat{P}^{2}}{2\mu}+\frac{1}{2}\mu\nu_{g}^{2}\hat{R}^{2},$
(2b)

with the reduced mass $\mu$, the momentum operator and position operators
$\hat{P}$ and $\hat{R}$ satisfying the commutation relation
$[\hat{R},\hat{P}]=i$. Notice that the matrix elements of the Hamiltonian
$\mathcal{H}$ can be written in a basis formed by
$\\{\ket{g;m_{g}}=\ket{g}\otimes\ket{m_{g}},\ket{e;m_{e}}=\ket{e}\otimes\ket{m_{e}}\\}$,
where the Fock states $\ket{m_{e}}$ and $\ket{m_{g}}$, respectively, refer to
the eigenstates of the vibrational Hamiltonian part contained in
$\mathcal{V}_{g}(\hat{R},\hat{P})$ and $\mathcal{V}_{e}(\hat{R},\hat{P})$,
respectively.

However, one can express the quadratures in terms of creation
$b^{\dagger}=(\hat{R}/{R_{\mathrm{zpm}}}-iR_{\mathrm{zpm}}\hat{P})/\sqrt{2}$
and annihilation
$b=(\hat{R}/{R_{\mathrm{zpm}}}+iR_{\mathrm{zpm}}\hat{P})/\sqrt{2}$ operators.
The operators fulfill the following commutation $\left[b,b^{\dagger}\right]=1$
and the zero-point motion is defined as
$R_{\mathrm{zpm}}=1/\sqrt{2\mu\nu_{g}}$. Notice that the definition of this
bosonic operator is performed with respect to the ground state such that it
diagonalizes the ground state vibrational problem. The Hamiltonian in Eq. (1)
can now be written as

$\displaystyle\mathcal{H}=$
$\displaystyle\nu_{g}b^{\dagger}b+\omega_{0}\sigma^{\dagger}\sigma+\lambda_{1}\nu_{g}(b+b^{\dagger})\sigma^{\dagger}\sigma$
(3)
$\displaystyle+\lambda_{2}\nu_{g}(b+b^{\dagger})^{2}\sigma^{\dagger}\sigma.$

The linear coupling parameter results from the mismatch in the positions of
the minima
$\lambda_{1}=-\mu\nu_{e}^{2}R_{\mathrm{eg}}R_{\mathrm{zpm}}/\nu_{g}$ while the
quadratic coupling parameter is proportional to the relative change in
vibrational frequencies
$\lambda_{2}=\left(\nu_{e}^{2}-\nu_{g}^{2}\right)/\left(4\nu_{g}^{2}\right)$.
The bare electronic frequency splitting is modified by the vibronic coupling
$\omega_{0}=\omega_{e}-\omega_{g}+\lambda_{1}^{2}\nu_{g}^{3}/\nu_{e}^{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 $\nu_{g}$ and $\nu_{e}$. Histogram
of the vibrational state occupancy in the electronic ground state upon
emission from $\ket{e,0_{e}}$ in (b) and in the electronic excited state upon
external drive from the $\ket{g;0_{g}}$ state in (c) for various values of
$\lambda_{2}$ at fixed $\lambda_{1}=1$.

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 $\nu_{g}b^{\dagger}b$ denoted by $\\{\ket{m_{g}}\\}$.
To this end, one can take the level-dependent unitary transformation
$\tilde{\mathcal{H}}=\mathcal{U}^{\dagger}\mathcal{H}\mathcal{U}$ with

$\mathcal{U}=\mathcal{D}(r_{d})^{\sigma^{\dagger}\sigma}\mathcal{S}(r_{s})^{\sigma^{\dagger}\sigma}=\sigma\sigma^{\dagger}+\mathcal{D}(r_{d})\mathcal{S}(r_{s})\sigma^{\dagger}\sigma.$
(4)

The definitions of the displacement and squeezing operators are the standard
ones employed in quantum optics

$\displaystyle\mathcal{D}(r_{d})=e^{r_{d}(b^{\dagger}-b)}\quad\mathrm{and}\quad\mathcal{S}(r_{s})=e^{\frac{1}{2}r_{s}\left(b^{2}-b^{\dagger
2}\right)}$ (5)

which employ the following displacement $r_{d}$ and squeezing $r_{s}$
parameters defined as

$r_{d}=-\lambda_{1}\frac{\nu_{g}^{2}}{\nu_{e}^{2}}\quad\mathrm{and}\quad
r_{s}=\frac{1}{2}\left(\ln\nu_{e}-\ln\nu_{g}\right).$ (6)

Finally, the Hamiltonian is expressed in diagonal form

$\tilde{\mathcal{H}}=\nu_{g}b^{\dagger}b\sigma\sigma^{\dagger}+\left(\nu_{e}b^{\dagger}b+\omega_{00}\right)\sigma^{\dagger}\sigma,$
(7)

where the effective frequency
$\omega_{00}=\omega_{e}-\omega_{g}+(\nu_{e}-\nu_{g})/2$ relates to the zero-
phonon line.

This is nothing more than a generalized polaron transformation where the
electronic coherence operator $\sigma$ is dressed by the vibrational modes as
$\sigma\mathcal{D}(r_{d})\mathcal{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 $\ket{e;0_{e}}$, the
probability of ending up in the state $\ket{g;m_{g}}$ is governed by the
overlap between the two vibrational wave functions [see Fig. 2(a)] as

$\displaystyle S_{m}^{\mathrm{em}}=$
$\displaystyle\left|\braket{m_{g}}{0_{e}}\right|^{2}=\left|\bra{m_{g}}\mathcal{S}(r_{s})\mathcal{D}(r_{d})\ket{0_{g}}\right|^{2}$
(8) $\displaystyle=$
$\displaystyle\frac{e^{-r_{d}^{2}\alpha}}{\cosh(r_{s})}\left[H_{m}\left(\frac{\alpha
r_{d}}{2\sqrt{\beta}}\right)\right]^{2}\frac{\beta^{m}}{m!},$

where $H_{m}(x)$ are Hermite polynomials, $\alpha=\tanh r_{s}+1$, and
$\beta=\left(\tanh r_{s}\right)/2$. Similarly, we can find the absorption
probability amplitude for the absorption transition
$\ket{g;0_{g}}\rightarrow\ket{e;m_{e}}$ via the Hermitian adjoint operator
$\sigma^{\dagger}\mathcal{S}^{\dagger}(r_{s})\mathcal{D}^{\dagger}(r_{d})$
such that

$S_{m}^{\mathrm{ab}}=\frac{e^{\alpha^{\prime}r_{d}^{2}\exp\left(2r_{s}\right)}}{\cosh(r_{s})}\left[H_{m}\left(-\frac{i\alpha^{\prime}r_{d}e^{r_{s}}}{2\sqrt{\beta}}\right)\right]^{2}\frac{\left(-\beta\right)^{m}}{m!},$
(9)

with $\alpha^{\prime}=\tanh r_{s}-1$.

We numerically illustrate the departure from such a statistics with various
values of $\lambda_{2}$ in Figs. 2(b)-(c). Given the commutator
$\left[\mathcal{D}(r_{s}),\mathcal{S}(r_{s})\right]\neq 0$, the presence of
the product $\mathcal{D}(r_{s})\mathcal{S}(r_{s})$ renders an asymmetry
between the emission event $\ket{e;0_{e}}\rightarrow\ket{g;m_{g}}$ and the
absorption event $\ket{g;0_{g}}\rightarrow\ket{e;m_{e}}$. Also, as a simple
check, in the limiting case where $\lambda_{2}=0$, i.e., $\nu_{e}=\nu_{g}$,
both transition strengths follow the same Poissonian distribution
$e^{-\lambda_{1}^{2}}\lambda_{1}^{2m}/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

$\mathcal{H}_{\ell}=i\eta_{\ell}\left(\sigma^{\dagger}e^{-i\omega_{\ell}t}-\sigma
e^{i\omega_{\ell}t}\right),$ (10)

with the Rabi frequency $\eta_{\ell}$ and laser frequency $\omega_{\ell}$.
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 $\sigma$ at
rate $\gamma$. In addition, as the electronic transition is modified by the
vibrational mode Mustroph _et al._ (2009); Miao _et al._ (2023); Zirkelbach
_et al._ (2022), the influence of the environment onto the dynamics of the
vibrational mode can be well described by a collapse operator
$\mathcal{U}b\mathcal{U}^{\dagger}$ at the rate $\Gamma$. For numerical
investigations, the master equation for the system is given

$\dot{\rho}=-i\left[\mathcal{H}+\mathcal{H}_{\ell},\rho\right]+\mathcal{L}_{\gamma}[\sigma]\rho+\mathcal{L}_{\Gamma}\left[\mathcal{U}b\mathcal{U}^{\dagger}\right]\rho,$
(11)

where the standard Lindblad superoperator is written as
$\mathcal{L}_{\gamma_{\mathcal{O}}}\cdot=\gamma_{\mathcal{O}}\left(2\mathcal{O}\cdot\mathcal{O}^{\dagger}-\mathcal{O}^{\dagger}\mathcal{O}\cdot-\cdot\mathcal{O}^{\dagger}\mathcal{O}\right)$
for a collapse operator $\mathcal{O}$ and a corresponding decay rate
$\gamma_{\mathcal{O}}$. In particular in the polaron transformation
$\tilde{\rho}=\mathcal{U}^{\dagger}\rho\mathcal{U}$, the last term in Eq. (11)
is going to the familiar form
$\mathcal{L}_{\Gamma}\left[b\right]\tilde{\rho}$. 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 $\mathcal{A}$ this can
be done as follows Reitz _et al._ (2019); Gardiner and Zoller (2004)

$\displaystyle\dot{\mathcal{A}}=$
$\displaystyle-i\left[\mathcal{A},\mathcal{H}+\mathcal{H}_{\ell}\right]-\left[\mathcal{A},\mathcal{O}^{\dagger}\right]\left(\gamma_{\mathcal{O}}\mathcal{O}-\sqrt{2\gamma_{\mathcal{O}}}\mathcal{O}_{\mathrm{in}}^{\phantom{\dagger}}\right)$
$\displaystyle+\left(\gamma_{\mathcal{O}}\mathcal{O}^{\dagger}-\sqrt{2\gamma_{\mathcal{O}}}\mathcal{O}_{\mathrm{in}}^{\dagger}\right)\left[\mathcal{A},\mathcal{O}\right],$
(12)

where $\mathcal{O}_{\mathrm{in}}$ is the zero-averaged and delta-correlated
input noise operator associated with the collapse operator $\mathcal{O}$ and
$\gamma_{\mathcal{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 $\gamma$ and
$\eta_{\ell}$. Therefore, fluorescence occurs preferentially from the state
$\ket{e,0_{e}}$, which lies at the bottom of the excited state manifold: this
is generally referred to as Kasha’s rule del Valle and Catalán (2019). The
same mechanism is valid for the absorption process, where absorption occurs
from the state $\ket{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 $\mathcal{P}_{g}=\sigma\sigma^{\dagger}$
and $\mathcal{P}_{e}=\sigma^{\dagger}\sigma$. Let us first pay attention to
the dynamical equation in the manifold of $\mathcal{P}_{e}$. For convenience
reasons, we introduce a projected bosonic operator
$b_{e}=\mathcal{U}b\mathcal{U}^{\dagger}\mathcal{P}_{e}$ acting only in this
manifold and more explicitly expressed as

$b_{e}=\left(\cosh r_{s}b+\sinh
r_{s}b^{\dagger}-r_{d}e^{r_{s}}\right)\mathcal{P}_{e}$ (13)

and obeying the relation
$b_{e}^{\dagger}b_{e}\ket{e;m_{e}}=m_{e}\ket{e;m_{e}}$. Meanwhile, we define a
time-dependent generalized polaron operator Reitz _et al._ (2019, 2020), by
the transformation
$\tilde{\sigma}_{e}=\sigma{\mathcal{S}}_{e}^{\dagger}{\mathcal{D}}_{e}^{\dagger}\exp\left[i(\nu_{e}-\nu_{g})b^{\dagger}_{e}b_{e}^{\phantom{\dagger}}t\right]$.
This allows the derivation of a set of effective QLEs in the rotating frame at
the driving frequency $\omega_{\ell}$ for the emission process (see Appendix B
for details)

$\displaystyle\dot{b}_{e}\approx-\left(i\nu_{e}+\Gamma\right)b_{e}+\sqrt{2\Gamma}\mathcal{B}^{\mathrm{in}}_{e}\mathcal{P}_{e},$
(14a) $\displaystyle\begin{aligned}
\dot{\tilde{\sigma}}_{e}\approx&-\left(i\Delta_{\ell}+\gamma\right)\tilde{\sigma}_{e}-\eta_{\ell}{\mathcal{S}}_{e}^{\dagger}{\mathcal{D}}_{e}^{\dagger}e^{i\left(\nu_{e}-\nu_{g}\right)b^{\dagger}_{e}b_{e}^{\phantom{\dagger}}t}\\\
&+\sqrt{2\gamma}\sigma_{\mathrm{in}}{\mathcal{S}}_{e}^{\dagger}{\mathcal{D}}_{e}^{\dagger}e^{i\left(\nu_{e}-\nu_{g}\right)b^{\dagger}_{e}b_{e}^{\phantom{\dagger}}t},\end{aligned}$
(14b)
$\displaystyle\dot{\mathcal{P}}_{e}=-2\gamma\mathcal{P}_{e}+\eta_{\ell}(\sigma+\sigma^{\dagger})+\sqrt{2\gamma}(\sigma^{\dagger}\sigma_{\mathrm{in}}^{\phantom{\dagger}}+\sigma_{\mathrm{in}}^{\dagger}\sigma),$
(14c)

with the detuning $\Delta_{\ell}=\omega_{00}-\omega_{\ell}$, the displacement
operator
$\mathcal{D}_{e}=\exp\left[r_{d}(b_{e}^{\dagger}-b_{e})\right]\mathcal{P}_{e}$
and the squeezing operator
$\mathcal{S}_{e}=\exp\left[r_{s}(b_{e}^{2}-b_{e}^{\dagger
2})/2\right]\mathcal{P}_{e}$. The input noises $\mathcal{B}_{e}^{\mathrm{in}}$
and $\sigma_{\mathrm{in}}$ are zero-averaged and have the following two-time
correlations
$\braket{\mathcal{B}_{e}^{\mathrm{in}}(t)\mathcal{B}_{e}^{\mathrm{in}\dagger}(t^{\prime})}=\delta(t-t^{\prime})$
and
$\braket{\sigma_{\mathrm{in}}^{\phantom{\dagger}}(t)\sigma^{\dagger}_{\mathrm{in}}(t^{\prime})}=\delta(t-t^{\prime})$.

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 transformation
$\tilde{\sigma}_{g}=\exp{\left[i\left(\nu_{e}-\nu_{g}\right)b^{\dagger}_{g}b_{g}^{\phantom{\dagger}}t\right]}\mathcal{S}_{g}^{\dagger}\mathcal{D}_{g}^{\dagger}\sigma$
which allows one to derive a similar set of QLEs

$\displaystyle\dot{b}_{g}\approx-(i\nu_{g}+\Gamma)b_{g}+\sqrt{2\Gamma}\mathcal{B}_{g}^{\mathrm{in}}\mathcal{P}_{g},$
(15a) $\displaystyle\begin{aligned}
\dot{\tilde{\sigma}}_{g}\approx&-(i\Delta_{\ell}+\gamma)\tilde{\sigma}_{g}+\eta_{\ell}e^{i\left(\nu_{e}-\nu_{g}\right)b^{\dagger}_{g}b_{g}^{\phantom{\dagger}}t}\mathcal{S}_{g}^{\dagger}\mathcal{D}_{g}^{\dagger}\\\
&+\sqrt{2\gamma}e^{i\left(\nu_{e}-\nu_{g}\right)b^{\dagger}_{g}b_{g}^{\phantom{\dagger}}t}\mathcal{S}_{g}^{\dagger}\mathcal{D}_{g}^{\dagger}\sigma_{\mathrm{in}}.\end{aligned}$
(15b)

As above, the new displacement operator is
$\mathcal{D}_{g}=\exp[r_{d}(b^{\dagger}_{g}-b_{g})]\mathcal{P}_{g}$, and the
new squeezing operator is $\mathcal{S}_{g}=\exp[r_{s}(b^{2}_{g}-b^{\dagger
2}_{g})/2]\mathcal{P}_{g}$. The nonvanishing correlation of the zero-average
noise operator is given by
$\braket{\mathcal{B}_{g}^{\mathrm{in}}(t)\mathcal{B}_{g}^{\mathrm{in}\dagger}(t^{\prime})}=\delta(t-t^{\prime})$.

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 $\braket{\sigma}$ expressed
as

$\displaystyle\braket{\sigma}=$
$\displaystyle-\eta_{\ell}\int_{0}^{\infty}d\tau\,\Theta(t-\tau)e^{-(i\Delta_{\ell}+\gamma)\left(t-\tau\right)}\braket{\mathcal{S}_{e}^{\dagger}\left(\tau\right)\mathcal{D}_{e}^{\dagger}\left(\tau\right)e^{i\left(\nu_{e}-\nu_{g}\right)b_{e}^{\dagger}b_{e}^{\phantom{\dagger}}\tau}e^{-i\left(\nu_{e}-\nu_{g}\right)b_{e}^{\dagger}b_{e}^{\phantom{\dagger}}t}\mathcal{D}_{e}\left(t\right)\mathcal{S}_{e}\left(t\right)}$
(16)
$\displaystyle+\eta_{\ell}\int_{0}^{\infty}d\tau\,\Theta(t-\tau)e^{-(i\Delta_{\ell}+\gamma)\left(t-\tau\right)}\braket{\mathcal{D}_{g}(t)\mathcal{S}_{g}(t)e^{-i\left(\nu_{e}-\nu_{g}\right)b_{g}^{\dagger}b_{g}^{\phantom{\dagger}}t}e^{i\left(\nu_{e}-\nu_{g}\right)b_{g}^{\dagger}b_{g}^{\phantom{\dagger}}\tau}\mathcal{S}_{g}^{\dagger}(\tau))\mathcal{D}^{\dagger}_{g}(\tau)}.$

Here, we have used the Heaviside step function $\Theta(t)$ and the initial
value $\braket{\sigma(0)}=0$. Considering that the vibrational mode has a
large relaxation rate (i.e., $\Gamma\gg\gamma$), we then decouple the vibronic
and electronic degrees of freedom. The two-time correlation functions on the
right side of above equation could be expressed as (see Appendix B for
details)

$\displaystyle\braket{\mathcal{D}_{e}(\tau)\mathcal{S}_{e}(\tau)e^{i(\nu_{e}-\nu_{g})b^{\dagger}_{e}b_{e}^{\phantom{\dagger}}\tau}e^{-i(\nu_{e}-\nu_{g})b^{\dagger}_{e}b_{e}^{\phantom{\dagger}}t}\mathcal{S}_{e}^{\dagger}(t)\mathcal{D}_{e}^{\dagger}(t)}=$
$\displaystyle\sum_{m=0}^{\infty}S_{m}^{\mathrm{em}}e^{-m(i\nu_{g}+\Gamma)(t-\tau)}\braket{\mathcal{P}_{e}\left(\tau\right)},$
(17a)
$\displaystyle\braket{\mathcal{D}_{g}(t)\mathcal{S}_{g}(t)e^{-i\left(\nu_{e}-\nu_{g}\right)b_{g}^{\dagger}b_{g}^{\phantom{\dagger}}t}e^{i\left(\nu_{e}-\nu_{g}\right)b_{g}^{\dagger}b_{g}^{\phantom{\dagger}}\tau}\mathcal{S}_{g}^{\dagger}(\tau))\mathcal{D}^{\dagger}_{g}(\tau)}=$
$\displaystyle\sum_{m=0}^{\infty}S_{m}^{\mathrm{ab}}e^{m(-i\nu_{e}+\Gamma)(t-\tau)}\braket{\mathcal{P}_{g}(\tau)}.$
(17b)

Figure 3: (a) Jablonski diagram illustrating possible emission and absorption
processes. (b) Comparison of analytical and numerical results of the excited
state population as a function of normalized detuning. Parameters are
$\lambda_{2}=1$ (i.e., $\nu_{e}/\nu_{g}=2$), $\Gamma/\nu_{g}=0.1$,
$\gamma/\Gamma=0.1$ and $\eta_{\ell}/\gamma=2$. (c) Comparison of analytical
results versus numerical simulations for the absorption (shaded in orange) and
emission (shaded in green) profiles. Parameters are $\eta_{\ell}/\gamma=0.1$
and $\Delta_{\ell}/\nu_{g}=0.1$.

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 $\hbox to0.0pt{\hskip
0.48611pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$f$}}(s)=\int_{0}^{\infty}dt\,f(t)\exp(-st)$
for a time-dependent function $f(t)$ at $t\geq 0$. In such a case, Eq. (16)
takes a much simpler form

$\hbox to0.0pt{\hskip
4.24596pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\braket{\sigma}$}}=\frac{\eta_{\ell}}{s}\hbox
to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathcal{G}$}}_{\mathrm{ab}}-\eta_{\ell}\hbox
to0.0pt{\hskip
6.09543pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\braket{\mathcal{P}_{e}}$}}(\hbox
to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color${\mathcal{G}}$}}_{\mathrm{em}}+\hbox
to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color${\mathcal{G}}$}}_{\mathrm{ab}}),$
(18)

with the following functions identified corresponding to emission and
absorption events, respectively

$\displaystyle\hbox to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathcal{G}$}}_{\mathrm{em}}=$
$\displaystyle\sum_{m=0}^{\infty}\frac{S_{m}^{\mathrm{em}}}{s+m\Gamma+\gamma+i\left(\Delta_{\ell}-m\nu_{g}\right)},$
(19) $\displaystyle\hbox to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathcal{G}$}}_{\mathrm{ab}}=$
$\displaystyle\sum_{m=0}^{\infty}\frac{S_{m}^{\mathrm{ab}}}{s+\gamma+m\Gamma+i\left(\Delta_{\ell}+m\nu_{e}\right)}.$
(20)

From these expressions, one can proceed in evaluating analytically the
population of the excited state
$p_{e}^{\mathrm{ss}}=\lim\limits_{t\to\infty}\braket{\sigma^{\dagger}(t)\sigma(t)}$
in steady state (as detailed in Appendix D)

$p_{e}^{\mathrm{ss}}=\frac{\sum\limits_{m=0}^{\infty}\gamma_{m}^{\uparrow}(\omega_{\ell})}{\gamma+\sum\limits_{m=0}^{\infty}\left[\gamma_{m}^{\uparrow}(\omega_{\ell})+\gamma_{m}^{\downarrow}(\omega_{\ell})\right]},$
(21)

The coefficients $\gamma_{m}^{\uparrow}(\omega)$ and
$\gamma_{m}^{\downarrow}(\omega)$ represent the dynamic equilibrium population
transfer rates for absorption from the ground state to the excited state
$\ket{g;0_{g}}\rightarrow\ket{e;m_{e}}$ and emission from the excited to the
ground state $\ket{e;0_{e}}\rightarrow\ket{g;m_{g}}$ as illustrated in Fig.
3(a). The rates are analytically expressed as

$\displaystyle\gamma^{\uparrow}_{m}(\omega)=$
$\displaystyle\frac{\eta_{\ell}^{2}S_{m}^{\mathrm{ab}}\left(m\Gamma+\gamma\right)}{(m\Gamma+\gamma)^{2}+(\omega_{00}+m\nu_{e}-\omega)^{2}},$
(22a) $\displaystyle\gamma^{\downarrow}_{m}(\omega)=$
$\displaystyle\frac{\eta_{\ell}^{2}S_{m}^{\mathrm{em}}\left(m\Gamma+\gamma\right)}{(m\Gamma+\gamma)^{2}+(\omega_{00}-m\nu_{g}-\omega)^{2}}.$
(22b)

Specifically, these rates contribute to the rate equation for the population
of the excited state, given by (see Appendix E for detailed derivations):

$\partial_{t}p_{e}=-2(\gamma+\sum\limits_{m=0}^{\infty}\gamma_{m}^{\downarrow})p_{e}+2\sum_{m=0}^{\infty}\gamma_{m}^{\uparrow}(1-p_{e}).$
(23)

This equation holds true under the condition $\eta_{\ell}\ll\Gamma$.
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 Chenu _et al._
(2019); Banerjee and Gangopadhyay (2002).

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 $\omega_{\ell}$, transitioning to the excited state $\ket{e;m_{e}}$
under the resonant condition $\omega_{\ell}=\omega_{00}+m\nu_{e}$.
Subsequently, after undergoing fast vibrational relaxation, the molecule emits
a photon centered around the frequency $\omega_{00}-m^{\prime}\nu_{g}$, which
can be detected with a modified linewidth $\gamma+m^{\prime}\Gamma$. The
absorption and emission profiles are then obtained by summing up the
contributions from all possible cases, resulting in Lorentzian profiles
represented by $\gamma_{m}^{\uparrow/\downarrow}$, as shown in Fig. 3(c):

$S_{\mathrm{Ab}}=\sum_{m=0}^{\infty}\gamma^{\uparrow}_{m}\quad\mathrm{and}\quad
S_{\mathrm{Em}}=\sum_{m=0}^{\infty}\gamma^{\downarrow}_{m}.$ (24)

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 electron-vibron
coupling under realistic conditions, is expected to only slightly break the
symmetry between the emission and absorption spectra, as the expected values
for $\lambda_{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 Gardiner and
Zoller (2004); Carmichael (1999).

Figure 4: Schematics of an ensemble of molecules inside a Fabry-Pérot
resonator. Cavity photon loss occurs at rates $\kappa_{1}$ and $\kappa_{2}$
via the mirrors M1 and M2, respectively. Light-molecule interactions occur at
rate $g$ while spontaneous emission and cavity driving are at rates $\gamma$
and $\eta_{c}$, respectively.

##  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 Wu _et al._ (2016); del Pino _et al._
(2018a, b); Strashko _et al._ (2018); Sentef _et al._ (2018); Herrera and
Owrutsky (2020); Neuman and Aizpurua (2018); Zeb _et al._ (2018); Kansanen _et
al._ (2019); Du _et al._ (2018); Ribeiro _et al._ (2018) 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. Reitz _et al._ (2019) 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 $\mathcal{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
$\mathcal{H}$ from Eq. (3). The interaction between the $\mathcal{N}$
molecules and the cavity field mode is characterized by the Tavis-Cummings
model,

$\displaystyle\mathcal{H}_{\mathrm{cav}}=$
$\displaystyle\omega_{c}a^{\dagger}a+g\sum_{n=1}^{\mathcal{N}}\left(a\sigma^{\dagger}_{n}+\mathrm{h.c.}\right)$
(25) $\displaystyle+i\eta_{c}(a^{\dagger}e^{-i\omega_{\ell}t}-\mathrm{h.c.}),$

consisting of the free cavity field at frequency $\omega_{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 $\eta_{c}$ and
frequency $\omega_{\ell}$.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
$\tilde{\sigma}_{e,n}$ and $\tilde{\sigma}_{g,n}$ for the $n$th molecule in
the rotating frame at the laser frequency $\omega_{\ell}$:

$\displaystyle\dot{a}=-(i\Delta_{c}+\kappa)a-ig\sum_{n=1}^{\mathcal{N}}\sigma_{n}+\sqrt{2\kappa_{1}}A_{\mathrm{1,in}}+\sqrt{2\kappa_{2}}a_{\mathrm{2,in}},$
(26a) $\displaystyle\begin{aligned}
\dot{\tilde{\sigma}}_{e,n}\approx&-\left(i\Delta_{\ell}+\gamma\right)\tilde{\sigma}_{e,n}+iga{\mathcal{S}}_{e,n}^{\dagger}{\mathcal{D}}_{e,n}^{\dagger}e^{i\left(\nu_{e}-\nu_{g}\right)b^{\dagger}_{e,n}b_{e,n}t}\\\
&+\sqrt{2\gamma}\sigma_{\mathrm{in,n}}{\mathcal{S}}_{e,n}^{\dagger}{\mathcal{D}}_{e,n}^{\dagger}e^{i\left(\nu_{e}-\nu_{g}\right)b^{\dagger}_{e,n}b_{e,n}t},\end{aligned}$
(26b) $\displaystyle\begin{aligned}
\dot{\tilde{\sigma}}_{g,n}\approx&-(i\Delta_{\ell}+\gamma)\tilde{\sigma}_{g,n}-igae^{i\left(\nu_{e}-\nu_{g}\right)b^{\dagger}_{g,n}b_{g,n}t}\mathcal{S}_{g,n}^{\dagger}\mathcal{D}_{g,n}^{\dagger}\\\
&+\sqrt{2\gamma}e^{i\left(\nu_{e}-\nu_{g}\right)b^{\dagger}_{g,n}b_{g,n}t}\mathcal{S}_{g,n}^{\dagger}\mathcal{D}_{g,n}^{\dagger}\sigma_{\mathrm{in,n}}^{\phantom{\dagger}}.\end{aligned}$
(26c)

Here, the total dissipation for the cavity field
$\kappa=\kappa_{1}+\kappa_{2}$ encompasses the losses via both mirrors. The
operators $A_{\mathrm{1,in}}=\eta_{c}/\sqrt{2\kappa_{1}}+a_{\mathrm{1,in}}$
describes the input classical field coming through the left mirror
$\eta_{c}/\sqrt{2\kappa_{1}}$ and the zero-average input noise with the only
non-vanishing two-time correlations
$\braket{a_{\mathrm{1,in}}(t)a_{\mathrm{1,in}}^{\dagger}(t^{\prime})}=\delta(t-t^{\prime})$.
Additionally, zero-average input noise comes through the right side mirror
with similar correlations
$\braket{a^{\phantom{\dagger}}_{\mathrm{2,in}}(t)a_{\mathrm{2,in}}^{\dagger}(t^{\prime})}=\delta(t-t^{\prime})$
and uncorrelated with the $a_{\mathrm{1,in}}(t)$.

The Markovian limit is achieved under the large relaxation rate condition for
vibrational mode, i.e., $\Gamma\gg\kappa$ and $\Gamma\gg\gamma$. 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

$\displaystyle\hbox to0.0pt{\hskip
5.9266pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\left\langle\sigma_{n}\right\rangle$}}=$
$\displaystyle ig\left(\hbox to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathcal{G}$}}_{\mathrm{em}}+\hbox
to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathcal{G}$}}_{\mathrm{ab}}\right)\hbox
to0.0pt{\hskip
11.19681pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\left\langle\mathcal{P}_{e,n}a\right\rangle$}}-ig\hbox
to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathcal{G}$}}_{\mathrm{ab}}\hbox
to0.0pt{\hskip
4.03183pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\left\langle
a\right\rangle$}}.$ (27)

Figure 5: Cavity transmission ($|\mathcal{T}|^{2}$) of the molecule with
various linear and quadratic electron-vibron couplings $\lambda_{1}$ and
$\lambda_{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 $\mathcal{N}=25$ and $100$, respectively. Parameters:
$\omega_{c}=\omega_{00},\,g=3\kappa,\,\nu_{g}=10\kappa,\,\gamma=0.01\kappa,\Gamma=20\kappa,\,2\kappa_{1}=2\kappa_{2}=\kappa$,
and the driving is assumed very weak $\eta_{c}/\kappa=0.001$.

The coupling between the cavity mode $a$ and the projection operator
$\mathcal{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
$\ket{e}$ is negligible (under the condition that $\eta_{c}\ll\kappa$). 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 $\sigma_{n}$ for the $n$-th molecule will be equivalent to
that of the other molecules, i.e.,
$\braket{\sigma}=\braket{\sigma_{n}}=\braket{\sigma_{m}}$ ($m\neq n$). Then,
the equations of motion are written in the vector form (in the Laplace
transform domain) as $\hbox to0.0pt{\hskip
2.08333pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathbf{M}$}}\hbox
to0.0pt{\hskip
0.13889pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathbf{v}$}}+\hbox
to0.0pt{\hskip
0.13889pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathbf{v}$}}_{c}=0$,
with the drift matrix

$\hbox to0.0pt{\hskip
2.08333pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathbf{M}$}}=\left(\begin{array}[]{ccc}-\left(i\Delta_{c}+\kappa\right)-s&-i\mathcal{N}g\\\
-ig&-1/\overline{\mathcal{G}}_{\mathrm{ab}}\end{array}\right),$ (28)

and the definitions $\mathbf{v}=\left(\braket{a},\braket{\sigma}\right)^{T}$
and $\mathbf{v}_{c}=\left(\eta_{c},0\right)^{T}$. The diagonalization of the
drift matrix (under resonance condition $\Delta_{c}=0$) yields the frequencies
$\omega_{\pm}$ and linewidths $\gamma_{\pm}$ Plankensteiner _et al._ (2019);
Herrera and Owrutsky (2020) of the two polaritons as

$\displaystyle\omega_{\pm}=$
$\displaystyle\frac{-\Delta_{\mathrm{eff}}}{2}\pm\frac{1}{2}\mathcal{I}\sqrt{(\Gamma_{\mathrm{eff}}-\kappa+i\Delta_{\mathrm{eff}})^{2}-\mathcal{N}g^{2}},$
(29a) $\displaystyle\gamma_{\pm}=$
$\displaystyle\frac{\Gamma_{\mathrm{eff}}+\kappa}{2}\pm\frac{1}{2}\mathcal{R}\sqrt{(\Gamma_{\mathrm{eff}}-\kappa+i\Delta_{\mathrm{eff}})^{2}-\mathcal{N}g^{2}},$
(29b)

with $\Gamma_{\mathrm{eff}}=\mathcal{R}\lim\limits_{s\rightarrow
0}1/\overline{\mathcal{G}}_{\mathrm{ab}}$ and
$\Delta_{\mathrm{eff}}=\mathcal{I}\lim\limits_{s\rightarrow
0}1/\overline{\mathcal{G}}_{\mathrm{ab}}$ denoting the effective decay rate
and additional frequency shift. 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 light-matter 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

$\mathcal{T}=\frac{\sqrt{2\kappa_{2}}\braket{a}_{\mathrm{ss}}}{\eta_{c}/\sqrt{2\kappa_{1}}},$
(30)

and illustrate its behavior with respect to the scanning laser frequency in
Fig. 5. The quantity $\braket{a}_{\mathrm{ss}}$ is the average value of the
cavity mode amplitude in steady state in the linear response regime

$\displaystyle\braket{a}_{\mathrm{ss}}=$
$\displaystyle\frac{\eta_{c}}{\mathcal{N}g^{2}\chi_{\mathrm{ab}}+\kappa+i(\omega_{c}-\omega_{\ell})},$
(31)

with $\chi_{\mathrm{ab}}=\lim\limits_{s\rightarrow
0}1/\overline{\mathcal{G}}_{\mathrm{ab}}.$

We illustrate numerical and analytical results in Fig. 5 where the profile of
the cavity transmission at $\omega_{c}=\omega_{00}$ is plotted. The presence
of the linear electron-vibron coupling scaling with $\lambda_{1}$ induces an
interaction between upper and lower polaritons already presented at the
theoretical level in a few treatments Neuman and Aizpurua (2018); Reitz _et
al._ (2019, 2020). 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. Reitz _et al._ (2019)
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.

## Acknowledgments

We acknowledge financial support from the Max Planck Society and from the
Deutsche Forschungsgemein- schaft (DFG, German Research Foundation) – Project-
ID 9429529648 – TRR 306 QuCoLiMa (”Quantum Cooperativity of Light and
Matter”).

## References

  * Neuman and Aizpurua (2018) T. Neuman and J. Aizpurua, “Origin of the asymmetric light emission from molecular exciton polaritons,” Optica 5, 1247 (2018).
  * Clark and Dines (1986) H. J. Clark and T. J. Dines, “Resonance Raman spectroscopy, and its application to inorganic chemistry,” Angew. Chem. Int. Ed. 25, 131 (1986).
  * Smith and Dent (2005) E. Smith and G. Dent, “Chapter 4: Resonance Raman Scattering,” in _Modern Raman Spectroscopy: A Practical Approach_ (John Wiley & Sons, Chichester, England, 2005).
  * Riede _et al._ (2018) M. Riede, B. Lüssem, K. Leo, and A. Z. M. S. Rahman, “Organic semiconductors,” in _Reference Module in Materials Science and Materials Engineering_ (Elsevier, 2018).
  * Holstein (1959) T. Holstein, “Studies of polaron motion: Part I. The molecular-crystal model,” Ann. Phys. 8, 325 (1959).
  * Spano (2006) F. C. Spano, “Excitons in conjugated oligomer aggregates, films, and crystals,” Ann. Rev. Phys. Chem. 57, 217 (2006).
  * Reitz _et al._ (2019) M. Reitz, C. Sommer, and C. Genes, “Langevin approach to quantum optics with molecules,” Phys. Rev. Lett. 122, 203602 (2019).
  * Reitz _et al._ (2020) M. Reitz, C. Sommer, B. Gurlek, V. Sandoghdar, D. Martin-Cano, and C. Genes, “Molecule-photon interactions in phononic environments,” Phys. Rev. Res. 2, 033270 (2020).
  * Zhang _et al._ (2023) Z. Zhang, X. Nie, D. Lei, and S. Mukamel, “Multidimensional coherent spectroscopy of molecular polaritons: Langevin approach,” Phys. Rev. Lett. 130, 103001 (2023).
  * Kansanen _et al._ (2021) K. S. U. Kansanen, J. J. Toppari, and T. T. Heikkilä, “Polariton response in the presence of Brownian dissipation from molecular vibrations,” J. Chem. Phys. 154, 044108 (2021).
  * Kansanen _et al._ (2019) K. S. U. Kansanen, A. Asikainen, J. J. Toppari, G. Groenhof, and T. T. Heikkilä, “Theory for the stationary polariton response in the presence of vibrations,” Phys. Rev. B 100, 245426 (2019).
  * Keeling and Kéna-Cohen (2020) J. Keeling and S. Kéna-Cohen, “Bose–Einstein condensation of exciton-polaritons in organic microcavities,” Ann. Rev. Phys. Chem. 71, 435 (2020).
  * Gilmore and McKenzie (2005) J. Gilmore and R. H. McKenzie, “Spin boson models for quantum decoherence of electronic excitations of biomolecules and quantum dots in a solvent,” J. Phys. Condens. Matter 17, 1735 (2005).
  * Aspelmeyer _et al._ (2014a) M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, eds., _Cavity Optomechanics: Nano- and Micromechanical Resonators Interacting with Light_ , 1st ed., Quantum Science and Technology (Springer Berlin, Heidelberg, Berlin, Heidelberg, 2014).
  * Aspelmeyer _et al._ (2014b) M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, “Cavity optomechanics,” Rev. Mod. Phys. 86, 1391 (2014b).
  * Rabl (2011) P. Rabl, “Photon blockade effect in optomechanical systems,” Phys. Rev. Lett. 107, 063601 (2011).
  * Nunnenkamp _et al._ (2011) A. Nunnenkamp, K. Børkje, and S. M. Girvin, “Single-photon optomechanics,” Phys. Rev. Lett. 107, 063602 (2011).
  * Jang (2022) S. J. Jang, “Partially polaron-transformed quantum master equation for exciton and charge transport dynamics,” J. Chem. Phys. 157, 104107 (2022).
  * Liao and Nori (2013) J.-Q. Liao and F. Nori, “Photon blockade in quadratically coupled optomechanical systems,” Phys. Rev. A 88, 023853 (2013).
  * Mustroph _et al._ (2009) H. Mustroph, J. Mistol, B. Senns, D. Keil, M. Findeisen, and L. Hennig, “Relationship between the molecular structure of merocyanine dyes and the vibrational fine structure of their electronic absorption spectra,” Angew. Chem. Int. Ed. 48, 8773 (2009).
  * Miao _et al._ (2023) R. Miao, Y. Fu, D. Lu, F. Liang, H. Yu, H. Zhang, and Y. Wu, “Deciphering the vibronic lasing performances in an electron-phonon-photon coupling system,” Opt. Express 31, 9790 (2023).
  * Zirkelbach _et al._ (2022) J. Zirkelbach, M. Mirzaei, I. Deperasińska, B. Kozankiewicz, B. Gurlek, A. Shkarin, and et al., “High-resolution vibronic spectroscopy of a single molecule embedded in a crystal,” J. Chem. Phys. 156, 104301 (2022).
  * Gardiner and Zoller (2004) C. Gardiner and P. Zoller, _Quantum noise: a handbook of Markovian and non-Markovian quantum stochastic methods with applications to quantum optics_ , Vol. 56 (Springer Science & Business Media, 2004).
  * del Valle and Catalán (2019) J. C. del Valle and J. Catalán, “Kasha’s rule: a reappraisal,” Phys. Chem. Chem. Phys. 21, 10061 (2019).
  * Chenu _et al._ (2019) A. Chenu, S.-Y. Shiau, and M. Combescot, “Two-level system coupled to phonons: Full analytical solution,” Phys. Rev. B 99, 014302 (2019).
  * Banerjee and Gangopadhyay (2002) S. Banerjee and G. Gangopadhyay, “Spectra of displaced distorted oscillator molecular system,” Chem. Phys. Lett. 359, 295 (2002).
  * Carmichael (1999) H. J. Carmichael, _Statistical Methods in Quantum Optics 1: Master Equations and Fokker-Planck Equations_ , 1st ed., Theoretical and Mathematical Physics (Springer Berlin, Heidelberg, Berlin, Heidelberg, 1999).
  * Wu _et al._ (2016) N. Wu, J. Feist, and F. J. Garcia-Vidal, “When polarons meet polaritons: Exciton-vibration interactions in organic molecules strongly coupled to confined light fields,” Phys. Rev. B 94, 195409 (2016).
  * del Pino _et al._ (2018a) J. del Pino, F. A. Y. N. Schröder, A. W. Chin, J. Feist, and F. J. Garcia-Vidal, “Tensor network simulation of non-markovian dynamics in organic polaritons,” Phys. Rev. Lett. 121, 227401 (2018a).
  * del Pino _et al._ (2018b) J. del Pino, F. A. Y. N. Schröder, A. W. Chin, J. Feist, and F. J. Garcia-Vidal, “Tensor network simulation of non-Markovian dynamics in organic polaritons,” Phys. Rev. Lett. 121, 227401 (2018b).
  * Strashko _et al._ (2018) A. Strashko, P. Kirton, and J. Keeling, “Organic polariton lasing and the weak to strong coupling crossover,” Phys. Rev. Lett. 121, 193601 (2018).
  * Sentef _et al._ (2018) M. A. Sentef, M. Ruggenthaler, and A. Rubio, “Cavity quantum-electrodynamical polaritonically enhanced electron-phonon coupling and its influence on superconductivity,” Sci. Adv. 4, 6969 (2018).
  * Herrera and Owrutsky (2020) F. Herrera and J. Owrutsky, “Molecular polaritons for controlling chemistry with quantum optics,” J. Chem. Phys. 152, 100902 (2020).
  * Zeb _et al._ (2018) M. A. Zeb, P. G. Kirton, and J. Keeling, “Exact states and spectra of vibrationally dressed polaritons,” ACS Photonics 5, 249 (2018).
  * Du _et al._ (2018) M. Du, L. A. Martinez-Martinez, R. F. Ribeiro, Z. Hu, V. M. Menon, and J. Yuen-Zhou, “Theory for polariton-assisted remote energy transfer,” Chem. Sci. 9, 6659 (2018).
  * Ribeiro _et al._ (2018) R. F. Ribeiro, L. A. Martínez-Martínez, M. Du, J. Campos-Gonzalez-Angulo, and J. Yuen-Zhou, “Polariton chemistry: controlling molecular dynamics with optical cavities,” Chem. Sci. 9, 6325 (2018).
  * Plankensteiner _et al._ (2019) D. Plankensteiner, C. Sommer, M. Reitz, H. Ritsch, and C. Genes, “Enhanced collective Purcell effect of coupled quantum emitter systems,” Phys. Rev. A 99, 043843 (2019).
  * Johansson _et al._ (2012) J. Johansson, P. Nation, and F. Nori, “Qutip: An open-source python framework for the dynamics of open quantum systems,” Comput. Phys. Commun. 183, 1760 (2012).

##  Appendix A The modified Holstein Hamiltonian

Let us illustrate how the modified Holstein Hamiltonian arises, how it can be
diagonalized and how the quantum Langevin equations for the squeezed and
displaced polaron operators can be derived.

### First principle derivation of the non-linear Hamiltonian

We consider a single molecule with ground $|g\rangle$ (frequency $\omega_{g}$)
and excited $|e\rangle$ (frequency $\omega_{e}$) electronic levels coupled to
the ground phonons ($\nu_{g}$) and the excited phonons ($\nu_{e}$) of a single
vibrational mode with mass $\mu$ respectively. Assuming that the ground and
excited electronic states have different parabolic shape around the minima.
Then the total Hamiltonian of the electron-phonon system reads as ($\hbar=1$).

$\displaystyle\mathcal{H}$ $\displaystyle=$
$\displaystyle\left[\omega_{e}+\dfrac{{\hat{P}}^{2}}{2\mu}+\dfrac{1}{2}\mu\nu^{2}_{e}(\hat{R}-R_{\mathrm{eg}})^{2}\right]\sigma^{\dagger}\sigma+\left(\omega_{g}+\dfrac{{\hat{P}}^{2}}{2\mu}+\dfrac{1}{2}\mu\nu^{2}_{g}{\hat{R}}^{2}\right)\sigma\sigma^{\dagger}$
(32) $\displaystyle=$
$\displaystyle\omega_{g}\sigma\sigma^{\dagger}+\left(\omega_{e}+\mu\nu^{2}_{e}R_{\mathrm{eg}}^{2}/2\right)\sigma^{\dagger}\sigma+\dfrac{{\hat{P}}^{2}}{2\mu}+\dfrac{1}{2}\mu\nu^{2}_{g}{\hat{R}}^{2}+\dfrac{1}{2}\mu\nu^{2}_{e}R_{\mathrm{eg}}\hat{R}\sigma^{\dagger}\sigma+\dfrac{1}{2}\mu(\nu^{2}_{e}-\nu^{2}_{g}){\hat{R}}^{2}\sigma^{\dagger}\sigma,$

where $\sigma=|g\rangle\langle e|$ is the Pauli lowering operator. By
rewriting position ($\hat{R}$) and momentum ($\hat{P}$) in terms of creation
$b^{\dagger}=(\hat{R}/{R_{\mathrm{zpm}}}-iR_{\mathrm{zpm}}\hat{P})/\sqrt{2}$
and annihilation
$b=(\hat{R}/{R_{\mathrm{zpm}}}+iR_{\mathrm{zpm}}\hat{P})/\sqrt{2}$ operators
that fulfill $\left[b,b^{\dagger}\right]=1$, the Hamiltonian in Eq. (32) can
be written as

$\displaystyle\mathcal{H}=\omega_{0}\sigma^{\dagger}\sigma+\nu_{g}b^{\dagger}b+\lambda_{1}\nu_{g}(b+b^{\dagger})\sigma^{\dagger}\sigma+\lambda_{2}\nu_{g}(b+b^{\dagger})^{2}\sigma^{\dagger}\sigma,$
(33)

where
${\omega}_{0}=\omega_{e}-\omega_{g}+\mu\nu^{2}_{e}R_{\text{eg}}^{2}/2=\omega_{e}-\omega_{g}+\lambda_{1}^{2}\nu_{g}^{3}/\nu_{e}^{2}$
is the modified frequency of electronic excited state.
$\lambda_{1}=\mu\nu_{e}^{2}R_{\mathrm{eg}}R_{\mathrm{zpm}}/\nu_{g}$ and
$\lambda_{2}=(\nu_{e}^{2}-\nu_{g}^{2})/4\nu^{2}_{g}$ are the linear and
quadratic coupling constants.

### Quadratic Holstein Hamiltonian Diagonalization

In the presence of both linear and quadratic couplings, the diagonalization of
Hamiltonian in Eq. (3) can be achieved by performing a sequence of unitary
transformations. This transformation could be accomplished by first removing
all the linear terms via the polaron transformation
$\mathcal{U}_{d}={[\mathcal{D}(r_{d})}]^{\sigma^{\dagger}\sigma}=\sigma\sigma^{\dagger}+\mathcal{D}(r_{d})\sigma^{\dagger}\sigma$,
where displacement operator is define as
${\mathcal{D}(r_{d})}=\exp[r_{d}(b^{\dagger}-b)]$. The polaron transformation
has the effect that
$b\rightarrow\mathcal{U}_{d}^{\dagger}b\mathcal{U}_{d}=b+r_{d}\sigma^{\dagger}\sigma$.
Specifically, when

$r_{d}=-\frac{\lambda_{1}}{1+4\lambda_{2}}=-\lambda_{1}\frac{\nu_{g}^{2}}{\nu_{e}^{2}},$
(34)

the resulting Hamiltonian
$\mathcal{H}_{1}=\mathcal{U}_{d}^{\dagger}\mathcal{H}\mathcal{U}_{d}$ can be
written as

$\displaystyle\mathcal{H}_{1}=$
$\displaystyle\left[{\omega}_{0}+r^{2}_{d}\nu_{g}(1+4\lambda_{2})+2\lambda_{1}r_{d}\right]\sigma^{\dagger}\sigma+\nu_{g}b^{\dagger}b+\lambda_{2}\nu_{g}(b+b^{\dagger})^{2}\sigma^{\dagger}\sigma.$
(35)

We now diagonalize this Hamiltonian via squeezing transformation
$\mathcal{U}_{s}={[\mathcal{S}(r_{s})]}^{\sigma^{\dagger}\sigma}=\sigma\sigma^{\dagger}+\mathcal{S}(r_{s})\sigma^{\dagger}\sigma$,
where ${\mathcal{S}(r_{s})}=\exp[r_{s}(b^{2}-b^{\dagger^{2}})/2]$ is a single
mode squeezing operator, so as to remove the quadratic terms. Under this
transformation
$b\rightarrow\mathcal{U}^{\dagger}_{s}b\mathcal{U}_{s}=b\sigma\sigma^{\dagger}+\left[b\cosh(r_{s})+b^{\dagger}\sinh(r_{s})\right]\sigma^{\dagger}\sigma$,
the resulting Hamiltonian
$\mathcal{H}_{2}=\mathcal{U}^{\dagger}_{s}\mathcal{H}_{1}\mathcal{U}_{s}$ can
be written as

$\displaystyle\mathcal{H}_{2}$ $\displaystyle=$
$\displaystyle\nu_{g}b^{\dagger}b\sigma\sigma^{\dagger}+\nu_{e}b^{\dagger}b\sigma^{\dagger}\sigma+{\omega}_{00}\sigma^{\dagger}\sigma$
(36)

under the condition of

$e^{4r_{s}}=1+4\lambda_{2}=\frac{\nu_{e}^{2}}{\nu_{g}^{2}}.$ (37)

Here $\omega_{00}=\omega_{0}+r_{d}\lambda_{1}\nu_{g}+(\nu_{e}-\nu_{g})/2$ is
the zero-phonon line. This transformation could be also accomplished in a
reverse order, by first removing quadratic terms under applying
$\mathcal{U}_{s}^{\prime}=[\mathcal{S}(r_{s})]^{\sigma^{\dagger}\sigma}$ and
then removing the all linear terms via polaron transformation
$\mathcal{U}_{d}^{\prime}=[\mathcal{D}(r_{d}e^{r_{s}})]^{\sigma^{\dagger}\sigma}$.
The results are, of course, identical.

##  Appendix B Derivation of the Effective Quantum Langevin Equations for the
Electronic Transition

### Effective Quantum Langevin Equation for the Vibrational Mode.

Let us consider one special case where the relaxation rate for the vibrations
is much larger than the rate of change for the population on the electric
excited state. Then the dynamical behaviors for vibrations with molecules
populating on the state $\ket{e}$ prefer to being different with that with
molecules occupying in the state $\ket{g}$. This fact motivates the
partitioning of the total Hilbert space into the orthogonal subspaces via the
following two projection operators

$\mathcal{P}_{e}=\sigma^{\dagger}\sigma\quad\text{and}\quad\mathcal{P}_{g}=\sigma\sigma^{\dagger}.$
(38)

The bosonic annihilation operator can thus be partitioned into
$b=b\sigma^{\dagger}\sigma+b\sigma\sigma^{\dagger}$, which gives two dynamical
equations corresponding to the operator $b_{1}=b\sigma^{\dagger}\sigma$ and
$b_{2}=b\sigma\sigma^{\dagger}$ as

$\displaystyle\dot{b}_{1}=$
$\displaystyle\,b\frac{d}{dt}(\sigma^{\dagger}\sigma)-\left[i\left(1+2\lambda_{2}\right)\nu_{g}+\Gamma\right]b_{1}-i\lambda_{1}\nu_{g}\sigma^{\dagger}\sigma-i2\lambda_{2}\nu_{g}b_{1}^{\dagger}+\Gamma
r_{d}\mathcal{P}_{e}+\sqrt{2\Gamma}\mathcal{B}_{1}^{\mathrm{in}}\mathcal{P}_{e},$
(39a) $\displaystyle\dot{b}_{2}=$
$\displaystyle-b\frac{d}{dt}(\sigma^{\dagger}\sigma)-\left(i\nu+\Gamma\right)b_{2}+\sqrt{2\Gamma}\mathcal{B}_{2}^{\mathrm{in}}\mathcal{P}_{g},$
(39b)

where $\mathcal{B}_{1}^{\mathrm{in}}$ and $\mathcal{B}_{2}^{\mathrm{in}}$ are
the noise operators.

Assuming a large relaxation rate for the vibrational mode, the evolution of
$\sigma^{\dagger}\sigma$ can be approximately negligible, while the
vibrational mode rapidly relaxes to the steady situation. Then we can get the
effective dynamical equations of $b_{1}$ and $b_{2}$

$\displaystyle\dot{b}_{1}$
$\displaystyle\approx-\left[i\left(1+2\lambda_{2}\right)\nu_{g}+\Gamma\right]b_{1}-i\lambda_{1}\nu_{g}\mathcal{P}_{e}-i2\lambda_{2}\nu_{g}b_{1}^{\dagger}+\Gamma
r_{d}\mathcal{P}_{e}+\sqrt{2\Gamma}\mathcal{B}_{1}^{\mathrm{in}}\mathcal{P}_{e},$
(40a) $\displaystyle\dot{b_{2}}$
$\displaystyle\approx-\left(i\nu+\Gamma\right)b_{2}+\sqrt{2\Gamma}\mathcal{B}_{2}^{\mathrm{in}}\mathcal{P}_{g}.$
(40b)

By introducing $b_{e}=\mathcal{U}b_{1}\mathcal{U}^{\dagger}=\cosh
r_{s}b_{1}+\sinh r_{s}b^{\dagger}_{1}-r_{d}\exp({r_{s}})\mathcal{P}_{e}$, one
can get the effective Langevin equation for $b_{e}$

$\frac{d}{dt}b_{e}=-(i\nu_{e}+\Gamma)b_{e}+\sqrt{2\Gamma}\mathcal{B}_{\mathrm{e}}^{\mathrm{in}}\mathcal{P}_{e}.$
(41)

Under the condition $\Gamma\gg\gamma$, and $\Gamma\gg\eta_{\ell}$, one can
assume the correlation between noise operators
$\mathcal{B}_{e}^{\mathrm{in}{}}$ and $\mathcal{B}_{g}^{\mathrm{in}{}}$ is
negligible as

$\displaystyle\braket{\mathcal{B}_{e}^{\mathrm{in}{}\dagger}(t)\mathcal{B}_{g}^{\mathrm{in}{}}(\tau)}\approx
0\quad\mathrm{and}\quad\braket{\mathcal{B}_{g}^{\mathrm{in}{}\dagger}(t)\mathcal{B}_{e}^{\mathrm{in}{}}(\tau)}\approx
0,$ (42)

and the nonvanishing correlation functions obeying the fluctuation-dissipation
relation

$\braket{\mathcal{B}_{e}^{\mathrm{in}{}\dagger}(t)\mathcal{B}_{e}^{\mathrm{in}{}}(\tau)}=\delta(t-\tau)\quad\mathrm{and}\quad\braket{\mathcal{B}_{g}^{\mathrm{in}{}\dagger}(t)\mathcal{B}_{g}^{\mathrm{in}{}}(\tau)}=\delta(t-\tau).$
(43)

### 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,
$\tilde{\sigma}_{e}^{\prime}=\sigma\mathcal{S}^{\dagger}_{1}\mathcal{D}^{\dagger}_{1}$.
The Langevin equation
$\dot{\tilde{\sigma}}_{e}^{\prime}=\dot{\sigma}\mathcal{S}^{\dagger}_{1}\mathcal{D}^{\dagger}_{1}+\sigma\partial_{t}(\mathcal{S}^{\dagger}_{1}\mathcal{D}^{\dagger}_{1})$
in a rotating frame at driving frequency $\omega_{\ell}$ can be expressed as

$\displaystyle\dot{\tilde{\sigma}}_{e}^{\prime}\approx$
$\displaystyle-\left[i(\omega_{00}-\omega_{\ell})+\gamma\right]\tilde{\sigma}_{e}^{\prime}-\eta_{\ell}\mathcal{S}^{\dagger}_{1}\mathcal{D}^{\dagger}_{1}\sigma^{\dagger}\sigma+\sqrt{2\gamma}\sigma_{\mathrm{in}}\mathcal{S}^{\dagger}_{1}\mathcal{D}^{\dagger}_{1}$
$\displaystyle-i\left(1-\frac{\nu_{g}}{\nu_{e}}\right)\tilde{\sigma}_{e}^{\prime}\left[\nu_{g}b_{1}^{\dagger}b_{1}+\lambda_{1}\nu_{g}(b_{1}+b_{1}^{\dagger})+\lambda_{2}\nu_{g}(b_{1}+b_{1}^{\dagger})^{2}\right],$
(44)

with
$\mathcal{S}_{1}=\exp[r_{s}(b_{1}^{2}-b_{1}^{\dagger^{2}})/2]\mathcal{P}_{e}$
and
$\mathcal{D}_{1}=\mathrm{exp}\left[r_{d}(b^{\dagger}_{1}-b_{1})\right]\mathcal{P}_{e}$.

Considering that the quadratic terms in the second line of the equation above,
i.e.,
$\nu_{g}b_{1}^{\dagger}b_{1}+\lambda_{1}\nu_{g}(b_{1}+b_{1})+\lambda_{2}\nu_{g}(b_{1}+b_{1}^{\dagger})^{2}$,
can be reformed as $\nu_{e}b_{e}^{\dagger}b_{e}$, where the definition of
$b_{e}$ is introduced in Eq. (41), we can thus reform the equation above into

$\dot{\tilde{\tilde{\sigma}}}_{e}\approx-\left[i(\omega_{00}-\omega_{\ell})+\gamma\right]\tilde{\tilde{\sigma}}_{e}-i(\nu_{e}-\nu_{g})\tilde{\tilde{\sigma}}_{e}b^{{}\dagger}_{e}b_{e}-\eta_{\ell}\mathcal{S}^{{}\dagger}_{e}\mathcal{D}^{{}\dagger}_{e}\sigma^{\dagger}\sigma+\sqrt{2\gamma}\sigma_{\mathrm{in}}\mathcal{S}^{{}\dagger}_{e}\mathcal{D}^{{}\dagger}_{e},$
(45)

with
$\tilde{\tilde{\sigma}}_{e}=\sigma\mathcal{S}^{{}\dagger}_{e}\mathcal{D}^{{}\dagger}_{e}$,
$\mathcal{S}_{e}=\exp\left[\left(r_{s}(b^{{}2}-b^{{}\dagger
2})/2\right)\sigma^{\dagger}\sigma\right]\mathcal{P}_{e}$, and
$\mathcal{D}_{e}=\exp\left[r_{d}(b^{{}\dagger}-b)\sigma^{\dagger}\sigma\right]\mathcal{P}_{e}$.

Taking account of its dynamical equation for
$\exp\left[i\left(\nu_{e}-\nu_{g}\right)b^{{}\dagger}_{e}b_{e}t\right]$ given
by

$\displaystyle\frac{d}{dt}e^{i\left(\nu_{e}-\nu_{g}\right)b^{{}\dagger}_{e}b_{e}t}=$
$\displaystyle\,i\left(\nu_{e}-\nu_{g}\right)\int_{0}^{1}e^{i\alpha\left(\nu_{e}-\nu_{g}\right)b^{{}\dagger}_{e}b_{e}t}\frac{d}{dt}\left(b^{{}\dagger}_{e}b_{e}t\right)e^{i\left(1-\alpha\right)i\left(\nu_{e}-\nu_{g}\right)b^{{}\dagger}_{e}b_{e}t}d\alpha$
$\displaystyle=$
$\displaystyle\,i\left(\nu_{e}-\nu_{g}\right)b^{{}\dagger}_{e}b_{e}e^{i\left(\nu_{e}-\nu_{g}\right)b^{{}\dagger}_{e}b_{e}t}+i\left(\nu_{e}-\nu_{g}\right)t\left[-2\Gamma
b^{{}\dagger}_{e}b_{e}+\sqrt{2\Gamma}\frac{e^{i(\nu_{e}-\nu_{g})t}-1}{i(\nu_{e}-\nu_{g})t}b^{{}\dagger}_{e}\mathcal{B}_{\mathrm{e}}^{\mathrm{in}{}}+\mathrm{h.c.}\right],$
(46)

one can recast Eq. (45) into

$\dot{\tilde{{\sigma}}}_{e}\approx-\left[i(\omega_{00}-\omega_{\ell})+\gamma\right]\tilde{{\sigma}}_{e}-\eta_{\ell}\mathcal{S}^{{}\dagger}_{e}\mathcal{D}^{{}\dagger}_{e}e^{i\left(\nu_{e}-\nu_{g}\right)b^{{}\dagger}_{e}b_{e}t}\sigma^{\dagger}\sigma+\sqrt{2\gamma}\sigma_{\mathrm{in}}\mathcal{S}^{{}\dagger}_{e}\mathcal{D}^{{}\dagger}_{e}e^{i\left(\nu_{e}-\nu_{g}\right)b^{{}\dagger}_{e}b_{e}t},$
(47)

with
$\tilde{{\sigma}}_{e}=\tilde{\tilde{\sigma}}_{e}\exp\left[\left(i\nu_{e}-i\nu_{g}\right)b^{{}\dagger}_{e}b_{e}t\right]$.
Here, we have dropped the second term on the right side of Eq. (46) to obtain
the equation above and to receive a sufficient approximation.

Notably, Eq. (45) only conclude the contribution of the vibrations projecting
to the manifold of $\mathcal{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
$\tilde{\tilde{\sigma}}_{g}=\mathcal{S}_{g}^{\dagger}\mathcal{D}_{g}^{\dagger}\sigma$
for the vibrational mode projecting to the manifold of $\mathcal{P}_{g}$ with
$\mathcal{S}_{g}=\exp\left[r_{s}(b_{g}-b_{g}^{\dagger
2})/2\right]\mathcal{P}_{g}$, and
$\mathcal{D}_{g}=\left[r_{d}(b_{g}-b_{g}^{\dagger})\right]\mathcal{P}_{g}$

$\dot{\tilde{\tilde{\sigma}}}_{g}\approx-\left[i(\omega_{00}-\omega_{\ell})+\gamma\right]\tilde{\tilde{\sigma}}_{g}-i(\nu_{e}-\nu_{g})b^{\dagger}_{g}b_{g}\tilde{\tilde{\sigma}}_{g}+\eta_{\ell}\mathcal{S}_{g}^{\dagger}\mathcal{D}_{g}^{\dagger}\sigma\sigma^{\dagger}+\sqrt{2\gamma}\sigma_{\text{in}}\mathcal{S}_{g}^{\dagger}\mathcal{D}_{g}^{\dagger}\sigma\sigma^{\dagger}.$
(48)

Repeat the process to deriving the dynamical equation of
$\tilde{{\sigma}}_{e}$, one can recast the above equation into

$\dot{\tilde{{\sigma}}}_{g}\approx-\left[i(\omega_{00}-\omega_{\ell})+\gamma\right]\tilde{{\sigma}}_{g}+e^{i(\nu_{e}-\nu_{g})b^{\dagger}_{g}b_{g}t}\eta_{\ell}\mathcal{S}_{g}^{\dagger}\mathcal{D}_{g}^{\dagger}\sigma\sigma^{\dagger}+\sqrt{2\gamma}\sigma_{\text{in}}e^{i(\nu_{e}-\nu_{g})b^{\dagger}_{g}b_{g}t}\mathcal{S}_{g}^{\dagger}\mathcal{D}_{g}^{\dagger}\sigma\sigma^{\dagger},$
(49)

with
$\tilde{{\sigma}}_{g}=\exp\left[\left(i\nu_{e}-i\nu_{g}\right)b^{\dagger}_{g}b_{g}t\right]\tilde{\tilde{\sigma}}_{g}$.

Meanwhile, one could also get the dynamical equation for the population of the
excited state given by

$\frac{d}{dt}\sigma^{\dagger}\sigma=-2\gamma\sigma^{\dagger}\sigma+\eta_{\ell}(\sigma+\sigma^{\dagger})+\sqrt{2\gamma}(\sigma^{\dagger}\sigma_{\mathrm{in}}+\sigma\sigma_{\mathrm{in}}^{\dagger}).$
(50)

##  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
$\mathcal{S}_{g}^{\dagger}\mathcal{D}_{g}^{\dagger}$,
$\mathcal{S}_{e}^{\dagger}\mathcal{D}_{e}^{\dagger}$.

### Nonlinear Vibrational Dynamics.

Figure 6: The two-time correlation function in the steady state for (a)
$\braket{b_{g}(\tau)b_{g}^{\dagger}(0)}$, (b)
$\braket{b_{e}(\tau)b_{e}^{\dagger}(0)}$, (c)
$\braket{b_{e}(\tau)b_{g}^{\dagger}(0)}$, and
$\braket{b_{g}(\tau)b_{e}^{\dagger}(0)}$ normalized in a unit factor. The
blue-solid lines in (a-d) are generated with the toolbox QuTip Johansson _et
al._ (2012) under the following parameters $\nu_{g}=1$, $\nu_{e}/\nu_{g}=2$,
$\Gamma/\nu_{g}=0.1$, $\gamma/\nu_{g}=0.01$ and $\eta_{\ell}/\gamma=2$. The
red-dashed lines in (a) and (b) are calculated via the first and the second
term in the right side of Eq. (55), respectively.

In the previous section, we have derived the effective Langevin equation (41)
describing the dynamics of the vibrational mode projected onto the manifold of
$\mathcal{P}_{e}$. One can easily obtain the exact solution of such an
equation, which is given by

$b_{e}=\sqrt{2\Gamma}\int_{-\infty}^{t}d\tau\,e^{-(i\nu_{e}+\Gamma)(t-\tau)}\mathcal{B}_{e}^{\mathrm{in}}(\tau)\mathcal{P}_{e}(\tau).$
(51)

Due to the correlation time for vibrations being much short than the
electrons, the quantity for $\sigma^{\dagger}\sigma(\tau)$ varies little
around $\sigma^{\dagger}\sigma(t)$. We can thus proceed via a Markov
approximation with taking it out of the integral, which yields

$b_{e}=\sqrt{2\Gamma}\mathcal{P}_{e}(t)\int_{-\infty}^{t}d\tau\,e^{-(i\nu_{e}+\Gamma)(t-\tau)}\mathcal{B}_{e}^{\mathrm{in}}(\tau).$
(52)

With this, we can derive the two-time correlation for $t>\tau$

$\displaystyle\braket{b_{e}(t)b_{e}^{\dagger}(\tau)}_{\mathrm{vib}}=$
$\displaystyle e^{-(i\nu_{e}+\Gamma)(t-\tau)}\mathcal{P}_{e}(t),$ (53a)
$\displaystyle\braket{b_{e}(\tau)b_{e}^{\dagger}(t)}_{\mathrm{vib}}=$
$\displaystyle e^{-(-i\nu_{e}+\Gamma)(t-\tau)}\mathcal{P}_{e}(t),$ (53b) where
$\braket{\cdot}_{\mathrm{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 $\mathcal{P}_{g}$ as

$\displaystyle\braket{b_{g}^{\phantom{\dagger}}(t)b_{g}^{\dagger}(\tau)}_{\mathrm{vib}}=$
$\displaystyle
e^{-\left(i\nu_{g}+\Gamma\right)\left(t-\tau\right)}\sigma\sigma^{\dagger}(t),$
(54a)
$\displaystyle\braket{b_{g}^{\phantom{\dagger}}(\tau)b_{g}^{\dagger}(t)}_{\mathrm{vib}}=$
$\displaystyle
e^{\left(i\nu_{g}-\Gamma\right)\left(t-\tau\right)}\sigma\sigma^{\dagger}(t).$
(54b)

Thus, the two-time correlation function for the vibrational mode in the whole
Hilbert space reads

$\displaystyle\braket{b(t+\tau)b^{\dagger}(t)}\approx$
$\displaystyle\braket{b_{g}^{\phantom{\dagger}}(\tau)b_{g}^{\dagger}(0)}+\braket{b_{e}^{\phantom{\dagger}}(\tau)b_{e}^{\dagger}(0)}$
$\displaystyle\approx$
$\displaystyle\,e^{-\left(i\nu_{g}+\Gamma\right)\left(\tau\right)}\braket{\sigma\sigma^{\dagger}(t)}+e^{-\left(i\nu_{e}+\Gamma\right)\left(\tau\right)}\braket{\sigma^{\dagger}\sigma(t)}.$
(55)

Here, the correlations between different manifolds have been properly dropped
because of the much smaller value for such terms
$\braket{b_{e}(\tau)b_{g}^{\dagger}(0)}$ and
$\braket{b_{g}(\tau)b_{e}^{\dagger}(0)}$ comparing with that in the same
manifold, as illustrate in Fig. (6).

### Two-time Correlation Function for the Displacement-Squeezing Operator.

The squeezing operator can be written in a disentangled form

$\mathcal{S}=\exp\left[-\frac{1}{2}\tanh r_{s}b^{\dagger
2}\right]\exp\left[-\ln\cosh
r_{s}\left(b^{\dagger}b+\frac{1}{2}\right)\right]\exp\left[\frac{1}{2}\tanh
r_{s}b^{2}\right],$ (56)

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 $\mathcal{P}_{e}$ is given by

$\displaystyle\left\langle\mathcal{S}_{e}^{{}\dagger}\left(\tau\right)\mathcal{D}_{e}^{{}\dagger}\left(\tau\right)\mathcal{D}_{e}\left(t\right)\mathcal{S}_{e}\left(t\right)\right\rangle_{\mathrm{vib}}=$
$\displaystyle\frac{e^{-r_{d}^{2}\alpha}}{\cosh(r_{s})}\left\langle e^{-\beta
b_{e}^{{}2}\left(\tau\right)}e^{\alpha r_{d}b_{e}\left(\tau\right)}e^{\alpha
r_{d}b_{e}^{{}\dagger}\left(t\right)}e^{-\beta b_{e}^{{}\dagger
2}\left(t\right)}\right\rangle_{\mathrm{vib}},$ (57)

with $\alpha=\tanh r_{s}+1=2\nu_{e}/(\nu_{e}+\nu_{g})$ and $\beta=\left(\tanh
r_{s}\right)/2$.

Using the generating function of the Hermite polynomials given by

$e^{2x\varphi-\varphi^{2}}=\sum_{n=0}^{\infty}H_{n}\left(x\right)\frac{\varphi^{n}}{n!},$
(58)

we can expand the first and second exponential in the equation above in terms
of Hermite polynomials with $x=\alpha r_{d}/(2\sqrt{\beta})$ and
$\varphi=\sqrt{\beta}b_{e}$. Thus, we finally obtain the expression for the
two-time correlation function under the Isserlis’ theorem

$\displaystyle\left\langle\mathcal{S}_{e}^{{}\dagger}\left(\tau\right)\mathcal{D}_{e}^{{}\dagger}\left(\tau\right)\mathcal{D}_{e}\left(t\right)\mathcal{S}_{e}\left(t\right)\right\rangle_{\mathrm{vib}}$
$\displaystyle=\mathcal{F}_{\mathrm{em}}(-\nu_{e},t-\tau)\mathcal{P}_{e}(t)$
(59)

where

$\displaystyle\mathcal{F}_{\mathrm{em}}(\nu,t)=$
$\displaystyle\sum_{m=0}^{\infty}S_{m}^{\mathrm{em}}e^{-m\left(i\nu+\Gamma\right)t},$
(60) $\displaystyle S_{m}^{\mathrm{em}}=$
$\displaystyle\frac{e^{-r_{d}^{2}\alpha}}{\cosh(r_{s})}\left[H_{m}\left(\frac{\alpha
r_{d}}{2\sqrt{\beta}}\right)\right]^{2}\frac{\beta^{m}}{m!}.$

Of course, we can obtain the result for the following equation in the same way

$\displaystyle\left\langle\mathcal{S}_{e}^{{}\dagger}\left(t\right)\mathcal{D}_{e}^{{}\dagger}\left(\tau\right)e^{i\left(\nu_{e}-\nu_{g}\right)b_{e}^{{}\dagger}b_{e}\tau}e^{-i\left(\nu_{e}-\nu_{g}\right)b_{e}^{{}\dagger}b_{e}t}\mathcal{D}_{e}\left(t\right)\mathcal{S}_{e}\left(t\right)\right\rangle_{\mathrm{vib}}=\mathcal{F}_{\mathrm{em}}\left(-\nu_{g},t-\tau\right)\mathcal{P}_{e}(t).$
(61)

The two-time correlation in the manifold of $\mathcal{P}_{g}$ reads

$\displaystyle\braket{\mathcal{D}_{g}(t)\mathcal{S}_{g}(t)\mathcal{S}_{g}^{\dagger}(\tau))\mathcal{D}^{\dagger}_{g}(\tau)}_{\mathrm{vib}}=$
$\displaystyle\mathcal{F}_{\mathrm{ab}}(\nu_{g},t-\tau)\mathcal{P}_{g}(t),$
(62)

and also

$\braket{\mathcal{D}_{g}(t)\mathcal{S}_{g}(t)e^{-i\left(\nu_{e}-\nu_{g}\right)b_{g}^{\dagger}b_{g}t}e^{i\left(\nu_{e}-\nu_{g}\right)b_{g}^{\dagger}b_{g}\tau}\mathcal{S}_{g}^{\dagger}(\tau))\mathcal{D}^{\dagger}_{g}(\tau)}_{\mathrm{vib}}=\mathcal{F}_{\mathrm{ab}}(\nu_{e},t-\tau)\mathcal{P}_{g}(t),$
(63)

where

$\displaystyle\mathcal{F}_{\mathrm{ab}}(\nu,t)=$
$\displaystyle\sum_{m=0}^{\infty}S_{m}^{\mathrm{ab}}e^{-m\left(i\nu+\Gamma\right)t},$
(64a) $\displaystyle S_{m}^{\mathrm{ab}}=$ $\displaystyle\frac{e^{\alpha
r_{d}^{2}\exp\left(2r_{s}\right)}}{\cosh(r_{s})}\left[H_{m}\left(-\frac{i\alpha
r_{d}\exp\left(r_{s}\right)}{2\sqrt{\beta}}\right)\right]^{2}\frac{\left(-\beta\right)^{m}}{m!},$
(64b)

with $\alpha=\tanh r_{s}-1=-2\nu_{g}/(\nu_{e}+\nu_{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. (47) and (49),

$\displaystyle\sigma(t)=$
$\displaystyle\int_{0}^{t}d\tau\,\left(-\eta_{\ell}+\sqrt{2\Gamma}\sigma_{\mathrm{in}}\left(\tau\right)\right)e^{\left[-i\left(\omega_{00}-\omega_{\ell}\right)+\gamma\right]\left(t-\tau\right)}\mathcal{S}_{e}^{{}\dagger}\left(t\right)\mathcal{D}_{e}^{{}\dagger}\left(\tau\right)e^{i\left(\nu_{e}-\nu_{g}\right)b_{e}^{{}\dagger}b_{e}\tau}e^{-i\left(\nu_{e}-\nu_{g}\right)b_{e}^{{}\dagger}b_{e}t}\mathcal{D}_{e}\left(t\right)\mathcal{S}_{e}\left(t\right)$
(65)
$\displaystyle+\int_{0}^{t}d\tau\,\left(\eta_{\ell}+\sqrt{2\Gamma}\sigma_{\mathrm{in}}\left(\tau\right)\right)e^{\left[-i\left(\omega_{00}-\omega_{\ell}\right)+\gamma\right]\left(t-\tau\right)}\mathcal{D}_{g}(t)\mathcal{S}_{g}(t)e^{-i\left(\nu_{e}-\nu_{g}\right)b_{g}^{\dagger}b_{g}t}e^{i\left(\nu_{e}-\nu_{g}\right)b_{g}^{\dagger}b_{g}\tau}\mathcal{S}_{g}^{\dagger}(\tau))\mathcal{D}^{\dagger}_{g}(\tau),$

with the initial value $\sigma(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. (61) as well as Eq. (63) into the equation above, we then
obtain

$\displaystyle\sigma(t)=$
$\displaystyle\int_{0}^{\infty}d\tau\,\left(-\eta_{\ell}+\sqrt{2\Gamma}\sigma_{\mathrm{in}}\left(\tau\right)\right)\mathcal{G}_{\mathrm{em}}\left(-\nu_{g},t-\tau\right)\mathcal{P}_{e}\left(\tau\right)+\int_{0}^{\infty}d\tau\,\left(\eta_{\ell}+\sqrt{2\Gamma}\sigma_{\mathrm{in}}\left(\tau\right)\right)\mathcal{G}_{\mathrm{ab}}\left(\nu_{e},t-\tau\right)\mathcal{P}_{g}\left(\tau\right),$
(66)

with
$\mathcal{G}_{Q}(\nu,\tau)=\exp\left[\left(i\omega_{\ell}-i\omega_{00}+\gamma\right)\left(t-\tau\right)\right]\mathcal{F}_{Q}(\nu,t-\tau)\Theta(t-\tau)$
and $Q\in\left\\{\mathrm{em},\mathrm{ab}\right\\}$, where $\Theta(t)$ is the
Heaviside step function.

Tracing over the electronic transition, we finally obtain the simplified
formal solution for $\braket{\sigma}$ as

$\braket{\sigma(t)}=-\eta_{\ell}\int_{0}^{\infty}d\tau\,\mathcal{G}_{\mathrm{em}}(-\nu_{g},t-\tau)P_{e}(\tau)+\eta_{\ell}\int_{0}^{\infty}d\tau\,\mathcal{G}_{\mathrm{ab}}(\nu_{e},t-\tau)[1-P_{e}(\tau)],$
(67)

with the population on the electronic excited state
$P_{e}=\braket{\sigma^{\dagger}\sigma}$.

We proceed our calculation via the Laplace transformation (defined as $\hbox
to0.0pt{\hskip
0.48611pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$f$}}(s)=\int_{0}^{\infty}dt\,f(t)\exp(-st)$
for a time-dependent function $f(t)$ at $t\geq 0$). Then Eq. (67) can be
written in the Laplace domain as

$\hbox to0.0pt{\hskip
4.24596pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\braket{\sigma}$}}=\frac{\eta_{\ell}}{s}\hbox
to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathcal{G}$}}_{\mathrm{ab}}-\eta_{\ell}\hbox
to0.0pt{\hskip
1.40451pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$P$}}_{e}(\hbox
to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color${\mathcal{G}}$}}_{\mathrm{ab}}+\hbox
to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color${\mathcal{G}}$}}_{\mathrm{em}}),$
(68)

where $\hbox to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathcal{G}$}}_{\mathrm{em}}$
and $\hbox to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathcal{G}$}}_{\mathrm{ab}}$
are the Laplace transform of $\mathcal{G}_{\mathrm{ab}}(\nu_{e},t)$ and
$\mathcal{G}_{\mathrm{em}}(-\nu_{g},t)$, expressed as

$\displaystyle\hbox to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathcal{G}$}}_{\mathrm{ab}}=$
$\displaystyle\sum_{m=0}^{\infty}\frac{S_{m}^{\mathrm{ab}}}{s+m\Gamma+\gamma+i\left(\omega_{00}-\omega_{\ell}+m\nu_{e}\right)},$
(69a) $\displaystyle\hbox to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathcal{G}$}}_{\mathrm{em}}=$
$\displaystyle\sum_{m=0}^{\infty}\frac{S_{m}^{\mathrm{em}}}{s+\gamma+m\Gamma+i\left(\omega_{00}-\omega_{\ell}-m\nu_{g}\right)}.$
(69b)

Assuming that the molecule is prepared in the electronic ground state
$\ket{g}$ followed by taking an average over the electronic transition on both
sides of Eq. (50), and applying the Laplace transformation, we finally obtain

$s\hbox to0.0pt{\hskip
1.40451pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$P$}}_{e}=-2\gamma\hbox
to0.0pt{\hskip
2.70827pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$P_{e}$}}+\eta_{\ell}(\hbox
to0.0pt{\hskip
4.24596pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\braket{\sigma}$}}+\hbox
to0.0pt{\hskip
4.24596pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\braket{\sigma}$}}^{*}).$
(70)

Plugging Eq. (70) into Eq. (67), we can get

$\displaystyle\hbox to0.0pt{\hskip
4.24596pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\braket{\sigma}$}}=$
$\displaystyle\frac{i2\eta_{\ell}^{3}\mathcal{I}[\hbox to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathcal{G}$}}_{\mathrm{ab}}\hbox
to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathcal{G}$}}_{\mathrm{em}}^{*}]+\eta_{\ell}(s+2\gamma)\hbox
to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathcal{G}$}}_{\mathrm{ab}}}{s(2\gamma+s+2\eta_{\ell}^{2}\mathcal{R}[\hbox
to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathcal{G}$}}_{\mathrm{ab}}+\hbox
to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathcal{G}$}}_{\mathrm{em}}])},$
(71a) $\displaystyle\hbox to0.0pt{\hskip
1.40451pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$P$}}_{e}=$
$\displaystyle\frac{2\eta_{\ell}^{2}}{s(2\gamma+s+2\eta_{\ell}^{2}\mathcal{R}[\hbox
to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathcal{G}$}}_{\mathrm{ab}}+\hbox
to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathcal{G}$}}_{\mathrm{em}}])}\mathcal{R}[\hbox
to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathcal{G}$}}_{\mathrm{ab}}],$
(71b)

where $\mathcal{I}[\cdot]$ and $\mathcal{R}[\cdot]$ denote taking the
imaginary and real part, respectively. According to the final value theorem,
we get the steady values

$\displaystyle\braket{\sigma}_{\mathrm{ss}}=$ $\displaystyle\lim_{s\rightarrow
0}s\hbox to0.0pt{\hskip
4.24596pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\braket{\sigma}$}}=\frac{i\eta_{\ell}^{3}\mathcal{I}[{\chi}_{\mathrm{ab}}{\chi}_{\mathrm{em}}^{*}]+\eta_{\ell}\gamma\chi_{\mathrm{ab}}}{\gamma+\eta_{\ell}^{2}\mathcal{R}[{\chi}_{\mathrm{ab}}+{\chi}_{\mathrm{em}}]},$
(72a) $\displaystyle P_{e}^{\mathrm{ss}}=$ $\displaystyle\lim_{s\rightarrow
0}s\hbox to0.0pt{\hskip
1.40451pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$P$}}_{e}=\frac{\eta_{\ell}^{2}}{\gamma+\eta_{\ell}^{2}\mathcal{R}[{\chi}_{\mathrm{ab}}+{\chi}_{\mathrm{em}}]}\mathcal{R}[{\chi}_{\mathrm{ab}}],$
(72b)

with $\chi_{Q}=\lim\limits_{s\to 0}\hbox to0.0pt{\hskip
1.4236pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathcal{G}$}}_{Q}$.
In the limit of weak driving $\eta_{\ell}\ll\gamma$, the equation above can be
simplified to

$\displaystyle\braket{\sigma}_{\mathrm{ss}}\rightarrow$
$\displaystyle\eta_{\ell}\chi_{\mathrm{ab}}=\sum_{m=0}^{\infty}\frac{S_{m}^{\mathrm{ab}}}{m\Gamma+\gamma+i\left(\omega_{00}-\omega_{\ell}+m\nu_{e}\right)},$
(73a) $\displaystyle P_{e}^{\mathrm{ss}}\rightarrow$
$\displaystyle\frac{\eta_{\ell}^{2}}{\gamma}\mathcal{R}[\chi_{\mathrm{ab}}]=\frac{\eta_{\ell}^{2}}{\gamma}\sum_{m=0}^{\infty}\frac{S_{m}^{\mathrm{ab}}\left(m\Gamma+\gamma\right)}{\left(m\Gamma+\gamma\right)^{2}+\left(\omega_{00}-\omega_{\ell}+m\nu_{e}\right)^{2}}.$
(73b)

##  Appendix E Rate Equation

For large vibrational relaxation rates $\Gamma\gg\eta_{\ell}$, the electronic
transition is usually going from the lowest vibrational state in both
electronic states $\ket{g}$ and $\ket{e}$. Thus, the motion of the population
$p^{e}_{m}$ on the state $\ket{e;m_{e}}$ as well as the population $p_{m}^{g}$
on the state $\ket{g;m_{g}}$ are given phenomenology by

$\displaystyle\begin{aligned}
\partial_{t}p_{m}^{e}=&2\gamma_{m}^{\uparrow}p_{0}^{g}+\Gamma
p_{m+1}^{e}-\Gamma p_{m}^{e},\\\
\partial_{t}p_{m-1}^{e}=&2\gamma_{m-1}^{\uparrow}p_{0}^{g}+\Gamma
p_{m}^{e}-\Gamma p_{m-1}^{e},\\\ &\vdots\\\
\partial_{t}p_{0}^{e}=&2\gamma_{0}^{\uparrow}p_{0}^{g}+\Gamma
p_{0}^{e}-2\gamma
p_{0}^{e}-2\sum_{m=0}^{\infty}\gamma_{m}^{\downarrow}p_{0}^{e},\end{aligned}$
(74a) $\displaystyle\begin{aligned}
\partial_{t}p_{m}^{g}=&2\gamma_{m}^{\downarrow}p_{0}^{e}+2\gamma_{m}p_{0}^{e}+\Gamma
p_{m+1}^{g}-\Gamma p_{m}^{g},\\\
\partial_{t}p_{m-1}^{e}=&2\gamma_{m-1}^{\downarrow}p_{0}^{e}+2\gamma_{m-1}p_{0}^{e}+\Gamma
p_{m}^{g}-\Gamma p_{m-1}^{g},\\\ &\vdots\\\
\partial_{t}p_{0}^{g}=2&\gamma_{0}^{\downarrow}p_{e}+2\gamma_{0}p_{0}^{e}+\Gamma
p_{0}^{e}-2\sum_{m=0}^{\infty}\gamma_{m}^{\uparrow}p_{0}^{e},\end{aligned}$
(74b)

where $\gamma_{m}$ describe the incoherent spontaneous emission progress
$\ket{e;0}\rightarrow\ket{g;m}$ satisfying
$\sum\limits_{m=0}^{\infty}\gamma_{m}=\gamma$. It is obvious that the total
population on the state $\ket{e}$ is the sum of all the occupations of its
sublevels as $p_{e}=\sum\limits_{m=0}^{\infty}p_{m}^{e}$. Thus, we have

$\dot{p}_{e}=2\sum\limits_{m=0}^{\infty}\gamma_{m}^{\uparrow}p_{0}^{g}-2\gamma
p_{0}^{e}-2\sum_{m=0}^{\infty}\gamma_{m}^{\downarrow}p_{0}^{e}.$ (75)

Considering less population on states $\ket{e;m_{e}>0}$ and $\ket{g;m_{g}>0}$
due to the large vibrational relaxation $\Gamma$, we can simplify Eq. (75)
into

$\dot{p}_{e}=2\sum\limits_{m=0}^{\infty}\gamma_{m}^{\uparrow}(1-p_{e})-2\left(\gamma+\sum_{m=0}^{\infty}\gamma_{m}^{\downarrow}\right)p_{e}$
(76)

by assuming $p_{e}\approx p_{0}^{e}$ and $p_{g}=1-p_{e}\approx p_{0}^{g}$.

##  Appendix F Absorption and Emission Spectra

### Effective Quantum Langevin Equation for the Commutator.

In this subsection, we will describe the process to find the absorption and
emission spectroscopic signal. Before that, let us introduce the following
correlators “dressed” by a vibrational mode

$\tilde{C}_{\mathrm{em}}^{e}(t+\tau)=\sigma^{\dagger}(t)\sigma(t+\tau)\mathcal{S}^{{}\dagger}_{e}(t+\tau)\mathcal{D}^{{}\dagger}_{e}(t+\tau)\quad\mathrm{and}\quad\tilde{C}_{\mathrm{em}}^{g}(t+\tau)=\mathcal{S}_{g}^{\dagger}(t+\tau)\mathcal{D}_{g}^{\dagger}(t+\tau)\sigma^{\dagger}(t)\sigma(t+\tau)$
(77)

for emission and

$\tilde{C}_{\mathrm{ab}}^{e}(t+\tau)=\sigma(t+\tau)\sigma^{\dagger}(t)\mathcal{S}^{{}\dagger}_{e}(t+\tau)\mathcal{D}^{{}\dagger}_{e}(t+\tau)\quad\mathrm{and}\quad\tilde{C}_{\mathrm{ab}}^{g}(t+\tau)=\mathcal{S}_{g}^{\dagger}(t+\tau)\mathcal{D}_{g}^{\dagger}(t+\tau)\sigma(t+\tau)\sigma^{\dagger}(t)$
(78)

for absorption. According to the discussion in the Subsec. (B), one will
obtain

$\displaystyle\begin{aligned}
\frac{d{\tilde{C}}_{\mathrm{em}}^{e}}{d\tau}\approx&-\left(i\omega_{00}+\gamma\right)\tilde{C}_{\mathrm{em}}^{e}-i(\nu_{e}-\nu_{g})\tilde{C}_{\mathrm{em}}^{e}b^{{}\dagger}_{e}b_{e}(t+\tau)-\eta_{\ell}e^{-i\omega
t}\sigma^{\dagger}(t)\mathcal{S}^{{}\dagger}_{e}(t+\tau)\mathcal{D}^{{}\dagger}_{e}(t+\tau)\mathcal{P}_{e}(t+\tau)\\\
&+\sqrt{2\gamma}\sigma^{\dagger}(t)\sigma_{\mathrm{in}}(t+\tau)\mathcal{S}^{{}\dagger}_{e}(t+\tau)\mathcal{D}^{{}\dagger}_{e}(t+\tau)\mathcal{P}_{e}(t+\tau),\end{aligned}$
(79a) $\displaystyle\begin{aligned}
\frac{d{\tilde{C}}_{\mathrm{em}}^{g}}{d\tau}\approx&-\left(i\omega_{00}+\gamma\right)\tilde{C}_{\mathrm{em}}^{g}-i(\nu_{e}-\nu_{g})b^{\dagger}_{g}b_{g}(t+\tau)\tilde{C}_{\mathrm{em}}^{e}+\eta_{\ell}e^{-i\omega
t}\mathcal{S}_{g}^{\dagger}(t+\tau)\mathcal{D}_{g}^{\dagger}(t+\tau)\mathcal{P}_{g}(t+\tau)\sigma^{\dagger}(t)\\\
&+\sqrt{2\gamma}\mathcal{S}_{g}^{\dagger}(t+\tau)\mathcal{D}_{g}^{\dagger}(t+\tau)\mathcal{P}_{e}(t+\tau)\sigma^{\dagger}(t)\sigma_{\mathrm{in}}(t+\tau),\end{aligned}$
(79b) $\displaystyle\begin{aligned}
\frac{d{\tilde{C}}_{\mathrm{ab}}^{e}}{d\tau}\approx&-\left(i\omega_{00}+\gamma\right)\tilde{C}_{\mathrm{ab}}^{e}-i(\nu_{e}-\nu_{g})\tilde{C}_{\mathrm{ab}}^{e}b^{{}\dagger}_{e}b_{e}(t+\tau)-\eta_{\ell}e^{-i\omega
t}\sigma^{\dagger}(t)\mathcal{S}^{{}\dagger}_{e}(t+\tau)\mathcal{D}^{{}\dagger}_{e}(t+\tau)\mathcal{P}_{e}(t+\tau)\\\
&+\sqrt{2\gamma}\sigma_{\mathrm{in}}(t+\tau)\sigma^{\dagger}(t)\mathcal{S}^{{}\dagger}_{e}(t+\tau)\mathcal{D}^{{}\dagger}_{e}(t+\tau)\mathcal{P}_{e}(t+\tau),\end{aligned}$
(79c) $\displaystyle\begin{aligned}
\frac{d{\tilde{C}}_{\mathrm{ab}}^{g}}{d\tau}\approx&-\left(i\omega_{00}+\gamma\right)\tilde{C}_{\mathrm{ab}}^{g}-i(\nu_{e}-\nu_{g})b^{\dagger}_{g}b_{g}(t+\tau)\tilde{C}_{\mathrm{ab}}^{e}+\eta_{\ell}e^{-i\omega
t}\mathcal{S}_{g}^{\dagger}(t+\tau)\mathcal{D}_{g}^{\dagger}(t+\tau)\mathcal{P}_{g}(t+\tau)\sigma^{\dagger}(t)\\\
&+\sqrt{2\gamma}\mathcal{S}_{g}^{\dagger}(t+\tau)\mathcal{D}_{g}^{\dagger}(t+\tau)\mathcal{P}_{e}(t+\tau)\sigma_{\mathrm{in}}(t+\tau)\sigma^{\dagger}(t).\end{aligned}$
(79d)

### Emission Spectra in the Transient Regime.

We now assume the molecule is initially prepared in the excited state
$\ket{e;\tilde{0}}$ to compute the spectrum of emission in the transient
regime. By setting $\eta_{\ell}=0$, one could obtain the expression of the
two-time correlation function $\braket{\sigma^{\dagger}(0)\sigma(\tau)}$
through Eq. (79)

$\braket{\sigma^{\dagger}(0)\sigma(t)}=\mathcal{F}_{\mathrm{em}}(\nu_{g},\tau)e^{-\left(i\omega_{00}+\gamma\right)t}.$
(80)

Taking the Fourier transformation gives the expression of the emission
spectrum

$S_{\mathrm{Em}}(\omega)=2\mathcal{R}\int_{0}^{\infty}d\tau\,\mathcal{F}_{\mathrm{em}}(\nu_{g},\tau)e^{\left[i\omega-\left(i\omega_{00}+\gamma\right)\right]\tau}=\sum_{m=0}^{\infty}\frac{S_{m}^{\mathrm{em}}\left(\gamma+m\Gamma\right)}{\left(\gamma+m\Gamma\right)^{2}+\left(\omega_{00}-\omega-m\nu_{g}\right)^{2}}.$
(81)

### Absorption Spectra in the Stability Regime.

The formal solution of the two-time correlation function
$\braket{\sigma(t+\tau)\sigma^{\dagger}(t)}$ is given

$\displaystyle\braket{\sigma(t+\tau)\sigma^{\dagger}(t)}=$
$\displaystyle-\eta_{\ell}\int_{t}^{\infty}d\tau\mathcal{G}_{\mathrm{em}}(-\nu_{g},t+\tau-\tau)\braket{\mathcal{P}_{e}\left(\tau\right)\sigma^{\dagger}(t)}+\eta_{\ell}\int_{t}^{\infty}d\tau\mathcal{G}_{\mathrm{ab}}(\nu_{e},t+\tau-\tau)\braket{\mathcal{P}_{g}\left(\tau\right)\sigma^{\dagger}(t)}$
(82)
$\displaystyle+\mathcal{P}_{g}(t)\mathcal{F}_{\mathrm{ab}}(\nu_{e},\tau)e^{-\left(i\omega_{00}+\gamma\right)t}.$

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.,
$\mathcal{P}_{e}(t)\ll 1$ and $\mathcal{P}_{g}(t)\approx 1$. The solution will
become

$\displaystyle\braket{\sigma(t+\tau)\sigma^{\dagger}(t)}\approx$
$\displaystyle\,\eta_{\ell}\int_{t}^{\infty}d\tau\,\mathcal{G}_{\mathrm{ab}}(\nu_{e},t+\tau-\tau)\braket{\sigma^{\dagger}(t)}+\mathcal{F}_{\mathrm{ab}}(\nu_{e},\tau)e^{-\left(i\omega_{00}+\gamma\right)\tau}$
(83) $\displaystyle=$
$\displaystyle\,\eta_{\ell}\int_{t}^{\infty}d\tau\,e^{-(i\omega_{00}-i\omega_{\ell}+\gamma)(t+\tau-\tau)}\mathcal{F}_{\mathrm{ab}}(\nu_{e},t+\tau-\tau)\Theta(t+\tau-\tau)\braket{\sigma^{\dagger}(t)}+\mathcal{F}_{\mathrm{ab}}(\nu_{e},\tau)e^{-\left(i\omega_{00}+\gamma\right)\tau}.$

Let us pay attention to the steady state regime by setting
$t\rightarrow\infty$. According to Eq. (73b), we can get the expectation value
of the transition dipole moment beyond the rotating frame

$\lim\limits_{t\to\infty}\braket{\sigma(t)}=\chi_{\alpha}e^{-i\omega_{\ell}t}.$
(84)

Inserting the expression of the function $\mathcal{F}_{\mathrm{ab}}(\nu,t)$
given by Eq. (64) into Eq. (83), we obtain

$\displaystyle\lim\limits_{t\to\infty}\braket{\sigma(t+\tau)\sigma^{\dagger}(t)}\approx$
$\displaystyle\sum_{m,n=0}^{\infty}\frac{\eta_{\ell}^{2}S_{m}^{\mathrm{ab}}S_{n}^{\mathrm{ab}}\left[e^{-i\omega_{\ell}\tau}-e^{-\left(i\omega_{00}+im\nu_{e}+\gamma+m\Gamma\right)\tau}\right]}{\left[\gamma+m\Gamma+i\left(\omega_{00}+m\nu-\omega_{\ell}\right)\right]\left[\gamma+n\Gamma-i\left(\omega_{00}+n\nu_{e}-\omega_{\ell}\right)\right]}$
(85a)
$\displaystyle+\sum_{m=0}^{\infty}S_{m}^{\mathrm{ab}}e^{-\left(i\omega_{00}+im\nu_{e}+\gamma+m\Gamma\right)\tau}.$
(85b)

As $\eta_{\ell}\ll\gamma$, the contribution of the term in Eq. (85a) can be
ignored. Finally, we can simplify the expression of
$\braket{\sigma(t+\tau)\sigma^{\dagger}(t)}$ into

$\lim\limits_{t\to\infty}\braket{\sigma(t+\tau)\sigma^{\dagger}(t)}=\sum_{m=0}^{\infty}S_{m}^{\mathrm{ab}}e^{-\left(i\omega_{00}-im\nu_{e}+\gamma+m\Gamma\right)\tau}.$
(86)

Performing the Fourier transform, the absorption spectrum is obtained as

$S_{\mathrm{Ab}}(\omega)=2\mathcal{R}\int_{0}^{\infty}d\tau\,\lim\limits_{t\to\infty}\braket{\sigma(t+\tau)\sigma^{\dagger}(t)}e^{i\omega\tau}=\sum_{m=0}^{\infty}\frac{S_{m}^{\mathrm{ab}}\left(\gamma+m\Gamma\right)}{\left(\gamma+m\Gamma\right)^{2}+\left(\omega+m\nu_{e}-\omega_{00}\right)^{2}}.$
(87)