# The universal model of strong coupling at the nonlinear parametric resonance
in open cavity-QED systems

Mikhail Tokman  Institute of Applied Physics, Russian Academy of Sciences,
Nizhny Novgorod, 603950, Russia  Maria Erukhimova  Institute of Applied
Physics, Russian Academy of Sciences, Nizhny Novgorod, 603950, Russia  Qianfan
Chen  Department of Physics and Astronomy, Texas A&M University, College
Station, TX, 77843 USA  Alexey Belyanin  Department of Physics and Astronomy,
Texas A&M University, College Station, TX, 77843 USA

###### Abstract

Many molecular, quantum-dot, and optomechanical nanocavity-QED systems
demonstrate strong nonlinear interactions between electrons, photons, and
phonon (vibrational) modes. We show that such systems can be described by a
universal model in the vicinity of the nonlinear parametric resonance
involving all three degrees of freedom. We solve the nonperturbative quantum
dynamics in the strong coupling regime of the nonlinear resonance, taking into
account quantization, dissipation, and fluctuations of all fields. We find
analytic solutions for quantum states in the rotating wave approximation which
demonstrate tripartite quantum entanglement once the strong coupling regime is
reached. We show how the strong coupling at the nonlinear resonance modifies
photon emission and vibrational spectra, and how the observed spectra can be
used to extract information about relaxation rates and the nonlinear coupling
strength in specific systems.

##  I Introduction

Nonlinear optical interactions acquire qualitatively new features in the
strong-coupling regime of cavity quantum electrodynamics (QED), especially
when utilizing an extreme field localization achievable in nanophotonic
cavities. Even in the standard cavity QED scenario of the strong coupling at
the two-wave resonance between only two degrees of freedom, e.g., between an
electronic or vibrational transition in a molecule and an electromagnetic (EM)
cavity mode, strong coupling has been shown to modify the properties of Raman
scattering, generation of harmonics, four-wave mixing, and nonlinear
parametric interactions, with applications in photochemistry, quantum
information, and quantum sensing; see, e.g., schmidt2016 ; pino2015 ; pino2016
; ramelow2019 ; neuman2020 ; ge2021 ; k-wang2021 ; narang2021 ; wang2021 ;
zasedatelev2021 ; pscherer2021 ; hughes2021 .

The dynamics becomes more complicated but also more interesting when the
strong coupling regime is realized at the nonlinear resonance between three or
more degrees of freedom. We will call this nonlinear resonance “parametric”
for lack of a better word, although this term is often applied to nonlinear
processes involving only the EM fields, in which there is no energy exchange
between the fields and the medium, for example parametric down-conversion in
transparent nonlinear materials LLcont . In contrast, for all examples in this
paper the energy exchange between the bosonic fields and fermionic quantum
emitters is of principal importance; moreover, we are interested in the
nonperturbative regime of strong coupling when the excitation of the fermionic
degree of freedom is not small. This is the regime most interesting for
quantum technology applications as it actively involves the fermionic “qubit”
in the processes of writing, reading, and transferring the information encoded
in a quantum state.

In this paper we will focus on the nonlinear three-wave resonance between
three degrees of freedom for the sake of simpler algebra, although the
formalism does not have this limitation. One possible example where it can be
realized is a molecule or an ensemble of molecules placed in a photonic or
plasmonic nanocavity; e.g. benz2016 ; park2016 ; chikkaraddy2016 ;
pscherer2021 ; hughes2021 . In this case the fermion system may comprise two
or more electron states forming an optical transition at frequency
$\omega_{e}$, and the nonlinear parametric process may involve, e.g., a decay
of the electron excitation into a cavity photon at frequency $\omega$ and a
phonon of a given vibrational mode of a molecule at frequency $\Omega$, under
the nonlinear resonance condition $\omega_{e}=\omega+\Omega$, or an absorption
of a photon with simultaneous creation of electron and phonon excitations,
given by the nonlinear resonance condition $\omega=\omega_{e}+\Omega$. When
the strength of such a nonlinear three-wave interaction is higher than the
dissipation rates, hybrid electron-photon-phonon states are formed. If the
phonon mode is classical, the parametric process is simply the modulation of
the electron-photon coupling by molecular vibrations which serve as an
external driving force for the electron-photon quantum dynamics. If the phonon
mode is quantized, the strong coupling between photon, phonon, and electron
degrees of freedom near parametric resonance $\omega_{e}=\omega\pm\Omega$
leads inevitably to the formation of tripartite entangled states belonging to
the family of Greenberger-Horne-Zeilinger (GHZ) states tokman2021 .

Another route to the parametric resonance is within the framework of cavity
optomechanics aspelmeyer2014 ; meystre2013 ; roelli2016 ; pirkkalainen2015
and quantum acoustics chu2017 ; hong2017 ; arriola2019 . It can occur in the
situations where mechanical oscillations of a cavity parameter at frequency
$\Omega$ modulate the resonance between the optical transition in a quantum
emitter and a photon cavity mode. Here again the parametric resonance
$\omega_{e}=\omega\pm\Omega$ at strong coupling should give rise to tripartite
entangled states of the electrons, photons, and mechanical vibrations liao2018
. While quantization of all three degrees of freedom in experiment remains an
unsolved challenge, strong coupling and entanglement of acoustic phonons
satzinger2018 ; bienfait2020 , resolving the energy levels of a nanomechanical
oscillator arriola2019 , or cooling a macroscopic system into its motional
ground state delic2020  have already been demonstrated.

Yet another situation leading to parametric resonance is when a quantum
emitter such as a quantum dot or an optically active defect in a solid matrix
is coupled to the EM cavity field and the phonon modes of a crystal lattice
wilson2002 . In this case the coupling to phonons can be introduced via the
same Huang-Rhys theory as in the molecular systems weiler2012 ; roy2015 ;
denning2020 ; see the Hamiltonian Eq. (14) below. The exciton-photon coupling
strength (Rabi frequency) in nanocavities can be even high enough to exceed
the phonon frequency, which brings the system to the ultrastrong coupling
regime denning2020 .

There is a great variety of models and formalisms describing these diverse
physical systems, and attempts have been made to establish connections between
different models. For example, it was shown in roelli2016  that plasmon-
enhanced Raman scattering on individual molecules can be mapped onto a cavity
optomechanics Hamiltonian when the plasmon frequency is far detuned from the
electronic transition in a molecule, i.e., no electrons are excited. In this
case the vibrational mode of a molecule is analogous to the mechanical
oscillations of a cavity parameter. In pino2016  it was argued that Stokes
Raman scattering in molecules under the condition of a strong coupling between
the vibrational mode and the cavity mode is equivalent to the parametric decay
of the pump photon into the Stokes photon and the vibrational quantum. In both
cases a simple linear resonance $\omega=\Omega$ was assumed and the EM field
was far detuned from any electronic transition in a molecule, thus excluding
any real electron excitation. In hughes2021  resonant Raman scattering of
single molecules when the two-wave exciton-photon coupling frequency is
comparable to the vibrational frequency was analyzed. The situation when the
Rabi frequency becomes comparable to the vibrational frequency was also
considered in tokman2021 .

In this paper we deal with the strong coupling at the nonlinear parametric
resonance $\omega_{e}=\omega\pm\Omega$ when the excitation of the electron
transition and energy exchange between all three degrees of freedom are of
principal importance. We show that in the RWA (i.e., excluding ultrastrong-
coupling regimes) all physical models of electron-photon-vibrational coupling
in molecular, optomechanical, and any other coupled three-mode system can be
mapped onto the universal “parametric” Hamiltonian, independently on the
specific physical mechanism of coupling. Note that the system must be in the
RWA regime for the nonlinear parametric resonance and all associated physics
to exist and make sense. Otherwise the three-wave and two-way resonances
overlap tokman2021  and the GHZ-like entangled states cannot be created. Also,
it becomes impossible to build a universal Hamiltonian. That is why the
ultrastrong-coupling regimes are not of interest to us in this paper.

When solving for the quantum dynamics of the resulting nonlinear coupled
system, we include the effects of decoherence and coupling of each dynamic
subsystem (electrons, photons, and phonons) to its own reservoir in the Markov
approximation. Previous works (see, e.g., roy2015 ; denning2020 ) included the
non-Markovian effects in the coupling of the two-wave exciton-photon resonance
to the phonon reservoir. In our case the phonon mode, which is strongly
coupled to the exciton and photon modes through the nonlinear resonance, is
part of the dynamical system. One could say that the phonon effect on the
dynamics is “extremely non-Markovian”, except that this terminology ceases to
have any meaning in this case. The Markov approximation is of course related
only to (weak) coupling of all components of the dynamic system to their
dissipative reservoirs.

Within the formalism of the stochastic equation of evolution for the state
vector we are able to find the general analytic solution for the
nonperturbative dynamics of the open quantum system. This approach is well
known zoller1997 , but it is usually applied for numerical Monte-Carlo
simulations scully1997 ; plenio1998 ; gisin1992 ; diosi1998 ; tannoudji1993 ;
molmer1993 ; gisin1992-2 ; raedt2017 ; nathan2020 . We recently developed a
version of stochastic Schrödinger equation suitable for analytic solutions in
open strongly-coupled cavity QED problems tokman2021 ; chen2021 . We calculate
the photon and phonon emission spectra to obtain the experimentally observable
signatures of the strong coupling regime and tripartite quantum entanglement.
It is remarkable that we are able to find analytic solutions for the quantum
dynamics in systems of coupled electron, photon, and phonon excitations even
including dissipation and fluctuation effects in all subsystems. This allows
one to find the expressions for the time evolution of the state vector and
observable quantities in the form which shows explicitly the dependence on all
experimental parameters: transition energies and frequencies, matrix elements
of the optical transitions, the spatial structure of the field modes,
relaxation rates for all constituent subsystems, ambient temperatures etc. We
believe that the results obtained in this paper will be useful for designing
and interpreting the experiments on a broad range of cavity QED systems.

The paper is structured as follows. In Section II we present the Hamiltonian
for coupled quantized fermion, photon, and phonon fields near the parametric
resonance for one particular mechanism of three-wave coupling. In Section III
we show that a large variety of different three-wave coupling mechanisms and
physical systems are reduced to the same Hamiltonian which therefore can serve
as a universal model of the parametric resonance. Section IV includes the
effects of dissipation, decoherence, and fluctuations within the stochastic
equation for the state vector which describes the evolution of an open system
in contact with dissipative reservoirs. As compared to our recent work
tokman2021 ; chen2021 , we develop a model of fluctuations and dissipative
processes which includes all effects of phonon dissipation on the dynamics of
the parametric process and the emission spectra. In Section V we describe the
formation of entangled electron-photon-phonon states for an open system.
Section VI calculates the emission spectra of photons and phonons resulting
from the nonlinear parametric decay of an electron excitation. Appendix A
contains the derivation details for some results in Sec. V.

##  II The universal Hamiltonian of the nonlinear parametric resonance

Consider a simple model of three interacting quantum subsystems which includes
(1) an electron transition in a quantum emitter such as an atom, molecule,
optically active impurity, quantum dot, etc., which we will model as a 2-level
system, (2) a single-mode electromagnetic (EM) field in a cavity, and (3) a
mode of mechanical, acoustic, or molecular vibrations (“phonons”). Although in
this section we write the Hamiltonian for a specific model, in the next
section we show that the same Hamiltonian describes the nonlinear parametric
coupling in a variety of physical systems.

The generalization to many bosonic modes or fermionic degrees of freedom is
straightforward and still allows analytic solution within the rotating wave
approximation (RWA), but it leads to more cumbersome algebra (see, e.g.,
tokman2021-3 ), so we will keep only three degrees of freedom for clarity.

Figure 1 shows a generic model of parametric decay of an electron excitation
in a quantum emitter (e.g., a molecule) into a photon of a cavity mode at
frequency $\omega$ and a phonon of a given vibrational mode at frequency
$\Omega$, under the condition of the parametric resonance
$\omega_{e}\approx\omega+\Omega$.

Figure 1:  A sketch of nonlinear parametric resonance for a molecule in a
cavity showing the decay of the electron excitation at frequency $\omega_{e}$
into a cavity mode photon at frequency $\omega$ and a phonon of a given
vibrational mode at frequency $\Omega$. The relaxation rates of the electron,
photon, and vibrational excitations are $\gamma$, $\mu_{\omega}$, and
$\mu_{\Omega}$, respectively.

In the absence of coupling, the partial Hamiltonians are:

###  II.1 The 2-level fermion system

It is described by a standard effective Hamiltonian

$\hat{H}_{e}=\hbar\omega_{e}\hat{\sigma}^{\dagger}\hat{\sigma}.$ (1)

Here $\hat{\sigma}=\left|0\right\rangle\left\langle 1\right|$,
$\hat{\sigma}^{\dagger}=\left|1\right\rangle\left\langle 0\right|$ ;
$\left|0\right\rangle$ and $\left|1\right\rangle$ are the eigenstates of an
“atom” with energies $0$ and $\hbar\omega_{e}$ respectively. The Hamiltonian
(1) corresponds to the dipole moment operator

$\mathbf{\hat{d}}=\mathbf{d}\left(\hat{\sigma}^{\dagger}+\hat{\sigma}\right),$
(2)

where $\mathbf{d=-e}\left\langle 1\right|\mathbf{r}\left|0\right\rangle$,
$\mathbf{r}$ is a coordinate for the finite motion of a bound electron.

###  II.2 The EM field

Here we consider a single-mode EM field for simplicity, although including
many bosonic field modes does not present any principal difficulties. Besides,
in a micro- or nanocavity other EM modes will be far detuned from the
nonlinear resonance. The Hamiltonian is

$\hat{H}_{em}=\hbar\omega\hat{c}^{\dagger}\hat{c}.$ (3)

Here $\hat{c}$ and $\hat{c}^{\dagger}$ are standard bosonic annihilation and
creation operators of photons or plasmons in the EM mode of frequency $\omega$
. The electric field operator is

$\mathbf{\hat{E}}=\mathbf{E}\left(\mathbf{r}\right)\hat{c}+\mathbf{E}^{\ast}\left(\mathbf{r}\right)\hat{c}^{\dagger}.$
(4)

The spatial structure of the normalization amplitude of the field
$\mathbf{E}\left(\mathbf{r}\right)$ is determined by solving the boundary
value problem. The normalization condition is

$\int_{V}\frac{\partial\left[\omega^{2}\varepsilon\left(\omega,\mathbf{r}\right)\right]}{\omega\partial\omega}\mathbf{E}^{\ast}\left(\mathbf{r}\right)\mathbf{E}\left(\mathbf{r}\right)d^{3}r=4\pi\hbar\omega.$
(5)

Here $V$ is the quantization volume,
$\varepsilon\left(\omega,\mathbf{r}\right)$ the dielectric function of the
dispersive medium which fills in the resonator. Eq. (5) is derived, e.g., in
tokman2016 ; tokman2013 ; tokman2015 ; tokman2018 .

###  II.3 The phonons

We again assume a single bosonic mode of a vibrational field for the same
reasons,

$\hat{H}_{p}=\hbar\Omega\hat{b}^{\dagger}\hat{b},$ (6)

where $\hat{b}$ and $\hat{b}^{\dagger}$ are phonon annihilation and creation
operators. Depending on the situation, they may define, e.g., the radius-
vector of oscillations of the center of mass of an atom tokman2021 ;
neuman2020 ; felipe2017  or a geometric parameter of the optomechanical cavity
roelli2016 ; aspelmeyer2014 ,

$\mathbf{\hat{R}}=\mathbf{Q}\hat{b}+\mathbf{Q}^{\ast}\hat{b}^{\dagger}.$ (7)

The normalization amplitude $\mathbf{Q}$ depends on the system; its absolute
value can be expressed through an effective mass of the quantum mechanical
oscillator aspelmeyer2014 :
$\left|\mathbf{Q}\right|^{2}=\frac{\hbar}{2m_{eff}\Omega}$.

###  II.4 The coupling

