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

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.

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 .

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}$.

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 .

$$\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 .

$$\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.

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.

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)

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.

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$.