The coupling between subsystems is strongest at resonance. Since usually
$\omega,\omega_{e}\gg\Omega$, two most relevant resonances are a two-wave
resonance,

$\omega_{e}\approx\omega$ (8)

and a three-wave (parametric) resonance,

$\omega_{e}\approx\omega\pm\Omega$ (9)

There could be also harmonic resonances at the harmonics of the phonon
frequency: $\omega_{e}\approx\omega\pm M\Omega$, where $M$ is integer. The
modulation of the system parameters by a classical phonon field at frequency
$\Omega$ was studied, e.g., in chen2021 , and we don’t consider the classical
field here.

The two-wave resonance in the RWA scully1997  is desrcibed by the Jaynes-
Cummings (JC) Hamiltonian jaynes1963 ,

$\hat{H}=\hbar\omega\hat{c}^{\dagger}\hat{c}+\hbar\omega_{e}\hat{\sigma}^{\dagger}\hat{\sigma}+\hbar\left(\Omega_{R}^{\left(2\right)}\hat{\sigma}^{\dagger}\hat{c}+h.c.\right).$
(10)

The electric dipole coupling between the electron and EM subsystems is
expressed here through the effective Rabi frequency for the two-wave coupling,
$\Omega_{R}^{\left(2\right)}=$
$-\frac{\mathbf{d}\cdot\mathbf{E}\left(\mathbf{r}_{0}\right)}{\hbar}$, where
$\mathbf{r}_{0}$ is the coordinate of a point-like atom.

A three-wave (parametric) resonance appears in many different scenarios. As we
show in the next section, various models existing in the literature can be
described with one universal Hamiltonian. One of the scenarios leading to the
universal three-wave coupling Hamiltonian is when a quantized phonon
(vibrational) mode modulates the coupling energy between the electron
transition and the EM field mode. In this case the JC Hamiltonian is
generalized to the following form tokman2021 ; pino2016 :

$\hat{H}=\hbar\omega\hat{c}^{\dagger}\hat{c}+\hbar\omega_{e}\hat{\sigma}^{\dagger}\hat{\sigma}+\hbar\Omega\hat{b}^{\dagger}\hat{b}+\hbar\left(\Omega_{R}^{\left(3\right)}\hat{\sigma}^{\dagger}\hat{c}\hat{b}+h.c.\right),$
(11)

where $\Omega_{R}^{\left(3\right)}$ is the coupling parameter for the three-
wave resonance, which depends on the specific coupling mechanism. The above
interaction term is written for the decay of an electron transition into the
photon and the phonon, i.e., assuming that the electron excitation energy is
the largest of the three. This decay process corresponds to the upper (plus)
sign in the resonant condition (9). If we choose the lower (minus) sign in Eq.
(9), the three-wave coupling Hamiltonian will become
$\hbar\left(\Omega_{R}^{\left(3\right)}\hat{c}^{\dagger}\hat{b}\hat{\sigma}+h.c.\right)$.

Note that the three-wave coupling term in Eq. (11) has the structure formally
equivalent to the parametric down-conversion (PDC) Hamiltonian describing the
parametric decay of the quantized pump field into quantized signal and idler
modes couteau2018 ; tokman2021-2 . Of course one difference is that all fields
in the photonic PDC process are described by bosonic operators whereas the
electron excitation in Eq. (11) is described by fermionic operators, giving
rise to its specific nonlinearities.

The strong coupling regime is realized when the three-wave coupling parameter
in Eq. (11) is larger than a certain combination of the relaxation constants
$\gamma,\mu_{\omega},\mu_{\Omega}$ of all subsystems. The exact criterion can
be retrieved from the analytic solution presented in Sections V and VI below.

The specific form of the parameter $\Omega_{R}^{\left(3\right)}$ depends on
the nonlinear coupling mechanism. For example, a phonon mode can modulate the
position of the center of mass of an “atom” within a spatially nonuniform
distribution of the EM field of a cavity mode. It can be realized for all
kinds of quantum emitters: an electron transition in a molecule, a quantum dot
or defect in a solid matrix, an optomechanical system with a varying cavity
parameter, etc. In this case, in the limit of a small amplitude of vibrations,
one can obtain that tokman2021

$\Omega_{R}^{\left(3\right)}=-\frac{1}{\hbar}\left[\mathbf{d}\left(\mathbf{Q\cdot\nabla}\right)\mathbf{E}\right]_{\mathbf{r}=\mathbf{r}_{0}}.$
(12)

The Hamiltonian in Eq. (11) is valid if the three-wave resonance is well
separated from the two-wave one. The conditions for that are tokman2021

$\left|\omega_{e}-\omega-\Omega\right|\ll\left|\omega_{e}-\omega\right|,\
\left|\Omega_{R}^{\left(2,3\right)}\right|\ll\Omega.$ (13)

In the next section we will see that if the conditions (13) are satisfied,
other models of three-wave coupling can be reduced to the universal parametric
Hamiltonian (11). Note that in plasmonic nanocavities the two-wave or/and
three-wave Rabi frequency $\Omega_{R}^{(2,3)}$ can become higher than the
vibrational or phonon frequency $\Omega$, which would violate the last of
inequalities (13); see hughes2021 ; denning2020 ; tokman2021 .

##  III The models described by the universal parametric Hamiltonian

In addition to the three-wave coupling mechanism considered in the previous
section, there are other ways for phonons or any mechanical oscillations to
affect the coupling between the EM cavity field and the quantum emitter. Here
we give several examples, assuming without loss of generality that the
amplitudes $\mathbf{Q}$ in the expression for the position displacement
operator in Eq. (7) are real functions.

###  III.1 Phonons modulate the energy of the electron transition (see, e.g.,
neuman2020 ; felipe2017 )

$\hat{H}=\hbar\omega\hat{c}^{\dagger}\hat{c}+\hbar\omega_{e}\hat{\sigma}^{\dagger}\hat{\sigma}+\hbar\Omega\hat{b}^{\dagger}\hat{b}+\hbar\left(\Omega_{R}^{\left(2\right)}\hat{\sigma}^{\dagger}\hat{c}+h.c.\right)+\hbar\sqrt{S}\Omega\hat{\sigma}^{\dagger}\hat{\sigma}\left(\hat{b}+\hat{b}^{\dagger}\right).$
(14)

Here $S$ is the Huang–Rhys factor, which determines the dependence of the
transition energy on the dimensionless amplitude $\hat{b}+\hat{b}^{\dagger}$
of the phonon oscillations. Let’s call Eq. (14) the “molecular” Hamiltonian,
although it can also describe the exciton-phonon coupling in quantum-dot
systems weiler2012 ; roy2015 ; denning2020 .

###  III.2 Phonons modulate some geometric parameter of the cavity (see, e.g.,
roelli2016 ; aspelmeyer2014 )

$\hat{H}=\hbar\omega\hat{c}^{\dagger}\hat{c}+\hbar\omega_{e}\hat{\sigma}^{\dagger}\hat{\sigma}+\hbar\Omega\hat{b}^{\dagger}\hat{b}+\hbar\left(\Omega_{R}^{\left(2\right)}\hat{\sigma}^{\dagger}\hat{c}+h.c.\right)-\hbar
g\hat{c}^{\dagger}\hat{c}\left(\hat{b}+\hat{b}^{\dagger}\right).$ (15)

Here the factor $g$ determines the linear dependence of the resonant frequency
of the cavity on some geometric parameter $\mathbf{G}$ modulated by mechanical
vibrations:

$g=-\mathbf{Q}\frac{\partial\omega}{\partial\mathbf{G}}.$

We will call Eq. (15) the “optomechanical” Hamiltonian.

The Hamiltonians (14) and (15) appear to be very different from Eq. (11).
Indeed, they both contain the standard two-wave resonance and their three-wave
coupling terms are different. We will show now that when the conditions (13)
are satisfied, the “molecular” and “optomechanical” Hamiltonians are
equivalent to the universal Hamiltonian in Eq. (11). To prove this statement,
we write the Hamiltonian (11) in the interaction representation:

$\hat{H}_{int}=e^{i\hat{H}_{0}t}\hat{V}e^{-i\hat{H}_{0}t},$ (16)

where

$\hat{H}_{0}=\hbar\omega\hat{c}^{\dagger}\hat{c}+\hbar\omega_{e}\hat{\sigma}^{\dagger}\hat{\sigma}+\hbar\Omega\hat{b}^{\dagger}\hat{b},\
\
\hat{V}=\hbar\left(\Omega_{R}^{\left(3\right)}\hat{\sigma}^{\dagger}\hat{c}\hat{b}+h.c.\right).$

This yields

$\hat{H}_{int}=\hbar\left(\Omega_{R}^{\left(3\right)}\hat{\sigma}^{\dagger}\hat{c}\hat{b}e^{i\left(\omega_{e}-\omega-\Omega\right)t}+h.c.\right).$
(17)

Next, we write the Hamiltonians in Eqs. (14) and (15) in the interaction
representation by defining the unperturbed Hamiltonian as

$\hat{H}_{0}=\hbar\omega\hat{c}^{\dagger}\hat{c}+\hbar\omega_{e}\hat{\sigma}^{\dagger}\hat{\sigma}+\hbar\Omega\hat{b}^{\dagger}\hat{b}+\hbar\left(\Omega_{R}^{\left(2\right)}\hat{\sigma}^{\dagger}\hat{c}+h.c.\right).$
(18)

This gives

$\hat{V}=\hat{V}_{mol}=\hbar\sqrt{S}\Omega\hat{\sigma}^{\dagger}\hat{\sigma}\left(\hat{b}+\hat{b}^{\dagger}\right)$
(19)

for the “molecular” Hamiltonian and

$\hat{V}=\hat{V}_{optom}=-\hbar
g\hat{c}^{\dagger}\hat{c}\left(\hat{b}+\hat{b}^{\dagger}\right)$ (20)

for the optomechanical one.

The operator $\hat{H}_{0}$ given by Eq. (18) is not diagonal. In this case one
should generally diagonalize $\hat{H}_{0}$. However, sometimes a slightly
different approach is simpler. Indeed, consider the Hamiltonian given by

$\hat{H}=\hat{H}_{0}\left(\hat{A}_{1,}\cdots\hat{A}_{N}\right)+\hat{V}\left(\hat{A}_{1,}\cdots\hat{A}_{N}\right),$
(21)

where $\hat{A}_{i}$ are certain operators related to coupled subsystems (here
we assume that Hermitian conjugated operators are assigned different numbers
“$i$” ). For any interaction operator $\hat{V}$, which can be expanded in a
series

$\hat{V}=\sum_{j}k_{j}\prod_{i=1}^{N_{j}}\left(\hat{A}_{i}\right)^{n_{ji}}$

(here $n_{ji}$ are positive integers), the following relationships are true:

$\displaystyle\hat{H}_{int}$ $\displaystyle=$ $\displaystyle
e^{i\hat{H}_{0}t}\hat{V}\left(\hat{A}_{1,}\cdots\hat{A}_{N}\right)e^{-i\hat{H}_{0}t}$
$\displaystyle=$
$\displaystyle\sum_{j}k_{j}\prod_{i=1}^{N_{j}}e^{i\hat{H}_{0}t}\left(\hat{A}_{i}\right)^{n_{ji}}e^{-i\hat{H}_{0}t}$
$\displaystyle=$
$\displaystyle\sum_{j}k_{j}\prod_{i=1}^{N_{j}}\left(e^{i\hat{H}_{0}t}\hat{A}_{i}e^{-i\hat{H}_{0}t}\right)^{n_{ji}},$

from which one obtains

$\hat{H}_{int}=\hat{V}\left(\widehat{\tilde{A}}_{1,}\cdots\widehat{\tilde{A}}_{N,}\right),$
(22)

where

$\widehat{\tilde{A}}_{i}\left(t,\hat{A}_{1,}\cdots\hat{A}_{N}\right)=e^{i\hat{H}_{0}t}\hat{A}_{i}e^{-i\hat{H}_{0}t}$
(23)

The operators $\widehat{\tilde{A}}_{i}$ satisfy the Heisenberg equations

$\frac{\partial\widehat{\tilde{A}}_{i}}{\partial
t}=\frac{i}{\hbar}\left[\hat{H}_{0},\widehat{\tilde{A}}_{i}\right]$ (24)

for the initial conditions
$\widehat{\tilde{A}}_{i}\left(t=0\right)=\hat{A}_{i}$. In particular, for the
Hamiltonian $\hat{H}_{0}$ given by Eq. (18) one obtains

$\frac{\partial\widehat{\tilde{\sigma}}}{\partial
t}=-i\omega_{e}\widehat{\tilde{\sigma}}-i\Omega_{R}^{\left(2\right)}\widehat{\tilde{c}}\left(1-2\widehat{\tilde{\sigma}}^{\dagger}\widehat{\tilde{\sigma}}\right),$
(25) $\frac{\partial\widehat{\tilde{c}}}{\partial
t}=-i\omega\widehat{\tilde{c}}-i\Omega_{R}^{\left(2\right)\ast}\widehat{\tilde{\sigma}},$
(26)

where we used the integral of motion
$\widehat{\tilde{\sigma}}^{\dagger}\widehat{\tilde{\sigma}}+\widehat{\tilde{\sigma}}\widehat{\tilde{\sigma}}^{\dagger}=1$.
Taking into account the condition
$\left|\Omega_{R}^{\left(2\right)}\right|\ll\left|\omega_{e}-\omega\right|$
which follows from Eqs. (13), when solving for the operators
$\widehat{\tilde{c}}$ and $\widehat{\tilde{\sigma}}$ one can neglect the terms
of the order of
$\left|\frac{\Omega_{R}^{\left(2\right)}}{\omega_{e}-\omega}\right|^{2}$. For
the initial conditions $\widehat{\tilde{\sigma}}$
$\left(t=0\right)=\hat{\sigma}$ and
$\widehat{\tilde{c}}\left(t=0\right)=\hat{c}$ the solution expanded in series
in powers of the small parameter
$\left|\frac{\Omega_{R}^{\left(2\right)}}{\omega_{e}-\omega}\right|$ takes the
form

$\left(\begin{array}[]{c}\widehat{\tilde{\sigma}}\\\
\widehat{\tilde{c}}\end{array}\right)=\widehat{\widehat{M}}\left(\begin{array}[]{c}\hat{\sigma}\\\
\hat{c}\end{array}\right),$ (27)

where

$\displaystyle\widehat{\widehat{M}}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}e^{-i\omega_{e}t}&0\\\ 0&e^{-i\omega
t}\end{array}\right)$ (33)
$\displaystyle-\left(\begin{array}[]{cc}0&\frac{\Omega_{R}^{\left(2\right)}}{\omega_{e}-\omega}\left(1-2\hat{\sigma}^{\dagger}\hat{\sigma}\right)\left(e^{-i\omega
t}-e^{-i\omega_{e}t}\right)\\\
\frac{\Omega_{R}^{\left(2\right)\ast}}{\omega_{e}-\omega}\left(e^{-i\omega
t}-e^{-i\omega_{e}t}\right)&0\end{array}\right)+o\left(\left|\frac{\Omega_{R}^{\left(2\right)}}{\omega_{e}-\omega}\right|^{2}\right)$

There is an exact solution for the operator $\widehat{\tilde{b}}$:

$\widehat{\tilde{b}}=e^{i\hat{H}_{0}t}\hat{b}e^{-i\hat{H}_{0}t}=\hat{b}e^{-i\Omega
t}.$ (34)

Now we substitute the three-wave coupling Hamiltonian (19) into Eq. (22) and
use Eqs. (27)-(34). The conditions (13) combined with the RWA allow one to
keep only slowly varying terms $\propto
e^{i\left(\omega_{e}-\omega-\Omega\right)t}$ in the final expression. Taking
into account that
$\hat{\sigma}^{\dagger}\hat{\sigma}^{\dagger}=\hat{\sigma}\hat{\sigma}$ $=0$
and taking $\frac{1}{\omega_{e}-\omega}\approx\frac{1}{\Omega}$ in Eq. (33),
we obtain the following expression for the “ molecular” Hamiltonian in the
interaction picture:

$\left(\hat{H}_{int}\right)_{mol}=-\hbar\left(\sqrt{S}\Omega_{R}^{\left(2\right)}\hat{\sigma}^{\dagger}\hat{c}\hat{b}e^{i\left(\omega_{e}-\omega-\Omega\right)t}+h.c.\right)$
(35)

A similar derivation for the “ optomechanical” Hamiltonian given by Eq. (20)
leads to the following result:

$\left(\hat{H}_{int}\right)_{optom}=-\hbar\left(\frac{g\Omega_{R}^{\left(2\right)}}{\Omega}\hat{\sigma}^{\dagger}\hat{c}\hat{b}e^{i\left(\omega_{e}-\omega-\Omega\right)t}+h.c.\right)$
(36)

Clearly, in both cases the Hamiltonian has the same structure as the
parametric Hamiltonian (17), in which

$\left(\Omega_{R}^{\left(3\right)}\right)_{mol}=-\sqrt{S}\Omega_{R}^{\left(2\right)}\
\ {\rm or}\ \
\left(\Omega_{R}^{\left(3\right)}\right)_{optom}=-\frac{g}{\Omega}\Omega_{R}^{\left(2\right)}.$
(37)

Therefore, one can use the universal parametric Hamiltonian (11) for all kinds
of three-wave couplings after choosing an appropriate expression for the
coupling parameter $\Omega_{R}^{\left(3\right)}$.

We emphasize again that the universal character of the parametric Hamiltonian
(11) holds as long as inequalities (13) are satisfied, which ensure that the
three-wave nonlinear resonance can be separated from the two-wave resonance.

##  IV Including dissipation and fluctuations within the stochastic equation
for the state vector

###  IV.1 Stochastic equation for the state vector

Many quantum information applications are based on the strong coupling regime
in which the relaxation times $\tau$ are much longer than the dynamical
coupling times $T$ between the subsystems. For two- and three-wave couplings
those times are determined by effective Rabi frequencies, $T$ -1
$\sim\left|\Omega_{R}^{\left(2,3\right)}\right|$. As shown in tokman2021 ;
chen2021 , the method of the stochastic equation for the state vector is often
the most convenient way to describe the nonperturbative dynamics of open
strongly coupled systems; it leads to simpler derivations for the observables
and characterization of entanglement than the operator-valued Heisenberg-
Langevin equations or the master equation for the density matrix.

Indeed, when applied to strongly coupled systems the Heisenberg approach leads
to the nonlinear operator-valued equations even in the simplest case of a
single two-level atom coupled to a single-mode field scully1997 . In contrast,
the equations for the state vector components are always linear; they contain
much fewer variables as compared to the density matrix equations and are split
into low-dimensional blocks in the RWA, leading to analytic solutions for both
two-wave scully1997  and three-wave tokman2021 ; chen2021  resonant coupling.

One potential difficulty with this approach is that dissipation and
fluctuations may lead to the coupling between different blocks of equations
for the state vector components that were uncoupled in a closed system.
However, in the strong coupling regime the coupling through dissipative
reservoirs is weak (scales as a small parameter $T/\tau$) and can be taken
into account perturbatively tokman2021 ; chen2021 .

Following tokman2021 ; chen2021 , we apply the stochastic equation for the
state vector to derive analytic solution for the parametric coupling of a two-
level fermionic quantum emitter to two boson fields. The stochastic equation
has the following general form,

$\frac{d}{dt}\left|\Psi\right\rangle=-\frac{i}{\hbar}\hat{H}_{eff}\left|\Psi\right\rangle-\frac{i}{\hbar}\left|\mathfrak{R}\right\rangle.$
(38)

Here $\left|\Psi\right\rangle$ is the state vector;
$\left|\mathfrak{R}\right\rangle$ is the noise term satisfying
$\overline{\left|\mathfrak{R}\right\rangle}=0$, where the bar means averaging
over the noise statistics; $\hat{H}_{eff}=$
$\hat{H}+\hat{H}^{\left(ah\right)}$ is an effective Hamiltonian which is a
non-Hermitian operator. Its non-Hermitian component
$\hat{H}^{\left(ah\right)}$ describes the effects of relaxation. The
expressions for $\hat{H}^{\left(ah\right)}$ and
$\left|\mathfrak{R}\right\rangle$ must be consistent with each other to
guarantee the conservation of the noise-averaged norm,
$\overline{\left\langle\Psi\left(t\right)\right.\left|\Psi\left(t\right)\right\rangle}=1$,
and ensure that the system reaches a physically meaningful steady state in the
absence of any external driving force. To calculate the observables from the
state vector given by Eq. (38) one should apply a standard procedure but with
an important extra step: averaging over the noise statistics, i.e.,
$q=\overline{\left\langle\Psi\right|\hat{q}\left|\Psi\right\rangle}$ , where
$\hat{q}$ is a quantum-mechanical operator corresponding to the observable
$q$.

Perhaps the most popular version of the stochastic approach to derive the
state vector, i.e., the stochastic Schrödinger equation (SSE), is its
application for numerical Monte-Carlo simulations within the method of quantum
jumps scully1997 ; plenio1998 ; gisin1992 ; diosi1998 ; tannoudji1993 ;
molmer1993 ; gisin1992-2 ; raedt2017 ; nathan2020 . The stochastic equation in
a different form, the Schrödinger-Langevin equation (SLE), was suggested to
describe the Brownian motion of a quantum particle in a constant field
kostin1972 ; katz2016 . Generally, using some version of the stochastic
equation fits within the narrative of the Langevin method WM94 . Within the
Langevin approach which describes the system with stochastic equations of
evolution, the averaging over the reservoir degrees of freedom is equivalent
to averaging over the statistics of the noise sources landau1965 . This
paradigm allows one to describe open systems without relying on the density
matrix.

###  IV.2 Comparison with the Lindblad formalism

It was shown in tokman2021  that one can choose the form of
$\hat{H}^{\left(ah\right)}$ and $\left|\mathfrak{R}\right\rangle$ in such a
way that the observables calculated with Eq. (38) will coincide with those
obtained by solving the master equation in the Lindblad approximation. The
corresponding master equation has the form scully1997

$\frac{d}{dt}\hat{\rho}=-\frac{i}{\hbar}\left[\hat{H},\hat{\rho}\right]+\hat{L}\left(\hat{\rho}\right)$
(39)

where $\hat{L}\left(\hat{\rho}\right)$ is the relaxation operator
(Lindbladian) which can be represented as

$\hat{L}\left(\hat{\rho}\right)=-\frac{i}{\hbar}\left(\hat{H}^{\left(ah\right)}\hat{\rho}-\hat{\rho}\hat{H}^{\left(ah\right)\dagger}\right)+\delta\hat{L}\left(\hat{\rho}\right).$
(40)

The equivalence (in the above sense) between the stochastic equation and the
Lindblad approach exists if we substitute the anti-Hermitian part of the
Hamiltonian from Eq. (40) into Eq. (38), and postulate the following
correlation properties for the noise source:

$\overline{\left|\mathfrak{R}\left(t^{\prime}\right)\right\rangle\left\langle\mathfrak{R}\left(t^{\prime\prime}\right)\right|}=\hbar^{2}\delta\left(t^{\prime}-t^{\prime\prime}\right)\delta\hat{L}\left(\hat{\rho}\right)_{\hat{\rho}\Longrightarrow\overline{\left|\Psi\right\rangle\left\langle\Psi\right|}}$
(41)

###  IV.3 Parametric decay of the electron excitation into a photon and a
phonon

We will describe the dynamics near the three-wave electron-photon-phonon
resonance by the stochastic equation (38) with the parametric Hamiltonian
(11). We will seek the state vector in the form

$\Psi=\sum_{n,\alpha=0}^{\infty,\infty}\left(C_{\alpha
n0}\left|\alpha\right\rangle\left|n\right\rangle\left|0\right\rangle+C_{\alpha
n1}\left|\alpha\right\rangle\left|n\right\rangle\left|1\right\rangle\right),$

where the order of indices corresponds to

$C_{\mathrm{phonon\,photon\,fermion}}\left|phonon\right\rangle\left|photon\right\rangle\left|fermion\right\rangle.$

Consider the initial state with an excited electron and no bosonic
excitations, of the type
$\left|\Psi\left(0\right)\right\rangle=\left|0\right\rangle\left|0\right\rangle\left|1\right\rangle$
. Near the resonance $\omega_{e}\approx\omega+\Omega$ the three-wave coupling
leads to the excitation of the state
$\left|1\right\rangle\left|1\right\rangle\left|0\right\rangle$. In the zero-
temperature limit, which is valid when the reservoir temperature is much lower
than the transition frequency (for optical frequencies it is satisfied even at
room temperature), relaxation processes could populate only the states with
lower energies:
$\left|0\right\rangle\left|1\right\rangle\left|0\right\rangle$,
$\left|1\right\rangle\left|0\right\rangle\left|0\right\rangle$ and
$\left|0\right\rangle\left|0\right\rangle\left|0\right\rangle$. Therefore, for
this initial condition the state vector will have 5 components:

$\displaystyle\left|\Psi\left(t\right)\right\rangle$ $\displaystyle=$
$\displaystyle
C_{000}\left(t\right)\left|0\right\rangle\left|0\right\rangle\left|0\right\rangle+C_{010}\left(t\right)\left|0\right\rangle\left|1\right\rangle\left|0\right\rangle+C_{100}\left(t\right)\left|1\right\rangle\left|0\right\rangle\left|0\right\rangle$
(42)
$\displaystyle+C_{110}\left(t\right)\left|1\right\rangle\left|1\right\rangle\left|0\right\rangle+C_{001}\left(t\right)\left|0\right\rangle\left|0\right\rangle\left|1\right\rangle.$

Note that we are not restricting the basis in any way and only considering the
initial conditions leading to single photon and phonon states to simplify
algebra: see, for example, the Appendix in tokman2021-3  where arbitrary
multiphoton states are considered in the same way, leading of course to more
cumbersome expressions. Another reason to consider such initial conditions is
that single-photon states are used in most applications, whereas generating
multiphoton number states remains a major challenge.

To determine the anti-Hermitian part $\hat{H}^{\left(ah\right)}$ of the
Hamiltonian which describes relaxation and the correlator
$\overline{\left|\mathfrak{R}\left(t^{\prime}\right)\right\rangle\left\langle\mathfrak{R}\left(t^{\prime\prime}\right)\right|}$
describing fluctuations, we will use the expression for the total Lindbladian
of the system including a 2-level “atom”, photons, and phonons scully1997 .
For simplicity we will assume zero temperature of dissipative reservoirs,
which is satisfied at $T\ll\hbar\omega$. The finite-temperature expressions
are given in tokman2021 ; chen2021  . Then the Lindbladian is

$L\left(\hat{\rho}\right)=L_{e}\left(\hat{\rho}\right)+L_{em}\left(\hat{\rho}\right)+L_{p}\left(\hat{\rho}\right)$
(43)
$L_{e}\left(\hat{\rho}\right)=-\frac{\gamma}{2}\left(\hat{\sigma}^{\dagger}\hat{\sigma}\hat{\rho}+\hat{\rho}\hat{\sigma}^{\dagger}\hat{\sigma}-2\hat{\sigma}\hat{\rho}\hat{\sigma}^{\dagger}\right)$
(44)
$L_{em}\left(\hat{\rho}\right)=-\frac{\mu_{\omega}}{2}\left(\hat{c}^{\dagger}\hat{c}\hat{\rho}+\hat{\rho}\hat{c}\hat{c}^{\dagger}-2\hat{c}\hat{\rho}\hat{c}^{\dagger}\right)$
(45)
$L_{p}\left(\hat{\rho}\right)=-\frac{\mu_{\Omega}}{2}\left(\hat{b}^{\dagger}\hat{b}\hat{\rho}+\hat{\rho}\hat{b}\hat{b}^{\dagger}-2\hat{b}\hat{\rho}\hat{b}^{\dagger}\right)$
(46)

where $\gamma$, $\mu_{\omega}$ and $\mu_{\Omega}$ are relaxation rates of
corresponding subsystems.

Then the stochastic equation for the state vector takes the form

$\left(\begin{array}[]{ccc}\frac{d}{dt}&0&0\\\
0&\frac{d}{dt}+i\omega_{010}+\gamma_{010}&0\\\
0&0&\frac{d}{dt}+i\omega_{100}+\gamma_{100}\end{array}\right)\times\left(\begin{array}[]{c}C_{000}\\\
C_{010}\\\
C_{100}\end{array}\right)=-\frac{i}{\hbar}\left(\begin{array}[]{c}\mathfrak{R}_{000}\\\
\mathfrak{R}_{010}\\\ \mathfrak{R}_{100}\end{array}\right),$ (47)

$\left(\begin{array}[]{cc}\frac{d}{dt}+i\omega_{110}+\gamma_{110}&i\Omega_{R}^{\left(3\right)\ast}\\\
i\Omega_{R}^{\left(3\right)}&\frac{d}{dt}+i\omega_{001}+\gamma_{001}\end{array}\right)\left(\begin{array}[]{c}C_{110}\\\
C_{001}\end{array}\right)=-\frac{i}{\hbar}\left(\begin{array}[]{c}\mathfrak{R}_{110}\\\
\mathfrak{R}_{001}\end{array}\right),$ (48)

where

$\mathfrak{R}_{\alpha ni}=\left\langle\alpha
ni\right.\left|\mathfrak{R}\right\rangle;$ $\omega_{010}=\omega,\
\omega_{100}=\Omega,\ \omega_{110}=\omega+\Omega,\ \omega_{001}=\omega_{e};$
$\gamma_{010}=\frac{1}{2}\mu_{\omega},\ \gamma_{100}=\frac{1}{2}\mu_{\Omega},\
\gamma_{110}=\frac{1}{2}\left(\mu_{\omega}+\mu_{\Omega}\right),\
\gamma_{001}=\frac{1}{2}\gamma.$

The correlators of the noise sources in Eqs. (47) and (48) are

$\overline{\mathfrak{R}_{\alpha
ni}^{\ast}\left(t^{\prime}\right)\mathfrak{R}_{\beta
mj}\left(t^{\prime\prime}\right)}=\hbar^{2}D_{\alpha ni,\beta
mj}(t^{\prime})\delta\left(t^{\prime}-t^{\prime\prime}\right),$ (49)

where the quantities $D_{\alpha ni,\beta mj}$ are determined using Eqs. (41)
and (43)-(46). For the diagonal elements of the correlators we obtain

$\displaystyle D_{110,110}$ $\displaystyle=$ $\displaystyle D_{001,001}=0$
$\displaystyle D_{100,100}$ $\displaystyle=$
$\displaystyle\mu_{\omega}\overline{\left|C_{110}\right|^{2}}$ (50)
$\displaystyle D_{010,010}$ $\displaystyle=$
$\displaystyle\mu_{\Omega}\overline{\left|C_{110}\right|^{2}}$ $\displaystyle
D_{000,000}$ $\displaystyle=$
$\displaystyle\gamma\overline{\left|C_{001}\right|^{2}}+\mu_{\omega}\overline{\left|C_{010}\right|^{2}}+\mu_{\Omega}\overline{\left|C_{100}\right|^{2}}.$

The off-diagonal elements are given by

$\displaystyle D_{\alpha ni,\beta mj}$ $\displaystyle=$ $\displaystyle
D_{\beta mj,\alpha ni}^{\ast}$ $\displaystyle D_{110,\alpha ni}$
$\displaystyle=$ $\displaystyle D_{001,\alpha ni}=D_{100,010}=0$ (51)
$\displaystyle D_{000,100}$ $\displaystyle=$
$\displaystyle\mu_{\omega}\overline{C_{010}^{\ast}C_{110}}$ $\displaystyle
D_{000,010}$ $\displaystyle=$
$\displaystyle\mu_{\Omega}\overline{C_{100}^{\ast}C_{110}}.$

The dependence of quantities $D_{\alpha ni,\beta mj}$ in the right-hand side
of Eq. (49) on time $t^{\prime}$ is due to the time dependence of amplitudes
$C_{\alpha ni}$ which enter Eqs. (50) and (51).

The derivation of the stochastic equation for the state vector including pure
dephasing processes and the finite temperature of the reservoirs has been
discussed in tokman2021 ; chen2021 .

Eqs. (48) describe the dynamic generation of an entangled state of the type
$\left|MIX\right\rangle=A\left(t\right)\left|0\right\rangle\left|0\right\rangle\left|1\right\rangle+B\left(t\right)\left|1\right\rangle\left|1\right\rangle\left|0\right\rangle$
whereas Eqs. (47) describe the relaxation dynamics leading to relaxation of
populations to states with lower energies. The quantities $D_{010,010}$ and
$D_{100,100}$ in Eqs. (50) are associated with processes of the relaxation to
states $\left|0\right\rangle\left|1\right\rangle\left|0\right\rangle$ and
$\left|1\right\rangle\left|0\right\rangle\left|0\right\rangle$ from the
entangled state $\left|MIX\right\rangle$ . The structure of the expression for
$D_{000,000}$ corresponds to the relaxation of the system from the entangled
state to the ground state via both “direct” pathway
$\left|0\right\rangle\left|0\right\rangle\left|1\right\rangle\rightarrow\left|0\right\rangle\left|0\right\rangle\left|0\right\rangle$
and multistep pathways
$\left|1\right\rangle\left|1\right\rangle\left|0\right\rangle\rightarrow\left|0\right\rangle\left|1\right\rangle\left|0\right\rangle\rightarrow\left|0\right\rangle\left|0\right\rangle\left|0\right\rangle$
and
$\left|1\right\rangle\left|1\right\rangle\left|0\right\rangle\rightarrow\left|1\right\rangle\left|0\right\rangle\left|0\right\rangle\rightarrow\left|0\right\rangle\left|0\right\rangle\left|0\right\rangle$,
as illustrated in Figs. 3 and 5 below.

###  IV.4 Expressions for noise sources in the stochastic equation

In addition to the expressions for noise correlators, it is convenient to know
more detailed expressions for the random functions describing the noise
sources $\mathfrak{R}_{\alpha ni}(t)$.

The effect of the reservoir on the dynamic system is characterized by the
matrix elements of the operator which determines coupling to the reservoir.
For a weak coupling these matrix elements are linear with respect to the
matrix elements of the operators describing the dynamical system. Therefore,
the functions $\mathfrak{R}_{\alpha ni}(t)$ should depend linearly on the
components of the state vector of the system.

Here we again consider the low-temperature case when the relaxation processes
can bring the populations only down, not up. When the reservoirs for each
subsystem are statistically independent, one can try the following ansatz,

$\mathfrak{R}_{\alpha ni}\left(t\right)=\hbar\sqrt{\gamma}C_{\alpha
n\left(i+1\right)}\left(t\right)f_{e}\left(t\right)+\hbar\sqrt{\mu_{\omega}}C_{\alpha\left(n+1\right)i}\left(t\right)f_{em}\left(t\right)+\hbar\sqrt{\mu_{\Omega}}C_{\left(\alpha+1\right)ni}\left(t\right)f_{p}\left(t\right),$
(52)

where $\overline{f_{e,em,p}}=0.$ Here $f_{e,em,p}\left(t\right)$ are random
functions that are determined by the statistics of noise in unperturbed
electron, photon, and phonon reservoirs. The functions $f_{e,em,p}$ should not
depend on the set of variables $C_{\alpha ni}$ within the above
approximations.

The linear dependence of the noise term on the state vector was also assumed
in SLE kostin1972 ; katz2016 , in which the noise terms had the form

$\left|\mathfrak{R}\right\rangle=\hat{U}\left(\mathbf{r},t\right)\left|\Psi\right\rangle,$
(53)

where $\hat{U}\left(\mathbf{r},t\right)$ is the fluctuating component of the
potential.

To ensure that Eq. (52) leads to the correlators Eqs. (49)-(51) that are
consistent with the Lindblad master equation, one needs first to define the
correlators of the random functions in Eq. (52) in the following way:

$\overline{f_{\kappa}^{\ast}\left(t^{\prime}\right)f_{\lambda}\left(t^{\prime\prime}\right)}=\delta_{\kappa\lambda}\delta\left(t^{\prime}-t^{\prime\prime}\right),$
(54)

where $\kappa,\lambda=e,em,p$, i.e. the fluctuations in different reservoirs
are independent and Markovian. Second, one has to assume that the correlations
are factorized when calculating the averages

$\overline{C_{\alpha ni}^{\ast}\left(t^{\prime}\right)C_{\beta
mj}\left(t^{\prime\prime}\right)f_{\kappa}^{\ast}\left(t^{\prime}\right)f_{\lambda}\left(t^{\prime\prime}\right)}=\overline{C_{\alpha
ni}^{\ast}\left(t^{\prime}\right)C_{\beta
mj}\left(t^{\prime\prime}\right)}\times\overline{f_{\kappa}^{\ast}\left(t^{\prime}\right)f_{\lambda}\left(t^{\prime\prime}\right)}.$
(55)

Eq. (54) looks obvious, whereas the factorization in Eq. (55) is valid only in
linear approximation with respect to relaxation constants $\gamma$ and
$\mu_{\omega,\Omega}$. However, the Lindbladian of the form given in Eqs.
(42)-(45) is itself valid within the same approximation. Therefore, Eqs. (52),
(54), and (55) lead to all expressions in Eqs. (50),(51). One also has to keep
in mind that with our choice of our initial conditions the amplitudes
$C_{\alpha ni}=0$ for all states with energies above those in states
$\left|1\right\rangle\left|1\right\rangle\left|0\right\rangle$ and
$\left|0\right\rangle\left|0\right\rangle\left|1\right\rangle$.

##  V Dynamics of entangled fermion-photon-phonon states in a dissipative
system

Here we write an explicit solution of Eqs (47) and (48) for the initial state
vector
$\left|\Psi\left(0\right)\right\rangle=\left|0\right\rangle\left|0\right\rangle\left|1\right\rangle$,
when $C_{001}\left(0\right)=1$,
$C_{000}\left(0\right)=C_{010}\left(0\right)=C_{100}\left(0\right)=C_{110}\left(0\right)=0$.
Assuming exact resonance at $\omega_{e}=\omega+\Omega$, we obtain

$\displaystyle C_{000}\left(t\right)$ $\displaystyle=$
$\displaystyle-\frac{i}{\hbar}\int_{0}^{t}\mathfrak{R}_{000}\left(t^{\prime}\right)dt^{\prime}$
$\displaystyle C_{010}\left(t\right)$ $\displaystyle=$
$\displaystyle-\frac{i}{\hbar}\int_{0}^{t}\mathfrak{R}_{010}\left(t^{\prime}\right)e^{-i\omega\left(t-t^{\prime}\right)-\gamma_{010}\left(t-t^{\prime}\right)}dt^{\prime}$
(56) $\displaystyle C_{100}\left(t\right)$ $\displaystyle=$
$\displaystyle-\frac{i}{\hbar}\int_{0}^{t}\mathfrak{R}_{100}\left(t^{\prime}\right)e^{-i\Omega\left(t-t^{\prime}\right)-\gamma_{100}\left(t-t^{\prime}\right)}dt^{\prime},$

$\displaystyle C_{001}\left(t\right)$ $\displaystyle=$ $\displaystyle
e^{-i\omega_{e}t}e^{-\frac{\gamma_{110}+\gamma_{001}}{2}t}\left[\cos\left(\tilde{\Omega}_{R}t\right)+\frac{\gamma_{110}-\gamma_{001}}{\tilde{\Omega}_{R}}\sin\left(\tilde{\Omega}_{R}t\right)\right]+\delta
C_{001}$ (57) $\displaystyle C_{110}\left(t\right)$ $\displaystyle=$
$\displaystyle
e^{-i\omega_{e}t}e^{-\frac{\gamma_{110}+\gamma_{001}}{2}t}\left(ie^{-i\theta}\right)\sin\left(\tilde{\Omega}_{R}t\right)+\delta
C_{110}$ (58)

where the effective Rabi frequency
$\tilde{\Omega}_{R}=\sqrt{\left|\Omega_{R}^{\left(3\right)}\right|^{2}-\frac{\left(\gamma_{110}-\gamma_{001}\right)^{2}}{4}}$,
$\theta=\mathrm{Arg}\left[\Omega_{R}^{\left(3\right)}\right]$, and the terms
$\delta C_{100,110}$ are linear with respect to random functions
$\mathfrak{R}_{001}$ and $\mathfrak{R}_{110}$.

The term proportional to
$\frac{\gamma_{110}-\gamma_{001}}{\tilde{\Omega}_{R}}$ in the first of Eqs.
(57) can be omitted when calculating most observables when the dissipation is
weak, $\tilde{\Omega}_{R}\gg$ $\gamma_{ani}$ (see, e.g., tokman2021 ).
However, one has to keep in mind that this term is needed for Eqs. (57) to
satisfy an exact integral of motion of Eqs. (48) :

$\frac{d}{dt}\left(\overline{\left|C_{001}\right|^{2}}+\overline{\left|C_{110}\right|^{2}}\right)=-2\gamma_{001}\overline{\left|C_{001}\right|^{2}}-2\gamma_{110}\overline{\left|C_{110}\right|^{2}}.$

Note that when averaged over the period $\frac{2\pi}{\tilde{\Omega}_{R}}$ the
integral is conserved even without this term.

The solution for the 5-component state vector can be written as

$\displaystyle\left|\Psi\right\rangle$ (59) $\displaystyle=$ $\displaystyle
e^{-i\omega_{e}t}e^{-\frac{\gamma_{110}+\gamma_{001}}{2}t}\left\\{\left[\cos\left(\tilde{\Omega}_{R}t\right)+\frac{\gamma_{110}-\gamma_{001}}{2\tilde{\Omega}_{R}}\sin\left(\tilde{\Omega}_{R}t\right)\right]\left|0\right\rangle\left|0\right\rangle\left|1\right\rangle+ie^{-i\theta}\sin\left(\tilde{\Omega}_{R}t\right)\left|1\right\rangle\left|1\right\rangle\left|0\right\rangle\right\\}$
$\displaystyle+\delta
C_{001}\left|0\right\rangle\left|0\right\rangle\left|1\right\rangle+\delta
C_{110}\left|1\right\rangle\left|1\right\rangle\left|0\right\rangle+C_{000}\left|0\right\rangle\left|0\right\rangle\left|0\right\rangle+C_{100}\left|1\right\rangle\left|0\right\rangle\left|0\right\rangle+C_{010}\left|0\right\rangle\left|1\right\rangle\left|0\right\rangle,$

where

$\overline{\delta C_{001}}=\overline{\delta
C_{110}}=\overline{C_{000}}=\overline{C_{100}}=\overline{C_{010}}=0,$ (60)
$\overline{\left|\delta C_{001}\right|^{2}}=\overline{\left|\delta
C_{110}\right|^{2}}=0.$ (61)

It follows from Eq. (59) that in the entangled state
$\left|MIX\right\rangle=A\left(t\right)\left|0\right\rangle\left|0\right\rangle\left|1\right\rangle+B\left(t\right)\left|1\right\rangle\left|1\right\rangle\left|0\right\rangle$
the amplitudes $A\left(t\right)$ and $B\left(t\right)$ oscillate at the
effective Rabi frequency and decay with the decay rate

$\gamma_{{}_{MIX}}=\frac{\gamma_{110}+\gamma_{001}}{2}=\frac{1}{4}\left(\mu_{\omega}+\mu_{\Omega}+\gamma\right).$

The occupation probabilities $|C_{001}|^{2}$ and $|C_{110}|^{2}$ are plotted
in Fig. 2 as a function of normalized time $\tilde{\Omega}_{R}t$, along with
the real parts of their eigenfrequencies obtained from Eqs. (48) as a function
of detuning from the parametric resonance $\omega+\Omega-\omega_{e}$. Although
the plots look like standard anticrossing behavior and decaying Rabi
oscillations, one should keep in mind that (1) the anticrossing occurs not at
the standard exciton-photon or phonon-photon resonance but at the nonlinear
parametric resonance, which is controlled by the nonlinear coupling strength
$\Omega_{R}^{(3)}$ and entangles three degrees of freedom; (2) the relaxation
rates of each individual subsystem enter the analytic expressions plotted in
Fig. 2 in a nontrivial way, as described above. This dependence will become
even more nontrivial at high temperature of the reservoir. The presented
solution provides the way to retrieve the analytic dependence of any
observable on the relaxation and coupling parameters and determine correctly
the criterion for observing the strong parametric coupling in frequency or
time domain. Two obvious examples for such observables are photon and phonon
emission spectra that are derived and plotted in Sec. VI, see Figs. 4 and 6.

(a)

(b)

Figure 2:  (a) Real parts of eigenstate frequencies of Eqs. (48), shifted by
the electron transition frequency $\omega_{e}$ and normalized by
$|\Omega_{R}^{(3)}|$, as a function of detuning from the parametric resonance
$\omega+\Omega-\omega_{e}$ normalized by $|\Omega_{R}^{(3)}|$. The relaxation
rates are $\mu_{\omega}=\mu_{\Omega}=0.3|\Omega_{R}^{(3)}|$ and
$\gamma=0.2|\Omega_{R}^{(3)}|$. (b) Occupation probabilities $|C_{001}|^{2}$
and $|C_{110}|^{2}$ from Eqs. (57) and (58) as a function of normalized time
$\tilde{\Omega}_{R}t$ for the same relaxation rates.

We also give the expressions for the occupation probabilities
$\overline{\left|C_{100}\right|^{2}}$, $\overline{\left|C_{010}\right|^{2}}$
and $\overline{\left|C_{000}\right|^{2}}$ that are valid under the condition
$\tilde{\Omega}_{R}\gg$ $\gamma_{ani}$:

$\overline{\left|C_{100}\right|^{2}}\approx\frac{\mu_{\omega}}{\mu_{\Omega}-\mu_{\omega}-\gamma}\left(e^{-\frac{\mu_{\Omega}+\mu_{\omega}+\gamma}{2}t}-e^{-\mu_{\Omega}t}\right)$
(62)
$\overline{\left|C_{010}\right|^{2}}\approx\frac{\mu_{\Omega}}{\mu_{\omega}-\mu_{\Omega}-\gamma}\left(e^{-\frac{\mu_{\Omega}+\mu_{\omega}+\gamma}{2}t}-e^{-\mu_{\omega}t}\right)$
(63)
$\overline{\left|C_{000}\right|^{2}}\approx\frac{\gamma\left(\gamma-\mu_{\Omega}-\mu_{\omega}\right)}{\gamma^{2}-\left(\mu_{\Omega}-\mu_{\omega}\right)^{2}}\left(1-e^{-\frac{\mu_{\Omega}+\mu_{\omega}+\gamma}{2}t}\right)-\left(\mu_{\Omega}\frac{1-e^{-\mu_{\omega}t}}{\mu_{\omega}-\mu_{\Omega}-\gamma}+\mu_{\omega}\frac{1-e^{-\mu_{\Omega}t}}{\mu_{\Omega}-\mu_{\omega}-\gamma}\right)$
(64)

It is easy to see that Eqs. (62)-(64) don’t have the divergence when
$\left[\pm\left(\mu_{\omega}-\mu_{\Omega}\right)-\gamma\right]\longrightarrow
0$; for example,

$\lim_{\left(\mu_{\Omega}-\mu_{\omega}-\gamma\right)\longrightarrow
0}\left[\frac{e^{-\frac{\mu_{\Omega}+\mu_{\omega}+\gamma}{2}t}-e^{-\mu_{\Omega}t}}{\mu_{\Omega}-\mu_{\omega}-\gamma}\right]=\frac{1}{2}te^{-\mu_{\Omega}t}$

and so on. In the limit $t\longrightarrow\infty$ Eqs. (62)-(64) lead to
$\overline{\left|C_{100}\right|^{2}}=\overline{\left|C_{010}\right|^{2}}=0$,
$\overline{\left|C_{000}\right|^{2}}=1$.

The detailed derivation of Eqs. (62)-(64) is given in Appendix A. It is also
shown there that in our case the quantities $\delta C_{100}$ and $\delta
C_{110}$ do not contribute to the observables and therefore can be omitted.

In tokman2021  we used a simplified model to analyze the fermion-photon-phonon
entanglement, in which all relaxation pathways to the ground state
$\left|0\right\rangle\left|0\right\rangle\left|0\right\rangle$ are assumed to
be “ direct”, which corresponds to taking all correlators equal to zero except
$D_{000,000}$. This approach is essentially the Weisskopf-Wigner method,
modified in order to conserve the norm of the state vector. It gives a correct
result for the decay rate $\gamma_{MIX}$ of the entangled state. At the same
time, including multistep decay pathways is of principal importance when
calculating and interpreting the emission spectra, as we will see in the next
section.

##  VI Emission spectra of photons and phonons from the parametric decay of
the electron excitation

###  VI.1  Derivation of the emission spectra from the solution of the
stochastic equation for the state vector: a general scheme

Consider for definiteness the EM radiation out of a cavity. Its power spectrum
received by the detector is given by scully1997 ; madsen2013

$P\left(\nu\right)=A\cdot S\left(\nu\right),$

where

$S\left(\nu\right)=\frac{1}{\pi}\mathrm{Re}\int_{0}^{\infty}d\tau
e^{i\nu\tau}\int_{0}^{\infty}dtK\left(t,\tau\right),$ (65)

$K=\left\langle\hat{c}^{\dagger}\left(t\right)\hat{c}\left(t+\tau\right)\right\rangle$,
$\left\langle\cdots\right\rangle$ is a quantum-mechanical averaging. The
coefficient A includes the $Q$-factor of a cavity, spatial structure of the
outgoing field, and the position and properties of the detector.

To calculate the power spectrum one needs to know the solution of the
Heisenberg-Langevin equations for the field operators $\hat{c}\left(t\right)$
and $\hat{c}^{\dagger}\left(t\right)$, then calculate the correlator, and
average it over the statistics of Langevin noise: $K\Rightarrow$
$\overline{\left\langle\hat{c}^{\dagger}\left(t\right)\hat{c}\left(t+\tau\right)\right\rangle}$.
However, as we already discussed, the Heisenberg-Langevin equations become
nonlinear in the strong-coupling regime. Therefore, it may be more convenient
to obtain the spectra from the solution of the stochastic equation Eq.(33) for
the state vector. The general procedure is as follows.

First, we need to transform the correlator

$K\left(t,\tau\right)=\left\langle\hat{c}^{\dagger}\left(t\right)\hat{c}\left(t+\tau\right)\right\rangle=\left\langle\Psi\left(0\right)\right|\hat{c}^{\dagger}\left(t\right)\hat{c}\left(t+\tau\right)\left|\Psi\left(0\right)\right\rangle$
(66)

to the Schrödinger picture without taking into account dissipation and
fluctuations. If $\hat{U}\left(t\right)$ is the unitary operator of evolution
of the system, one can write

$K=\left\langle\Psi\left(0\right)\right|\hat{U}^{\dagger}\left(t\right)\hat{c}^{\dagger}\,\hat{U}\left(t\right)\hat{U}^{\dagger}\left(t+\tau\right)\hat{c}\,\hat{U}\left(t+\tau\right)\left|\Psi\left(0\right)\right\rangle=\left\langle\hat{c}\,\Psi\left(t\right)\right|\hat{U}\left(t\right)\hat{U}^{\dagger}\left(t+\tau\right)\left|\hat{c}\,\Psi\left(t+\tau\right)\right\rangle,$
(67)

where $\hat{c}$ is the Schrödinger’s (constant) operator which we will treat
as an initial condition for the Heisenberg operator $\hat{c}\left(t\right)$ at
$t=0$. We will use the notation
$\hat{U}\left(t\right)\equiv\hat{U}_{t_{0}}\left(t^{\prime}\right)$, where we
indicate explicitly the initial moment of time $t_{0}$ and the duration of
evolution $t^{\prime}=t-t_{0}$. This will lead to the following replacements
in Eq. (67): $\hat{U}\left(t\right)\Longrightarrow\hat{U}_{0}\left(t\right)$,
$\hat{U}\left(t+\tau\right)\Longrightarrow\hat{U}_{0}\left(t+\tau\right)$.
Furthermore, we obviously have

$\hat{U}_{0}\left(t+\tau\right)=\hat{U}_{t}\left(\tau\right)\hat{U}_{0}\left(t\right)$
(68)

which gives

$K=\left\langle\hat{c}\,\Psi\left(t\right)\right|\hat{U}_{0}\left(t\right)\left(\hat{U}_{t}\left(\tau\right)\hat{U}_{0}\left(t\right)\right)^{\dagger}\left|\hat{c}\,\Psi\left(t+\tau\right)\right\rangle=\left\langle\hat{c}\,\Psi\left(t\right)\right|\hat{U}_{0}\left(t\right)\hat{U}_{0}^{\dagger}\left(t\right)\hat{U}_{t}^{\dagger}\left(\tau\right)\left|\hat{c}\,\Psi\left(t+\tau\right)\right\rangle.$

Taking into account
$\hat{U}_{0}\left(t^{\prime}\right)\hat{U}_{0}^{\dagger}\left(t^{\prime}\right)=1$,
we obtain

$K=\left\langle\hat{U}_{t}\left(\tau\right)\hat{c}\,\Psi\left(t\right)\right.\left|\hat{c}\,\Psi\left(t+\tau\right)\right\rangle.$
(69)

Second, introducing the notations

$\Psi_{\hat{C}}\left(t\right)=\hat{c}\,\Psi\left(t\right),\ \
\Phi\left(t,\tau\right)=\hat{U}_{t}\left(\tau\right)\Psi_{\hat{C}}\left(t\right),$
(70)

we arrive at

$K\left(t,\tau\right)=\left\langle\Phi\left(t,\tau\right)\right.\left|\Psi_{\hat{C}}\left(t+\tau\right)\right\rangle.$
(71)

Therefore, in order to calculate the correlator $K$ through the solution of
the equation for the state vector, one has to perform the following steps:

a) Find vector $\left|\Psi_{\hat{C}}\left(t+\tau\right)\right\rangle=$
$\hat{c}\Psi\left(t+\tau\right)$, where $\Psi\left(t+\tau\right)$ is the
solution for the state vector $\left|\Psi\right\rangle$ at the time interval
$\left[0,t+\tau\right]$ with initial condition
$\left|\Psi\left(0\right)\right\rangle$.

b) Find vector $\left|\Phi\left(t,\tau\right)\right\rangle$. To do that, one
has to solve for the state vector $\left|\Psi\right\rangle$ at the time
interval $\left[t,t+\tau\right]$ with initial condition
$\left|\Psi_{\hat{C}}\left(t\right)\right\rangle$. The vector
$\left|\Psi_{\hat{C}}\left(t\right)\right\rangle$ is the same as in part a),
but instead of the time interval $\left[0,t+\tau\right]$ one has to take the
time interval $\left[0,t\right]$.

To check Eq. (71) for consistency, we note that in the absence of dissipation
one can go back from this equation to the standard expression which follows
directly from the initial equation (66):

$K\left(t,\tau\right)=\left\langle\Psi\left(0\right)\right|e^{i\frac{\hat{H}}{\hbar}t}\hat{c}^{\dagger}e^{i\frac{\hat{H}}{\hbar}\tau}\hat{c}e^{-i\frac{\hat{H}}{\hbar}\left(t+\tau\right)}\left|\Psi\left(0\right)\right\rangle.$

If the dynamic system is open, then a complete closed system “the dynamic
system + reservoir” has its own unitary operator of evolution
$\hat{U}_{t_{0}}\left(t^{\prime}\right)$. Therefore, Eq. (71) should be valid
for a complete system as well which includes the reservoir variables. Now we
apply the Langevin method which assumes that the averaging over the statistics
of noise sources entering a stochastic equation (in this case Eq. (38)) is
equivalent to averaging over the reservoir variables. Therefore, we can solve
Eq. (38) and, following the above steps, find the functions
$\Psi_{\hat{C}}\left(t\right)$, $\Psi_{\hat{C}}\left(t+\tau\right)$ and
$\Phi\left(t,\tau\right)$ which are now dependent on the noise sources. Then
we substitute the latter two functions into Eq. (71) and perform averaging
over the noise statistics. As a result, we obtain

$K\left(t,\tau\right)=\overline{\left\langle\
\Phi\left(t,\tau\right)\right.\left|\Psi_{\hat{C}}\left(t+\tau\right)\right\rangle}.$
(72)

###  VI.2 Photon emission spectra for the parametric decay of an excited
electron

Here we apply the general recipe of calculating $K\left(t,\tau\right)$
formulated in the previous section to a particular example of the parametric
decay of an initially excited fermionic two-level system under strong coupling
to a parametric electron-photon-phonon resonance.

a) Use the expressions (56)-(59) to find the vector
$\left|\Psi\left(t\right)\right\rangle$.

b) Determine vectors $\left|\Psi_{\hat{C}}\left(t\right)\right\rangle$ and
$\left|\Psi_{\hat{C}}\left(t+\tau\right)\right\rangle$, resulting in

$\left|\Psi_{\hat{C}}\left(t\right)\right\rangle=\hat{c}\left|\Psi\left(t\right)\right\rangle=C_{110}\left(t\right)\left|1\right\rangle\left|0\right\rangle\left|0\right\rangle+C_{010}\left(t\right)\left|0\right\rangle\left|0\right\rangle\left|0\right\rangle,$
(73)
$\left|\Psi_{\hat{C}}\left(t+\tau\right)\right\rangle=\hat{c}\left|\Psi\left(t+\tau\right)\right\rangle=C_{110}\left(t+\tau\right)\left|1\right\rangle\left|0\right\rangle\left|0\right\rangle+C_{010}\left(t+\tau\right)\left|0\right\rangle\left|0\right\rangle\left|0\right\rangle.$
(74)

c) To determine the vector $\left|\Phi\left(t,\tau\right)\right\rangle$ we
will use the solution of Eqs. (47),(48) at the time interval
$\left[t,t+\tau\right]$ where the initial condition
$\left|\Psi_{\hat{C}}\left(t\right)\right\rangle$ is given by Eq. (73). In our
case Eq. (73) determines the initial value of the state vector; its subsequent
evolution is determined by a simple Eq. (47) for the amplitudes of states
$\left|1\right\rangle\left|0\right\rangle\left|0\right\rangle$ and
$\left|0\right\rangle\left|0\right\rangle\left|0\right\rangle$. As a result,
vector $\left|\Phi\left(t,\tau\right)\right\rangle$ is given by

$\displaystyle\left|\Phi\left(t,\tau\right)\right\rangle$ $\displaystyle=$
$\displaystyle
C_{100}^{\left(\Phi\right)}\left(t,\tau\right)\left|1\right\rangle\left|0\right\rangle\left|0\right\rangle+C_{000}^{\left(\Phi\right)}\left(t,\tau\right)\left|0\right\rangle\left|0\right\rangle\left|0\right\rangle$
(75) $\displaystyle=$
$\displaystyle\left(e^{-i\omega\tau-\gamma_{100}\tau}C_{110}\left(t\right)-\frac{i}{\hbar}\int_{t}^{t+\tau}\mathfrak{R}_{100}^{\left(\Phi\right)}\left(t,t^{\prime}\right)e^{-i\omega\left(\tau+t-t^{\prime}\right)-\gamma_{100}\left(\tau+t-t^{\prime}\right)}dt^{\prime}\right)\left|1\right\rangle\left|0\right\rangle\left|0\right\rangle$
$\displaystyle+\left(C_{010}\left(t\right)-\frac{i}{\hbar}\int_{t}^{t+\tau}\mathfrak{R}_{000}^{\left(\Phi\right)}\left(t,t^{\prime}\right)dt^{\prime}\right)\left|0\right\rangle\left|0\right\rangle\left|0\right\rangle,$

where functions $C_{110}\left(t\right)$ and $C_{010}\left(t\right)$ are
determined by Eqs. (56), (57).

The superscript $\left(\Phi\right)$ in the terms $\mathfrak{R}_{\alpha ni}$ in
Eq. (75) means that the correlators of these noise terms correspond to the
state vector $\left|\Phi\right\rangle$. The dependence on the initial time
moment $t$ of the evolution in $\mathfrak{R}_{\alpha
ni}^{\left(\Phi\right)}\left(t,t^{\prime}\right)$ takes into account that the
correlators of these random functions may depend on the value of $t$ as a
parameter, because complex amplitudes $C_{\alpha
ni}^{\left(\Phi\right)}\left(t,\tau\right)\equiv C_{\alpha
ni}^{\left(\Phi\right)}\left(t,t^{\prime}-t\right)$ depend on this parameter.

Next, we substitute the expressions for $C_{110}\left(t\right)$,
$C_{010}\left(t\right)$, $C_{110}\left(t+\tau\right)$ determined by Eqs. (56)
and (57), into Eqs. (74) and (75); after that we substitute Eqs. (74) and (75)
into Eq. (72). In the resulting expression we average over the noise
statistics taking into account that the noise sources are delta-correlated.
Omitting the terms that become zero after averaging, we obtain

$\displaystyle K\left(t,\tau\right)$ $\displaystyle=$ $\displaystyle
e^{-i\omega\tau-\left(\gamma_{100}+\frac{\gamma_{110}+\gamma_{001}}{2}\right)\tau-\left(\gamma_{110}+\gamma_{001}\right)t}\sin\left(\tilde{\Omega}_{R}t\right)\sin\left[\tilde{\Omega}_{R}\left(t+\tau\right)\right]$
(76)
$\displaystyle+e^{-i\omega\tau-\gamma_{010}\tau-2\gamma_{010}t}\int_{0}^{t}D_{010,010}\left(t^{\prime}\right)e^{2\gamma_{010}t^{\prime}}dt^{\prime}$
$\displaystyle+e^{-i\omega\left(t+\tau\right)-\gamma_{010}\left(t+\tau\right)-2\gamma_{010}t}\int_{t}^{t+\tau}\tilde{D}_{000,010}\left(t,t^{\prime}\right)e^{\left(i\omega+\gamma_{010}\right)t^{\prime}}dt^{\prime}.$

Here the quantity $D_{010,010}$ is determined by Eqs. (50):

$D_{010,010}=2\gamma_{010}\overline{\left|C_{110}\left(t^{\prime}\right)\right|^{2}}=\mu_{\Omega}\overline{\left|C_{110}\left(t^{\prime}\right)\right|^{2}}.$
(77)

The function $\tilde{D}_{000,010}\left(t,t^{\prime}\right)$ corresponds to the
following correlator:

$\overline{\mathfrak{R}_{000}^{\left(\Phi\right)\ast}\left(t,t^{\prime}\right)\mathfrak{R}_{010}\left(t^{\prime\prime}\right)}=\hbar^{2}\tilde{D}_{000,010}\left(t,t^{\prime}\right)\delta\left(t^{\prime}-t^{\prime\prime}\right).$
(78)

To calculate the value of $\tilde{D}_{000,010}\left(t,t^{\prime}\right)$ it is
not enough to have expressions (50) and (51) because it is determined by
correlations between the noise terms for different state vectors
$\left|\Phi\right\rangle$ and $\left|\Psi\right\rangle$, which correspond to
the solutions of the equation for the state vector with different initial
conditions. We need to use the expression for the noise source obtained in
Sec. IV.D. From Eqs. (52), (54), and (55) we obtain

$\tilde{D}_{000,010}\left(t,t^{\prime}\right)=\mu_{\Omega}\overline{C_{100}^{\left(\Phi\right)\ast}\left(t,t^{\prime}-t\right)C_{110}\left(t^{\prime}\right)}.$
(79)

Substituting here the appropriate term from Eq. (75) gives

$\tilde{D}_{000,010}\left(t,t^{\prime}\right)=e^{i\omega\left(t^{\prime}-t\right)-\gamma_{100}\left(t^{\prime}-t\right)-2\gamma_{010}t}\overline{C_{110}^{\ast}\left(t\right)C_{110}\left(t^{\prime}\right)}.$
(80)

Substituting Eqs. (77) and (80) into Eq. (76), we arrive at

$\displaystyle K\left(t,\tau\right)$ $\displaystyle=$ $\displaystyle
e^{-i\omega\tau-\left(\gamma_{100}+\frac{\Gamma}{2}\right)\tau-\Gamma
t}\sin\left(\tilde{\Omega}_{R}t\right)\sin\left[\tilde{\Omega}_{R}\left(t+\tau\right)\right]$
(81)
$\displaystyle+e^{-i\omega\tau-\gamma_{010}\tau}\mu_{\Omega}\left[\frac{2\tilde{\Omega}_{R}^{2}}{\gamma_{n}(4\tilde{\Omega}_{R}^{2}+\gamma_{n}^{2})}e^{-2\gamma_{010}t}-\frac{1}{2\gamma_{n}}e^{-\Gamma
t}\right.$
$\displaystyle\left.-\frac{1}{4\left(2i\tilde{\Omega}_{R}-\gamma_{n}\right)}e^{2i\tilde{\Omega}_{R}t-\Gamma
t}+\frac{1}{4\left(2i\tilde{\Omega}_{R}+\gamma_{n}\right)}e^{-2i\tilde{\Omega}_{R}t-\Gamma
t}\right]$ $\displaystyle+e^{-i\omega\tau-\gamma_{010}\tau-\Gamma
t}\frac{\mu_{\Omega}}{4}\left[\frac{e^{\left(-\gamma_{d}+i\tilde{\Omega}_{R}\right)\tau}-1}{-\gamma_{d}+i\tilde{\Omega}_{R}}\left(1-e^{2i\tilde{\Omega}_{R}t}\right)+c.c.\right].$

where we denoted $\Gamma=\gamma_{110}+\gamma_{001}$, $\gamma_{n}=$
$\Gamma-2\gamma_{010}$,
$\gamma_{d}=\gamma_{100}+\frac{\Gamma}{2}-\gamma_{010}$.

Now we have everything to determine the emission spectra given by Eq. (65).
Using the values of $\gamma_{100}=\frac{\mu_{\Omega}}{2}$ and
$\gamma_{010}=\frac{\mu_{\omega}}{2}$, we obtain

$S\left(\nu\right)=\frac{1}{\pi}\mathrm{Re}\int_{0}^{\infty}d\tau
e^{i\nu\tau}\int_{0}^{\infty}dtK\left(t,\tau\right)=S_{1}\left(\nu\right)+S_{2}\left(\nu\right)+S_{3}\left(\nu\right),$
(82)

where

$\displaystyle S_{1}\left(\nu\right)$ $\displaystyle=$
$\displaystyle\frac{2\tilde{\Omega}_{R}^{2}}{\pi\Gamma(4\tilde{\Omega}_{R}^{2}+\Gamma^{2})}\mathrm{Re}\frac{\Gamma+\frac{\mu_{\Omega}}{2}-i\left(\nu-\omega\right)}{[\gamma_{\mathrm{ac}}-i\left(\nu-\omega\right)]^{2}+\tilde{\Omega}_{R}^{2}}$
$\displaystyle S_{2}\left(\nu\right)$ $\displaystyle=$
$\displaystyle\frac{\tilde{\Omega}_{R}^{2}}{\pi\Gamma(4\tilde{\Omega}_{R}^{2}+\Gamma^{2})}\frac{\mu_{\Omega}}{\frac{\mu_{\omega}^{2}}{4}+\left(\nu-\omega\right)^{2}}$
(83) $\displaystyle S_{3}\left(\nu\right)$ $\displaystyle=$
$\displaystyle\frac{\mu_{\Omega}\tilde{\Omega}_{R}^{2}}{\pi\Gamma\left(\gamma_{d}^{2}+\tilde{\Omega}_{R}^{2}\right)(4\tilde{\Omega}_{R}^{2}+\Gamma^{2})}\times$
$\displaystyle\left\\{-\mathrm{Re}\frac{\Gamma_{d}\left[\gamma_{\mathrm{ac}}-i\left(\nu-\omega\right)\right]+2\tilde{\Omega}_{R}^{2}-\gamma_{d}\Gamma}{\left[\gamma_{\mathrm{ac}}-i\left(\nu-\omega\right)\right]^{2}+\tilde{\Omega}_{R}^{2}}+\frac{\Gamma_{d}\frac{\mu_{\omega}}{2}}{\frac{\mu_{\omega}^{2}}{4}+\left(\nu-\omega\right)^{2}}\right\\}$

The parameters $\Gamma$, $\gamma_{\mathrm{ac}}$, $\gamma_{d}$, and
$\Gamma_{d}$ are expressed through the relaxation rates of the electron,
photon, and phonon subsystems $\mu_{\omega}$, $\mu_{\Omega}$ and $\gamma$ as

$\displaystyle\Gamma$ $\displaystyle=$
$\displaystyle\gamma_{110}+\gamma_{001}=\frac{1}{2}\left(\mu_{\omega}+\mu_{\Omega}+\gamma\right),\
\gamma_{\mathrm{ac}}=\gamma_{100}+\frac{\Gamma}{2}=\frac{1}{4}\left(\mu_{\omega}+\gamma\right)+\frac{3}{4}\mu_{\Omega},$
$\displaystyle\gamma_{d}$ $\displaystyle=$
$\displaystyle\gamma_{100}+\frac{\Gamma}{2}-\gamma_{010}=\frac{1}{4}\left(\gamma-\mu_{\omega}\right)+\frac{3}{4}\mu_{\Omega},\
\Gamma_{d}=\Gamma+2\gamma_{d}=2\mu_{\Omega}+\gamma.$

The expression for the power spectrum $S\left(\nu\right)$ contains the terms
$S_{2}$ and $S_{3}$ which originate from the correlators $D_{010,010}$ and
$\tilde{D}_{000,010}$ respectively. The term $S_{3}\left(\nu\right)$ consists
of two terms which have the same spectral shapes as the functions
$S_{1}\left(\nu\right)$ and $S_{2}\left(\nu\right)$ respectively. Therefore,
including the term $S_{3}\left(\nu\right)$ in Eq. (82) leads only to
corrections to the amplitudes of the functions $S_{1,2}\left(\nu\right)$;
moreover, under the strong coupling conditions $\tilde{\Omega}_{R}\gg$
$\gamma_{ani}$ these corrections are small: of the order of
$\sim\frac{\mu_{\Omega}\left(\Gamma+\gamma_{\mathrm{ac}}\right)}{\tilde{\Omega}_{R}^{2}}$
for the function $S_{1}\left(\nu\right)$ and of the order of
$\sim\frac{\Gamma_{d}\mu_{\omega}}{\tilde{\Omega}_{R}^{2}}$ for the function
$S_{2}\left(\nu\right)$. For qualitative discussion we will neglect the
contribution of $S_{3}\left(\nu\right)$ and keep only the terms
$S_{1}\left(\nu\right)$ and $S_{2}\left(\nu\right)$, although all terms are
included in the spectra plotted in Fig. 4.

Figure 3:  Energy levels of
$\left|phonon\right\rangle\left|photon\right\rangle\left|electron\right\rangle$
states involved into the photon emission in the parametric decay of a single-
electron excitation in a coupled phonon-photon-electron system. Bold red
arrows indicate photon emission transitions with their peak frequencies
labeled. Wavy purple arrows indicate various relaxation pathways.  Figure 4:
Normalized photon emission spectra $\tilde{\Omega}_{R}^{2}S(\nu)$ as a
function of normalized detuning $\frac{\nu-\omega}{\tilde{\Omega}_{R}}$ from
the cavity mode frequency $\omega$, for two different values of the phonon
relaxation rate: $\mu_{\Omega}=0.3$ (red solid line) and $\mu_{\Omega}=0.02$
(blue dashed line). Other relaxation rates are $\mu_{\omega}=0.2$ and
$\gamma=0.1$. All relaxation constants are in units of $\tilde{\Omega}_{R}$.

Figure 3 indicates all transitions that give contributions to the photon
emission. The function $S_{1}\left(\nu\right)$ describes the emission spectrum
at the transition
$\left|1\right\rangle\left|1\right\rangle\left|0\right\rangle\rightarrow\left|1\right\rangle\left|0\right\rangle\left|0\right\rangle$,
which is split due to Rabi oscillations. The width of the peaks located at
frequencies $\nu=\omega\pm\tilde{\Omega}_{R}$ is equal to
$\gamma_{\mathrm{ac}}=\frac{\mu_{\Omega}}{2}+\frac{\Gamma}{2}$. Here
$\frac{\Gamma}{2}=\frac{1}{4}\left(\mu_{\omega}+\mu_{\Omega}+\gamma\right)=\gamma_{MIX}$
is the decay rate of the entangled state
$\left|MIX\right\rangle=A\left(t\right)\left|0\right\rangle\left|0\right\rangle\left|1\right\rangle+B\left(t\right)\left|1\right\rangle\left|1\right\rangle\left|0\right\rangle$
(see Eq. (59)), and $\frac{\mu_{\Omega}}{2}$ is the broadening of the state
$\left|1\right\rangle\left|0\right\rangle\left|0\right\rangle$ due to
relaxation. The spectrum given by $S_{1}\left(\nu\right)$ agrees with the one
obtained in tokman2021 .

The function $S_{2}\left(\nu\right)$ describes the emission due to a two-step
relaxation process described in section IV.C:
$\left|1\right\rangle\left|1\right\rangle\left|0\right\rangle\rightarrow\left|0\right\rangle\left|1\right\rangle\left|0\right\rangle\rightarrow\left|0\right\rangle\left|0\right\rangle\left|0\right\rangle$
. The photons are emitted at the transition
$\left|0\right\rangle\left|1\right\rangle\left|0\right\rangle\rightarrow\left|0\right\rangle\left|0\right\rangle\left|0\right\rangle$,
which is not affected by Rabi oscillations; see Fig. 3. Therefore, this
contribution has a standard Lorentzian shape of an emitter at frequency
$\omega$:

$S_{2}\left(\nu\right)\propto\frac{1}{\gamma_{010}^{2}+\left(\nu-\omega\right)^{2}}=\frac{1}{\frac{\mu_{\omega}^{2}}{4}+\left(\nu-\omega\right)^{2}}.$

In the strong-coupling regime $\tilde{\Omega}_{R}\gg$ $\gamma_{ani}$ the ratio
of the amplitude of this central peak at frequency $\nu=\omega$ to the
amplitudes of the split peaks at frequencies $\nu=\omega\pm\tilde{\Omega}_{R}$
is given by

$\frac{S_{2}\left(\omega\right)}{S_{1}\left(\omega\pm\tilde{\Omega}_{R}\right)}\approx\frac{\mu_{\Omega}\left(\mu_{\omega}+\gamma+3\mu_{\Omega}\right)}{\mu_{\omega}^{2}}.$
(84)

When $\mu_{\Omega}\rightarrow 0$, Eq. (84) gives
$\frac{S_{2}\left(\omega\right)}{S_{1}\left(\omega\pm\tilde{\Omega}_{R}\right)}\rightarrow
0$: indeed without phonon relaxation the state
$\left|0\right\rangle\left|1\right\rangle\left|0\right\rangle$ cannot be
populated from $\left|1\right\rangle\left|1\right\rangle\left|0\right\rangle$;
therefore, the two-step radiation channel is suppressed. In the opposite limit
of a fast phonon relaxation, when $\mu_{\Omega}\gg$ $\mu_{\omega},\gamma$, we
obtain
$\frac{S_{2}\left(\omega\right)}{S_{1}\left(\omega\pm\tilde{\Omega}_{R}\right)}\gg
1$, i.e., the side peaks are weaker and more broadened than the central peak.
The latter statement is true despite the large Rabi frequency
$\tilde{\Omega}_{R}\gg$ $\mu_{\Omega}$. Indeed, when $\tilde{\Omega}_{R}\gg$
$\mu_{\Omega}$ the interaction has the time to mix the states
$\left|1\right\rangle\left|1\right\rangle\left|0\right\rangle$ and
$\left|0\right\rangle\left|0\right\rangle\left|1\right\rangle$ before the
phonon relaxation kicks in. Nevertheless, if $\mu_{\Omega}\gg$
$\mu_{\omega},\gamma$ the phonon relaxation is able to transfer population to
state $\left|0\right\rangle\left|1\right\rangle\left|0\right\rangle$ faster
than the radiative transition
$\left|1\right\rangle\left|1\right\rangle\left|0\right\rangle\rightarrow\left|1\right\rangle\left|0\right\rangle\left|0\right\rangle$
corresponding to the Rabi-split spectrum.

This behavior is illustrated in Fig. 4 which shows photon emission spectra
given by Eqs. (82) and (83) as a function of frequency detuning $\nu-\omega$
from the cavity mode resonance, for two different values of the phonon
relaxation rate: $\mu_{\Omega}=1.5\mu_{\omega}$ (red solid line) and
$\mu_{\Omega}=0.1\mu_{\omega}$ (blue dashed line). The electron relaxation
rate is kept at $\gamma=0.1$ and its exact value is not important for the
overall shape of the spectra, although it affects absolute values and widths
of the peaks. All quantities are normalized by $\tilde{\Omega}_{R}$.

The relative magnitudes of the peaks and their widths depend sensitively on
different combinations of the relaxation rates
$\gamma,\mu_{\omega},\mu_{\Omega}$. The onset of the strong coupling regime in
the frequency domain is determined by the visibility of nonlinear Rabi
splitting between the side peaks in Fig. 4, i.e., the condition
$\tilde{\Omega}_{R}>\frac{\gamma_{ac}}{2}$. We point out again that the
relaxation rates of the individual subsystems enter the quantum dynamics in a
very nontrivial way, and one need to know all of them to evaluate the
feasibility of strong coupling in any particular system. The reverse is also
true: once the strong coupling regime is reached, measurements of the
photoluminescence spectra yield both the relaxation rates and the nonlinear
coupling strength in the system.

A more detailed discussion of the feasibility of strong coupling at the
nonlinear resonance in particular systems can be found in tokman2021 . Here we
only point out that in dielectric microcavities the photon relaxation rates
can be very low, in the $\mu$eV range, and the strong coupling threshold is
likely to be determined by relaxation of the electron or vibrational
transitions. In plasmonic nanocavities the photon relaxation rate can easily
be tens of meV and will likely dominate the strong coupling threshold. On the
other hand, the nonlinear coupling strength $\Omega_{R}^{(3)}$ is much higher
in plasmonic nanocavities because of greatly enhanced electric field
localization and electric field gradient. One can obtain the magnitude of
$\Omega_{R}^{(3)}$ of the order of 100 meV for the field localization in the
few nm range, which is now routinely demonstrated in plasmonic nanocavities.

###  VI.3 The phonon emission spectra

There is a complete symmetry for the two bosonic fields in the decay process
close to the parametric resonance $\omega_{e}=\omega+\Omega$. Therefore, we
can obtain the phonon emission spectrum from the expressions for the photon
emission spectrum Eqs. (82), (83), after replacing

$\omega\Longleftrightarrow\Omega,\
\mu_{\omega}\Longleftrightarrow\mu_{\Omega},\
\gamma_{010}\Longleftrightarrow\gamma_{100}.$

This results in

$S_{p}\left(\nu\right)=S_{1p}\left(\nu\right)+S_{2p}\left(\nu\right)+S_{3p}\left(\nu\right),$
(85)

where

$\displaystyle S_{1p}\left(\nu\right)$ $\displaystyle=$
$\displaystyle\frac{2\tilde{\Omega}_{R}^{2}}{\pi\Gamma(4\tilde{\Omega}_{R}^{2}+\Gamma^{2})}\mathrm{Re}\frac{\Gamma+\frac{\mu_{\Omega}}{2}-i\left(\nu-\Omega\right)}{[\tilde{\gamma}_{\mathrm{ac}}-i\left(\nu-\Omega\right)]^{2}+\tilde{\Omega}_{R}^{2}}$
$\displaystyle S_{2p}\left(\nu\right)$ $\displaystyle=$
$\displaystyle\frac{\tilde{\Omega}_{R}^{2}}{\pi\Gamma(4\tilde{\Omega}_{R}^{2}+\Gamma^{2})}\frac{\mu_{\omega}}{\frac{\mu_{\Omega}^{2}}{4}+\left(\nu-\Omega\right)^{2}}$
(86) $\displaystyle S_{3p}\left(\nu\right)$ $\displaystyle=$
$\displaystyle\frac{\mu_{\omega}\tilde{\Omega}_{R}^{2}}{\pi\Gamma\left(\tilde{\gamma}_{d}^{2}+\tilde{\Omega}_{R}^{2}\right)(4\tilde{\Omega}_{R}^{2}+\Gamma^{2})}\times$
$\displaystyle\left\\{-\mathrm{Re}\frac{\tilde{\Gamma}_{d}\left[\tilde{\gamma}_{\mathrm{ac}}-i\left(\nu-\Omega\right)\right]+2\tilde{\Omega}_{R}^{2}-\tilde{\gamma}_{d}\Gamma}{\left[\tilde{\gamma}_{\mathrm{ac}}-i\left(\nu-\Omega\right)\right]^{2}+\tilde{\Omega}_{R}^{2}}+\frac{\tilde{\Gamma}_{d}\frac{\mu_{\Omega}}{2}}{\frac{\mu_{\Omega}^{2}}{4}+\left(\nu-\Omega\right)^{2}}\right\\}$
$\displaystyle\tilde{\gamma}_{\mathrm{ac}}$ $\displaystyle=$
$\displaystyle\gamma_{010}+\frac{\Gamma}{2}=\frac{1}{4}\left(\mu_{\Omega}+\gamma\right)+\frac{3}{4}\mu_{\omega},\
\tilde{\gamma}_{d}=\gamma_{010}+\frac{\Gamma}{2}-\gamma_{100}=\frac{1}{4}\left(\gamma-\mu_{\Omega}\right)+\frac{3}{4}\mu_{\omega}$
$\displaystyle\ \tilde{\Gamma}_{d}$ $\displaystyle=$
$\displaystyle\Gamma+2\tilde{\gamma}_{d}=2\mu_{\omega}+\gamma.$

Similarly to the photon spectrum, the term $S_{3p}\left(\nu\right)$ is the sum
of two terms which have the same spectral shape as $S_{1p}\left(\nu\right)$
and $S_{2p}\left(\nu\right)$, but much smaller magnitudes if
$\tilde{\Omega}_{R}\gg$ $\gamma_{ani}$. Therefore we will again include only
the first two terms in qualitative discussion, but include all terms when
plotting the spectra.

Figure 5:  Energy levels of
$\left|phonon\right\rangle\left|photon\right\rangle\left|electron\right\rangle$
states involved into the phonon emission in the parametric decay of a single-
electron excitation in a coupled phonon-photon-electron system. Bold red
arrows indicate phonon emission transitions with their peak frequencies
labeled. Wavy purple arrows indicate various relaxation pathways.  Figure 6:
Normalized phonon emission spectra $\tilde{\Omega}_{R}^{2}S_{p}(\nu)$ as a
function of normalized detuning $\frac{\nu-\Omega}{\tilde{\Omega}_{R}}$ from
the vibrational mode frequency $\Omega$, for two different values of the
photon relaxation rate: $\mu_{\omega}=0.3$ (red solid line) and
$\mu_{\omega}=0.02$ (blue dashed line). Other relaxation rates are
$\mu_{\Omega}=0.2$ and $\gamma=0.1$. All relaxation constants are in units of
$\tilde{\Omega}_{R}$.

Figure 5 shows all transitions giving contributions to the phonon emission
spectrum, with their peak frequencies indicated.

The function $S_{1p}\left(\nu\right)$ describes the phonon emission spectrum
due to the transition
$\left|1\right\rangle\left|1\right\rangle\left|0\right\rangle\rightarrow\left|0\right\rangle\left|1\right\rangle\left|0\right\rangle$,
which demonstrates Rabi splitting. The width of the peaks centered at
frequencies $\nu=\Omega\pm\tilde{\Omega}_{R}$ is equal to
$\tilde{\gamma}_{\mathrm{ac}}=\frac{\mu_{\omega}}{2}+\frac{\Gamma}{2}$. The
function $S_{2p}\left(\nu\right)$ describes phonon emission due to a two-step
relaxation
$\left|1\right\rangle\left|1\right\rangle\left|0\right\rangle\rightarrow\left|1\right\rangle\left|0\right\rangle\left|0\right\rangle\rightarrow\left|0\right\rangle\left|0\right\rangle\left|0\right\rangle$.
The phonons are emitted at the second step, i.e. the transition
$\left|1\right\rangle\left|0\right\rangle\left|0\right\rangle\rightarrow\left|0\right\rangle\left|0\right\rangle\left|0\right\rangle$,
which is not affected by Rabi oscillations. Therefore, the spectrum due to
this contribution is a standard Lorenzian line, similarly to the case of
photons:

$S_{2}\left(\nu\right)\propto\frac{1}{\gamma_{100}^{2}+\left(\nu-\Omega\right)^{2}}=\frac{1}{\frac{\mu_{\Omega}^{2}}{4}+\left(\nu-\Omega\right)^{2}}.$

The ratio of the amplitude of the central peak at $\nu=\Omega$ to those of the
side peaks at frequencies $\nu=\Omega\pm\tilde{\Omega}_{R}$ at
$\tilde{\Omega}_{R}\gg$ $\gamma_{ani}$ is given by the expression equivalent
to Eq. (84) after substituting $\omega\Longleftrightarrow\Omega\ $and
$\mu_{\omega}\Longleftrightarrow\mu_{\Omega}$:

$\frac{S_{2p}\left(\Omega\right)}{S_{1p}\left(\Omega\pm\tilde{\Omega}_{R}\right)}\approx\frac{\mu_{\omega}\left(\mu_{\Omega}+\gamma+3\mu_{\omega}\right)}{\mu_{\Omega}^{2}}.$
(87)

The phonon emission spectra are plotted in Fig. 6 as a function of frequency
detuning $\nu-\Omega$ from the cavity mode resonance. This time we keep the
phonon relaxation rate fixed at $\mu_{\Omega}=0.2$ and plot the spectra for
two values of the photon relaxation rate, greater and smaller than
$\mu_{\Omega}$: $\mu_{\omega}=0.3$ (red solid line) and $\mu_{\omega}=0.02$
(blue dashed line). All quantities are normalized by $\tilde{\Omega}_{R}$. The
numbers are chosen to prove the point that the phonon and photon spectra are
symmetric with respect to replacement indicated in the beginning of this
section.

In experiment, measuring the ratios given by Eq. (84) and (87) allows one to
determine the relationships between all relaxation rates $\mu_{\Omega}$,
$\gamma$, and $\mu_{\omega}$. Indeed, one can obtain from Eqs. (84) and (87)
that

$\xi_{\Omega}x^{3}+2x^{2}-2x-\xi_{\omega}=0,\ y=\xi_{\Omega}x^{2}-3-x;$

where $x=\frac{\mu_{\Omega}}{\mu_{\omega}}$, $y=$
$\frac{\gamma}{\mu_{\omega}}$,
$\xi_{\omega}=\frac{S_{2}\left(\omega\right)}{S_{1}\left(\omega\pm\tilde{\Omega}_{R}\right)}$,
$\xi_{\Omega}=\frac{S_{2p}\left(\Omega\right)}{S_{1p}\left(\Omega\pm\tilde{\Omega}_{R}\right)}$.

##  VII Conclusions

We developed a universal model of three-wave nonlinear resonance which is
applicable to the systems with coupled electron, photon, and vibrational
degrees of freedom, for example in cavity QED with molecular quantum emitters
or quantum dots and in cavity optomechanics systems. We obtained the analytic
solution for the nonperturbative quantum dynamics of such systems in the
vicinity of the nonlinear resonance, taking into account dissipation and
fluctuations for all degrees of freedom in Markov approximation. We showed
that the strong coupling regime can be reached at the nonlinear resonance when
the nonlinear coupling strength exceeds a certain threshold which includes all
relaxation rates. When all three degrees of freedom are quantized, strong
coupling leads to the formation of tripartite entangled states. The presented
solution can be used to derive the explicit analytic expression for any
observable. As an example, we calculated photon and phonon emission spectra
which have a characteristic three-peak form once the strong coupling is
reached. We showed how the relative heights and widths of the peaks can be
used to extract information about all relaxation rates in the system and the
nonlinear coupling strength, or to establish the threshold for reaching the
strong coupling regime.

###### Acknowledgements.

This work has been supported in part by the Air Force Office for Scientific
Research Grant No. FA9550-21-1-0272, National Science Foundation Award No.
1936276, and Texas A&M University through STRP, X-grant and T3-grant programs.
M.T. and M.E. acknowledge the support by the Center of Excellence “Center of
Photonics” funded by The Ministry of Science and Higher Education of the
Russian Federation, Contract No. 075-15-2020-906.

##  Appendix A The derivation of Eqs. (62)-(64)

A number of correlators of the random functions in Eq. (59) is zero due to
Eqs. (50),(51):

$\overline{\delta C_{110}^{\ast}\left(t^{\prime}\right)\mathfrak{R}_{\alpha
mi}\left(t^{\prime\prime}\right)}=\overline{\delta
C_{001}^{\ast}\left(t^{\prime}\right)\mathfrak{R}_{\alpha
mi}\left(t^{\prime\prime}\right)}=0$ (88) $\overline{\delta
C_{001}^{\ast}\left(t^{\prime}\right)\delta
C_{001}\left(t^{\prime\prime}\right)}=\overline{\delta
C_{110}^{\ast}\left(t^{\prime}\right)\delta
C_{110}\left(t^{\prime\prime}\right)}=\overline{\delta
C_{001}^{\ast}\left(t^{\prime}\right)\delta
C_{110}\left(t^{\prime\prime}\right)}=0$ (89) $\displaystyle\overline{\delta
C_{001}^{\ast}\left(t^{\prime}\right)C_{000}\left(t^{\prime\prime}\right)}$
$\displaystyle=$ $\displaystyle\overline{\delta
C_{110}^{\ast}\left(t^{\prime}\right)C_{100}\left(t^{\prime\prime}\right)}=\overline{\delta
C_{001}^{\ast}\left(t^{\prime}\right)C_{010}\left(t^{\prime\prime}\right)}=\overline{\delta
C_{110}^{\ast}\left(t^{\prime}\right)C_{000}\left(t^{\prime\prime}\right)}$
(90) $\displaystyle=$ $\displaystyle\overline{\delta
C_{110}^{\ast}\left(t^{\prime}\right)C_{100}\left(t^{\prime\prime}\right)}=\overline{\delta
C_{110}^{\ast}\left(t^{\prime}\right)C_{010}\left(t^{\prime\prime}\right)}=0.$

Eqs. (88)-(90) ensure that the variables $\delta C_{100}$ and $\delta C_{110}$
cannot contribute to the values of any observables and therefore can be
omitted.

The other correlators are given by the equations that follow from Eqs. (56):

$\frac{d}{dt}\overline{C_{100}^{\ast}C_{010}}=-\left(\gamma_{100}+\gamma_{010}\right)\overline{C_{100}^{\ast}C_{010}}+D_{100,010},$
$\frac{d}{dt}\overline{C_{100}^{\ast}C_{000}}=-\gamma_{100}\overline{C_{100}^{\ast}C_{000}}+D_{100,000},$
$\frac{d}{dt}\overline{C_{010}^{\ast}C_{000}}=-\gamma_{010}\overline{C_{010}^{\ast}C_{000}}+D_{010,000},$
$\frac{d}{dt}\overline{\left|C_{000}\right|^{2}}=D_{000,000},$
$\frac{d}{dt}\overline{\left|C_{010}\right|^{2}}=-2\gamma_{010}\overline{\left|C_{010}\right|^{2}}+D_{010,010},$
$\frac{d}{dt}\overline{\left|C_{100}\right|^{2}}=-2\gamma_{100}\overline{\left|C_{100}\right|^{2}}+D_{100,100}.$

Using Eqs. (50),(51), we arrive at

$\displaystyle\frac{d}{dt}\overline{C_{100}^{\ast}C_{010}}$ $\displaystyle=$
$\displaystyle-\frac{\mu_{\omega}+\mu_{\Omega}}{2}\overline{C_{100}^{\ast}C_{010}}$
$\displaystyle\frac{d}{dt}\overline{C_{100}^{\ast}C_{000}}$ $\displaystyle=$
$\displaystyle-\frac{\mu_{\Omega}}{2}\overline{C_{100}^{\ast}C_{000}}+\mu_{\omega}\overline{C_{110}^{\ast}C_{010}}$
(91) $\displaystyle\frac{d}{dt}\overline{C_{010}^{\ast}C_{000}}$
$\displaystyle=$
$\displaystyle-\frac{\mu_{\omega}}{2}\overline{C_{010}^{\ast}C_{000}}+\mu_{\Omega}\overline{C_{110}^{\ast}C_{100}}$

$\displaystyle\frac{d}{dt}\overline{\left|C_{000}\right|^{2}}$
$\displaystyle=$
$\displaystyle\gamma\overline{\left|C_{001}\right|^{2}}+\mu_{\omega}\overline{\left|C_{010}\right|^{2}}+\mu_{\Omega}\overline{\left|C_{100}\right|^{2}}$
$\displaystyle\frac{d}{dt}\overline{\left|C_{010}\right|^{2}}$
$\displaystyle=$
$\displaystyle-\mu_{\omega}\overline{\left|C_{010}\right|^{2}}+\mu_{\Omega}\overline{\left|C_{110}\right|^{2}}$
(92) $\displaystyle\frac{d}{dt}\overline{\left|C_{100}\right|^{2}}$
$\displaystyle=$
$\displaystyle-\mu_{\Omega}\overline{\left|C_{100}\right|^{2}}+\mu_{\omega}\overline{\left|C_{110}\right|^{2}}$

Taking into account Eqs. (90), one can obtain that
$\overline{C_{110}^{\ast}C_{010}}=\overline{C_{110}^{\ast}C_{100}}$ $=0$ in
Eqs. (91); as a result, for our initial conditions Eqs. (91) yield

$\overline{C_{100}^{\ast}C_{010}}=\overline{C_{100}^{\ast}C_{000}}=\overline{C_{010}^{\ast}C_{000}}=0.$
(93)

The last two equations in Eqs. (92) give

$\overline{\left|C_{010}\right|^{2}}=\mu_{\Omega}e^{-\mu_{\omega}t}\int_{0}^{t}e^{\mu_{\omega}\tau}\overline{\left|C_{110}\right|^{2}}d\tau,\
\
\overline{\left|C_{100}\right|^{2}}=\mu_{\omega}e^{-\mu_{\Omega}t}\int_{0}^{t}e^{\mu_{\Omega}\tau}\overline{\left|C_{110}\right|^{2}}d\tau.$

Substituting here the function from Eq. (57) and taking into account Eqs.
(89),(90) results in

$\overline{\left|C_{100}\right|^{2}}=\mu_{\omega}e^{-\mu_{\Omega}t}\int_{0}^{t}e^{\frac{\mu_{\Omega}-\mu_{\omega}-\gamma}{2}\tau}\sin^{2}\left(\tilde{\Omega}_{R}\tau\right)d\tau,\
\
\overline{\left|C_{010}\right|^{2}}=\mu_{\Omega}e^{-\mu_{\omega}t}\int_{0}^{t}e^{\frac{\mu_{\omega}-\mu_{\Omega}-\gamma}{2}\tau}\sin^{2}\left(\tilde{\Omega}_{R}\tau\right)d\tau;$

For further integration we use the limit $\frac{\gamma_{\alpha
ni}}{\tilde{\Omega}_{R}}\ll 1$, leading to

$\overline{\left|C_{100}\right|^{2}}\approx\frac{\mu_{\omega}}{\mu_{\Omega}-\mu_{\omega}-\gamma}\left(e^{-\frac{\mu_{\Omega}+\mu_{\omega}+\gamma}{2}t}-e^{-\mu_{\Omega}t}\right)$
(94)
$\overline{\left|C_{010}\right|^{2}}\approx\frac{\mu_{\Omega}}{\mu_{\omega}-\mu_{\Omega}-\gamma}\left(e^{-\frac{\mu_{\Omega}+\mu_{\omega}+\gamma}{2}t}-e^{-\mu_{\omega}t}\right)$
(95)

Note that Eqs. (94),(95) do not contain any divergence when
$\left[\pm\left(\mu_{\omega}-\mu_{\Omega}\right)-\gamma\right]\longrightarrow
0$. Indeed,

$\lim_{\left(\mu_{\Omega}-\mu_{\omega}-\gamma\right)\longrightarrow
0}\left[\frac{e^{-\frac{\mu_{\Omega}+\mu_{\omega}+\gamma}{2}t}-e^{-\mu_{\Omega}t}}{\mu_{\Omega}-\mu_{\omega}-\gamma}\right]=\frac{1}{2}te^{-\mu_{\Omega}t},\
\lim_{\left(\mu_{\omega}-\mu_{\Omega}-\gamma\right)\longrightarrow
0}\left[\frac{e^{-\frac{\mu_{\Omega}+\mu_{\omega}+\gamma}{2}t}-e^{-\mu_{\omega}t}}{\mu_{\omega}-\mu_{\Omega}-\gamma}\right]=\frac{1}{2}te^{-\mu_{\omega}t}.$

Now we return to the first of Eqs. (92), which yields

$\overline{\left|C_{000}\right|^{2}}=\int_{0}^{t}\left(\gamma\overline{\left|C_{001}\right|^{2}}+\mu_{\omega}\overline{\left|C_{010}\right|^{2}}+\mu_{\Omega}\overline{\left|C_{100}\right|^{2}}\right)d\tau.$

We substitute Eqs. (94) and (95) into the second and third terms in the
integrand and substitute the expression $\overline{\left|C_{001}\right|^{2}}$
which follows from Eqs. (57) into the first term in the integrand. Neglecting
the small terms $\propto$
$\frac{\gamma_{110}-\gamma_{001}}{\tilde{\Omega}_{R}}$and
$\frac{\gamma_{\alpha ni}}{\tilde{\Omega}_{R}}$ the integration results in

$\overline{\left|C_{000}\right|^{2}}=\frac{\gamma\left(\gamma-\mu_{\Omega}-\mu_{\omega}\right)}{\gamma^{2}-\left(\mu_{\Omega}-\mu_{\omega}\right)^{2}}\left(1-e^{-\frac{\mu_{\Omega}+\mu_{\omega}+\gamma}{2}t}\right)-\left(\mu_{\Omega}\frac{1-e^{-\mu_{\omega}t}}{\mu_{\omega}-\mu_{\Omega}-\gamma}+\mu_{\omega}\frac{1-e^{-\mu_{\Omega}t}}{\mu_{\Omega}-\mu_{\omega}-\gamma}\right)$
(96)

## References

  * (1) M. K. Schmidt, R. Esteban, A. Gonzalez-Tudela, G. Giedke, and J. Aizpurua, Quantum Mechanical Description of Raman Scattering from Molecules in Plasmonic Cavities, ACS Nano 10, 6291-6298 (2016).
  * (2) J. del Pino, J. Feist, and F. J. Garcia-Vidal, Signatures of Vibrational Strong Coupling in Raman Scattering, J. Phys. Chem. 119, 29132-29137 (2015).
  * (3) J. del Pino, F. J. Garcia-Vidal, and J. Feist, Exploiting vibrational strong coupling to make an optical parametric oscillator out of a Raman laser, Phys. Rev. Lett. 117, 277401 (2016).
  * (4) S. Ramelow, A. Farsi, Z. Vernon, S. Clemmen, X. Ji, J. E. Sipe, M. Liscidini, M. Lipson, and A. L. Gaeta, Strong Nonlinear Coupling in a Si3N4 Ring Resonator, Phys. Rev. Lett. 122, 153906 (2019).
  * (5) T. Neuman, J. Aizpurua and R. Esteban, Quantum theory of surface-enhanced resonant Raman scattering (SERRS) of molecules in strongly coupled plasmon-exciton systems, Nanophotonics 9, 295-308 (2020).
  * (6) F. Ge, X. Han, and J. Xu, Strongly coupled systems for nonlinear optics, Laser Photonics Rev. 15, 2000514 (2021).
  * (7) K. Wang, M. Seidel, K. Nagarajan, T. Chervy, C. Genet, and T. Ebbesen, Large optical nonlinearity enhancement under electronic strong coupling, Nat. Commun. 12, 1486 (2021).
  * (8) J. Flick, N. Rivera and P. Narang, Strong light-matter coupling in quantum chemistry and quantum photonics, Nanophotonics 7, 1479-1501 (2018).
  * (9) D. Wang, J. Flick, and P. Narang, Light-matter interaction of a molecule in a dissipative cavity from first principles, J. Chem. Phys. 154, 104109 (2021).
  * (10) A. V. Zasedatelev, A. V. Baranikov, D. Sannikov, D. Urbonas, F. Scafirimuto, V. Yu. Shishkov, E. S. Andrianov, Y. E. Lozovik, U. Scherf, T. Stoferle, R. F. Mahrt and P. G. Lagoudakis, Single-photon nonlinearity at room temperature, Nature 597, 493 (2021).
  * (11) A. Pscherer, M.Meierhofer, D. Wang, Hr. Kelkar, D. Martin-Cano, T. Utikal ,S. Gotzinger, and V. Sandoghdar, Single-Molecule Vacuum Rabi Splitting: Four-Wave Mixing and Optical Switching at the Single-Photon Level, Phys. Rev. Lett. 127, 133603 (2021).
  * (12) S. Hughes, A. Settineri, S. Savasta, and F. Nori, Resonant Raman scattering of single molecules under simultaneous strong cavity coupling and ultrastrong optomechanical coupling in plasmonic resonators: Phonon-dressed polaritons, Phys. Rev. B 104, 045431 (2021).
  * (13) L.D. Landau, E.M. Lifshitz, Electrodynamics of continuous media (Elsevier, Oxford, 1984).
  * (14) F. Benz, M. K. Schmidt, A. Dreismann, R. Chikkaraddy, Y. Zhang, A. Demetriadou, C. Carnegie, H. Ohadi, B. de Nijs, R. Esteban, J. Aizpurua, and J. J. Baumberg, Single-molecule optomechanics in picocavities, Science 354, 726 (2016).
  * (15) K.-D. Park, E. A. Muller, V. Kravtsov, P. M. Sass, J. Dreyer, J. M. Atkin, and M. B. Raschke, Variable-temperature tip-enhanced Raman spectroscopy of single-molecule fluctuations and dynamics, Nano Lett. 16, 479 (2016).
  * (16) R. Chikkaraddy, B. de Nijs, F. Benz, S. J. Barrow, O. A. Scherman, E. Rosta, A. Demetriadou, P. Fox, O. Hess, and J. J. Baumberg, Single-molecule strong coupling at room temperature in plasmonic nanocavities, Nature 535, 127 (2016).
  * (17) M. Tokman, M. Erukhimova, Y. Wang, Q. Chen, and A. Belyanin, Generation and dynamics of entangled fermion-photon-phonon states in nanocavities, Nanophotonics 10, 491 (2021).
  * (18) M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Cavity optomechanics, Rev. Mod. Phys. 86, 1391 (2014).
  * (19) P. Meystre, A short walk through quantum optomechanics, Ann. Phys. 525, 215 (2013).
  * (20) P. Roelli, C. Galland, N. Piro and T. J. Kippenberg. Molecular cavity optomechanics as a theory of plasmon-enhanced Raman scattering, Nature Nano. 11, 164 (2016).
  * (21) J.-M. Pirkkalainen, S.U. Cho, F. Massel, J. Tuorila, T.T. Heikkila, P.J. Hakonen, and M.A. Sillanpaa, Cavity optomechanics mediated by a quantum two-level system, Nat. Comm. 6, 6981 (2015).
  * (22) Y. Chu, P. Kharel, W. H. Renninger, L. D. Burkhart, L. Frunzio, P. T. Rakich, R. J. Schoelkopf, Quantum acoustics with superconducting qubits, Science 358, 199 (2017).
  * (23) S. Hong, R. Riedinger, I. Marinkovic, A. Wallucks, S. G. Hofer, R. A. Norte, M. Aspelmeyer, and S. Groblacher, Hanbury Brown and Twiss interferometry of single phonons from an optomechanical resonator, Science 358, 203 (2017).
  * (24) P. Arrangoiz-Arriola, E. A. Wollack, Z. Wang, M. Pechal, W. Jiang, T. P. McKenna, J. D. Witmer, R. Van Laer, and A.H. Safavi-Naeini, Resolving the energy levels of a nanomechanical oscillator, Nature 571, 537 (2019).
  * (25) Q. Liao, Y. Ye, P. Lin, N. Zhou, and W. Nie, Tripartite entanglement in an atom-cavity-optomechanical system, Int. J. Theor. Phys. 57, 1319 (2018).
  * (26) K. J. Satzinger, Y. P. Zhong, H. Chang, et al., Quantum control of surface acoustic-wave phonons, Nature 563, 661 (2018).
  * (27) A. Bienfait, Y. P. Zhong, H.-S. Chang, M.-H. Chou, C. R. Conner, E. Dumur, J. Grebel, G. A. Peairs, R. G. Povey, K. J. Satzinger, and A. N. Cleland, Quantum Erasure Using Entangled Surface Acoustic Phonons, Phys. Rev. X 10, 021055 (2020).
  * (28) U. Delic, M. Reisenbauer, K. Dare, D. Grass, V. Vuletic, N. Kiesel, and M. Aspelmeyer, Cooling of a levitated nanoparticle to the motional quantum ground state, Science 367, 892 (2020).
  * (29) I. Wilson-Rae and A. Imamoglu, Quantum dot cavity-QED in the presence of strong electron-phonon interactions, Phys. Rev. B 65, 235311 (2002).
  * (30) S. Weiler, A. Ulhaq, S. M. Ulrich, D. Richter, M. Jetter, P. Michler, C. Roy, and S. Hughes, Phonon-assisted incoherent excitation of a quantum dot and its emission properties, Phys. Rev. B 86, 241304(R) (2012).
  * (31) K, Roy-Choudhury and S. Hughes, Quantum theory of the emission spectrum from quantum dots coupled to structured photonic reservoirs and acoustic phonons, Phys. Rev. B 92, 205406 (2015).
  * (32) E. V. Denning, M. Bundgaard-Nielsen, and J. Mork, Optical signatures of electron-phonon decoupling due to strong light-matter interactions, Phys. Rev. B 102, 235303 (2020).
  * (33) P. Zoller and C. W. Gardiner, Quantum Noise in Quantum Optics: the Stochastic Schrödinger Equation. Lecture Notes for the Les Houches Summer School LXIII on Quantum Fluctuations in July 1995, Elsevier Science Publishers B.V. 1997, edited by E. Giacobino and S. Reynaud; arXiv:quant-ph/9702030v1.
  * (34) M. O. Scully and M. S. Zubairy, Quantum Optics (Cambridge University Press, 1997).
  * (35) M. B. Plenio and P. L. Knight, The quantum-jump approach to dissipative dynamics in quantum optics, Rev. Mod. Phys. 70, 101 (1998)
  * (36) N. Gisin and I. C. Percival, The quantum-state diffusion model applied to open systems, Journal of Physics A 25 , 5677 (1992).
  * (37) L. Diosi, N. Gisin, and W. T. Strunz, Non-Markovian quantum state diffusion, Phys. Rev. A 58, 1699-1712 (1998).
  * (38) C. Cohen-Tannoudji, B. Zambon and E. Arimondo, Quantum-jump approach to dissipative processes: application to amplification without inversion, J. Opt. Soc. Am. B 10, 2107-2120 (1993).
  * (39) K. Molmer, Y. Castin and J. Dalibard, Monte Carlo wave-function method in quantum optics, J. Opt. Soc. Am. B 10, 524-538 (1993).
  * (40) N. Gisin and I. C. Percival, Wave-function approach to dissipative processes: are there quantum jumps? Phys. Lett. A 167, 315-318 (1992).
  * (41) H. De Raedt, F. Jin, M. I. Katsnelson, and K. Michielsen, Relaxation, thermalization, and Markovian dynamics of two spins coupled to a spin bath, Phys. Rev. E 96, 053306 (2017).
  * (42) F. Nathan and M. S. Rudner, Universal Lindblad equation for open quantum systems, Phys. Rev. B 102, 115109 (2020).
  * (43) Q. Chen, Y. Wang, S. Almutairi, M. Erukhimova, M. Tokman, and A. Belyanin. Dynamics and control of entangled electron-photon states in nanophotonic systems with time-variable parameters, Phys. Rev. A 103, 013708 (2021).
  * (44) M. Tokman, Q. Chen, M. Erukhimova, Y. Wang, and A. Belyanin, Quantum dynamics of open many-qubit systems strongly coupled to a quantized electromagnetic field in dissipative cavities, https://arxiv.org/abs/2105.14674.
  * (45) M. Tokman, Y. Wang, I. Oladyshkin, A. R. Kutayiah, and A. Belyanin, Laser-driven parametric instability and generation of entangled photon-plasmon states in graphene, Phys. Rev. B 93, 235422 (2016).
  * (46) M. Tokman, X. Yao, and A. Belyanin, Generation of entangled photons in graphene in a strong magnetic field, Phys. Rev. Lett. 110, 077404 (2013).
  * (47) M. D. Tokman, M. A. Erukhimova, and V. V. Vdovin, The features of a quantum description of radiation in an optically dense medium, Ann. Phys. 360, 571 (2015).
  * (48) M. Tokman, Z. Long, S. AlMutairi, Y. Wang, M. Belkin, and A. Belyanin, Enhancement of the spontaneous emission in subwavelength quasi-two-dimensional waveguides and resonators, Phys. Rev. A. 97, 043801 (2018).
  * (49) F. Herrera, and F. C. Spano. Absorption and photoluminescence in organic cavity QED, Phys. Rev. A 95, 053867 (2017).
  * (50) E. T. Jaynes and F. W. Cummings, Comparison of quantum and semiclassical radiation theories with application to the beam maser, Proc. IEEE 51, 89 (1963).
  * (51) C. Couteau. Spontaneous parametric down-conversion, Cont. Phys. 59, 291-304 (2018).
  * (52) M. Tokman, Y. Wang, Q. Chen, L. Shterengas, and A. Belyanin, Generation of entangled photons via parametric down-conversion of lasing modes in semiconductor lasers and waveguides, https://arxiv.org/abs/2108.03528.
  * (53) M. D. Kostin, On the Schrödinger-Langevin Equation, J. Chem. Phys. 57, 3589 (1972).
  * (54) R. Katz, P.B. Gossiaux. The Schrödinger–Langevin equation with and without thermal fluctuations, Annals of Physics 368, 267-295 (2016).
  * (55) D.F. Walls and G. J. Milburn, Quantum Optics (Springer, Berlin, 1994).
  * (56) L.D. Landau, E.M. Lifshitz, Statistical Physics, Part 1 (Pergamon, Oxford, 1965)
  * (57) K. H. Madsen and P. Lodahl, Quantitative analysis of quantum dot dynamics and emission spectra in cavity quantum electrodynamics, New J. Phys. 15, 025013 (2013).