Introduction

Electromagnetically induced transparency discovered in quantum optics has long been an important effect in physics (see, e.g. ref. 1, for a review). This phenomenon of absorption cancelation is interpreted as the appearance of dark state or coherent population trapping. In addition to the atomic systems, dark state has also been observed in a number of solid-state systems including quantum dots2, 3, nitrogen-vacancy center4 and silicon-vacancy center in diamond5, 6. In fact, dark state can have different applications in physics. The atomic clocks based on coherent population trapping7,8,9,10 make the high-precision time estimation possible using the chip-scale and low-power devices. The state transfer can be done with the adiabatic passage of the dark state11,12,13. Operations on the quantum states like squeezing14 or decay suppression15, 16 can also be conducted with the help of dark state. Moreover, dark state can have important applications to the slow light17 and photocell18.

For the Λ-type three-level system, one of the advantages of the dark state is that it is composed of two lowest states without dipole transition between them. Within the framework of rotating-wave approximation (RWA) for the interaction between the system and the environment, the dark state is dissipation-free at a low enough temperature. For instance, the dark state is not influenced by the spontaneous emission1. Studies in two-level systems have shown that the counter-rotating terms can change the ground state19,20,21, but there were few studies regarding the influence of the counter-rotating terms on the dark state22. When the coupling between the system and the environment becomes strong, the counter-rotating terms cannot be neglected. Thus, interesting phenomena with the quantum dynamics of the dark state are expected even at zero temperature, because now the ground state within the framework of RWA is no longer the ground state of the system when including the counter-rotating terms.

In this article, we study the quantum dynamics of the dark state beyond the RWA, where the Λ-type three-level system couples to two bosonic baths at zero temperature and the couplings between the system and the two baths contain both rotating and counter-rotating terms. We derive a non-Markovoian quantum Bloch equation for the dark state using a quantum Langevin approach. In contrast to the dark state within the RWA, leakage of the dark state occurs due to the counter-rotating terms in the system-bath interaction, revealing the breakdown of the dissipation-free dark state at the zero temperature. To suppress the leakage, we apply a leakage elimination operator to the system, which plays the role of keeping the upper level of the Λ-type three-level system empty. Indeed, the leakage of the dark state can be much reduced when applying the elimination operator, as shown in our numerical results. This study provides a method to restore the quantum coherence of the dark state with counter-rotating dissipative channels.

Results

The model Hamiltonian

We study a Λ-type three-level system driven by two pump fields [see Fig. 1(a)]. The Hamiltonian of this three-level system is given by (setting ħ = 1)

$$\begin{array}{rcl}{H}_{{\rm{sys}}} & = & {\omega }_{1}\mathrm{|1}\rangle \langle \mathrm{1|}+{\omega }_{2}\mathrm{|2}\rangle \langle \mathrm{2|}+{\omega }_{3}\mathrm{|3}\rangle \langle \mathrm{3|,}\\ & & +({{\rm{\Omega }}}_{1}{e}^{-i{\omega }_{a}t}\mathrm{|3}\rangle \langle \mathrm{1|}+{{\rm{\Omega }}}_{2}{e}^{-i{\omega }_{b}t}\mathrm{|3}\rangle \langle \mathrm{2|}+{\rm{H}}{\rm{.c}}\mathrm{.),}\end{array}$$
(1)

where ω 1, ω 2, and ω 3 are the three energy levels of the system and the frequencies of the two pumping fields are ω a and ω b , respectively, with Ω1 and Ω2 characterizing their coupling strengths to the three-level system. To focus on the effect of the counter-rotating dissipative channels, we take Ω1 and Ω2 to be independent of time. Also, the effects like dephasing channels or non-adiabatic transitions are not included here because they were well studied13. Here we consider the resonant pumping case with \({\omega }_{a}={\omega }_{3}-{\omega }_{1}\) and \({\omega }_{b}={\omega }_{3}-{\omega }_{2}\). This Λ-type three-level system has a dark state1,

$$|D(t)\rangle =\frac{{{\rm{\Omega }}}_{2}}{{\rm{\Omega }}}{e}^{-i{\omega }_{1}t}|1\rangle -\frac{{{\rm{\Omega }}}_{1}}{{\rm{\Omega }}}{e}^{-i{\omega }_{2}t}|2\rangle ,$$
(2)

with \({\rm{\Omega }}=\sqrt{|{{\rm{\Omega }}}_{1}{|}^{2}+|{{\rm{\Omega }}}_{2}{|}^{2}}\), which is a solution of the Schrödinger equation, \(\frac{\partial }{\partial t}|D(t)\rangle =-i{H}_{{\rm{s}}ys}|D(t)\rangle \). When the three-level system is in this state, it remains in the subspace spanned by \(\{\mathrm{|1}\rangle \mathrm{,|2}\rangle \}\), even in the presence of the two pumping fields.

Figure 1
figure 1

(a) A Λ-type three-level system driven by two fields with frequencies ω a and ω b , respectively. The field on the left (right), which drives the transition between |1〉 (|2〉) and |3〉, has a coupling strength Ω1 (Ω2) with this transition. The three-level system is also coupled to two bosonic baths, with the coupling strengths characterized by \({{\rm{\Gamma }}}_{a}\) and \({{\rm{\Gamma }}}_{b}\). (b) Schematic illustration of the applied control pulses, where τ is the period of the pulses, Δ is the duration of each pulse, and h is the strength of the pulse.

Full size image

To study the effect of the environments on the dark state, we use two independent bosonic baths modeled by \({H}_{{\rm{bath}}}={\sum }_{k}{\omega }_{a,k}{a}_{k}^{\dagger }{a}_{k}+{\sum }_{k}{\omega }_{b,k}{b}_{k}^{\dagger }{b}_{k}\). The interaction between the system and the two baths is

$$\begin{array}{rcl}{H}_{{\rm{int}}} & = & (|1\rangle \langle 3|+|3\rangle \langle 1|)\sum _{k}({g}_{a,k}{a}_{k}+{g}_{a,k}^{\ast }{a}_{k}^{\dagger })\\ & & +\,(|2\rangle =\langle 3|+|3\rangle \langle 2|)\sum _{k}({g}_{b,k}{b}_{k}+{g}_{b,k}^{\ast }{b}_{k}^{\dagger }),\end{array}$$
(3)

where \({a}_{k}\) and \({b}_{k}\) are annihilation operators of the bosonic modes of the two baths. Note that both rotating and counter-rotating terms are included in this interaction Hamiltonian. The dipole transition between \(\mathrm{|1}\rangle \) and \(\mathrm{|2}\rangle \) is forbidden in the considered Λ-type three-level system, so we omit this channel here. Under the RWA (i.e., only the rotating terms are considered), Eq. (3) is reduced to

$${H}_{{\rm{RWA}}}=\sum _{k}({g}_{a,k}{a}_{k}|3\rangle \langle 1|+{g}_{b,k}{b}_{k}|3\rangle \langle 2|)+{\rm{H}}{\rm{.c}}{\rm{.}}$$
(4)

Let the two bosonic baths be in the vacuum state \(\mathrm{|0}\rangle \equiv \mathrm{|0}{\rangle }_{a}\otimes \mathrm{|0}{\rangle }_{b}\), which corresponds to the zero temperature for each bath. It is easy to check that \(({H}_{{\rm{RWA}}}+{H}_{{\rm{bath}}})|D(t)\rangle \otimes \mathrm{|0}\rangle =0\). Thus,

$$\frac{\partial }{\partial t}|D(t)\rangle \otimes |0\rangle =-i({H}_{{\rm{sys}}}+{H}_{{\rm{RWA}}}+{H}_{{\rm{bath}}})|D(t)\rangle \otimes |0\rangle ,$$
(5)

i.e., the dark state can persist even when the three-level system couples to the two baths. However, the counter-rotating terms in Eq. (3) can make the dark state transition to the state \(\mathrm{|3}\rangle \), so \(|D(t)\rangle \otimes \mathrm{|0}\rangle \) is not a solution of the Schrödinger equation, \(\frac{\partial }{\partial t}|{\rm{\Psi }}\rangle =-i{H}_{{\rm{t}}ot}|{\rm{\Psi }}\rangle \), where \({H}_{{\rm{tot}}}={H}_{{\rm{sys}}}+{H}_{{\rm{int}}}+{H}_{{\rm{bath}}}\) is the total Hamiltonian of the system.

Below we show how the two baths affect the quantum coherence of the dark state when the interaction Hamiltonian is H int, so as to reveal the effect of the counter-rotating terms in the interaction Hamiltonian. Then, we present a method to restore the quantum coherence of the dark state by applying a leakage elimination operator to the system.

Non-Markovian quantum Bloch equation

The reduced density operator of the system can be written as

$${\rho }_{{\rm{sys}}}(t)=\sum _{i,j}{\rho }_{ij}^{(\mathrm{sys})}(t)|i\rangle \langle j|,$$
(6)

where \(i,j=1,2,3\), and the reduced density matrix elements are given by

$${\rho }_{ij}^{(\mathrm{sys})}(t)=\langle i|{\rm{T}}{r}_{{\rm{e}}nv}\rho (t)|j\rangle \mathrm{.}$$
(7)

Here \(\rho (t)\) is the density operator of the total system and \({{\rm{Tr}}}_{{\rm{env}}}\) denotes the trace over the degrees of freedom of the environments. These reduced density matrix elements can be rewritten as

$$\begin{array}{rcl}{\rho }_{ij}^{(\mathrm{sys})}(t) & = & \sum _{n\mathrm{=1}}^{3}\langle i|n\rangle \langle n|{{\rm{Tr}}}_{{\rm{env}}}\rho (t)|j\rangle \\ & = & \sum _{n\mathrm{=1}}^{3}\langle n|{{\rm{Tr}}}_{{\rm{env}}}\rho (t)|j\rangle \langle i|n\rangle \\ & = & {\rm{Tr}}[\rho (t)|j\rangle \langle i|],\end{array}$$
(8)

where Tr denotes the trace over the degrees of freedom of the total system. Thus, \({\rho }_{ij}^{(\mathrm{sys})}(t)\) is just the expectation value of the system operator \(|j\rangle \langle i|\) and it can also be written, in the Heisenberg picture, as

$$\begin{array}{rcl}{\rho }_{ij}^{(\mathrm{sys})}(t) & = & \langle {\rm{\Psi }}(t)|j\rangle \langle i|{\rm{\Psi }}(t)\rangle \\ & = & \langle {{\rm{\Psi }}}_{0}|j\rangle \langle i|(t)|{{\rm{\Psi }}}_{0}\rangle ,\end{array}$$
(9)

with

$$|j\rangle \langle i|(t)\equiv {U}^{\dagger }(t)|j\rangle \langle i|U(t),$$
(10)

where \(|j\rangle \langle i|(t)\) is a system operator represented in the Heisenberg picture and \(|j\rangle \langle i|\) is this operator at the initial time t = 0, while U(t) is the evolution operator,

$$U(t)=T{e}^{-i{\int }_{0}^{t}{H}_{{\rm{tot}}}(s)ds},$$
(11)

with T being the time-ordering operator. Here we study the case with the initial state of the total system given \(|{{\rm{\Psi }}}_{0}\rangle \equiv |D\mathrm{(0)}\rangle \otimes \mathrm{|0}\rangle \).

To conveniently see the dynamical behavior of the system from non-Markovian to Markovian, we choose the correlation functions \({\alpha }_{i}(t,s)\equiv {\sum }_{k}{|{g}_{i,k}|}^{2}{e}^{-i{\omega }_{k}(t-s)}\) of the two baths as the typical Ornstein-Uhlenbeck correlation functions, \({\alpha }_{i}(t,s)=\frac{{{\rm{\Gamma }}}_{i}{\gamma }_{i}}{2}{e}^{-{\gamma }_{i}|t-s|}\), where i = a, b. The non-Markovian to Markovian transition can be demonstrated by tuning the parameters γ i , i.e., the inverse of the correlation times of the two baths. The coupling strength between the system and the ith bath is given by Γ i which corresponds to the decay rate under Markovian approximation23. Using the Heisenberg equation, we can derive the following non-Markovian quantum Bloch equation for the expectation values of the system operators (see method):

$$\frac{\partial }{\partial t}{\mathscr{A}}(t)= {\mathcal H} (t){\mathscr{A}}(t)+{ {\mathcal L} }_{a}{{\mathscr{A}}}^{\mathrm{(1,0)}}(t)+{ {\mathcal L} }_{b}{{\mathscr{A}}}^{\mathrm{(0,1)}}(t),$$
(12)

where \({\mathscr{A}}(t)\equiv ({{\mathscr{A}}}_{ij}(t))\), with \(i,j=1,2,3\), and \({{\mathscr{A}}}_{ij}(t)=\langle {{\rm{\Psi }}}_{0}|i\rangle \langle j|(t)|{{\rm{\Psi }}}_{0}\rangle \equiv {\rho }_{ij}^{(\mathrm{sys})}(t)\). In Eq. (12), \( {\mathcal H} (t)\), \({ {\mathcal L} }_{a}\) and \({ {\mathcal L} }_{b}\) are given by

$${\mathscr{H}}(t)=(\begin{array}{ccccccccc}0 & 0 & 0 & i{{\rm{\Omega }}}_{1}(t) & 0 & 0 & -i{{\rm{\Omega }}}_{1}^{\ast }(t) & 0 & 0\\ 0 & 0 & 0 & 0 & i{{\rm{\Omega }}}_{2}(t) & 0 & 0 & -i{{\rm{\Omega }}}_{2}^{\ast }(t) & 0\\ 0 & 0 & 0 & -i{{\rm{\Omega }}}_{1}(t) & -i{{\rm{\Omega }}}_{2}(t) & 0 & i{{\rm{\Omega }}}_{1}^{\ast }(t) & i{{\rm{\Omega }}}_{2}^{\ast }(t) & 0\\ i{{\rm{\Omega }}}_{1}^{\ast }(t) & 0 & -i{{\rm{\Omega }}}_{1}^{\ast }(t) & i{{\rm{\Delta }}}_{31} & 0 & i{{\rm{\Omega }}}_{2}^{\ast }(t) & 0 & 0 & 0\\ 0 & i{{\rm{\Omega }}}_{2}^{\ast }(t) & -i{{\rm{\Omega }}}_{2}^{\ast }(t) & 0 & i{{\rm{\Delta }}}_{32} & 0 & 0 & 0 & i{{\rm{\Omega }}}_{1}^{\ast }(t)\\ 0 & 0 & 0 & i{{\rm{\Omega }}}_{2}(t) & 0 & i{{\rm{\Delta }}}_{21} & 0 & -i{{\rm{\Omega }}}_{1}^{\ast }(t) & 0\\ -i{{\rm{\Omega }}}_{1}(t) & 0 & i{{\rm{\Omega }}}_{1}(t) & 0 & 0 & 0 & -i{{\rm{\Delta }}}_{31} & 0 & -i{{\rm{\Omega }}}_{2}(t)\\ 0 & -i{{\rm{\Omega }}}_{2}(t) & i{{\rm{\Omega }}}_{2}(t) & 0 & 0 & -i{{\rm{\Omega }}}_{1}(t) & 0 & -i{{\rm{\Delta }}}_{32} & 0\\ 0 & 0 & 0 & 0 & i{{\rm{\Omega }}}_{1}(t) & 0 & -i{{\rm{\Omega }}}_{2}^{\ast }(t) & 0 & -i{{\rm{\Delta }}}_{21}\end{array}),$$
(13)
$$\begin{array}{c}{{\mathscr{L}}}_{a}=(\begin{array}{ccccccccc}0 & 0 & 0 & i & 0 & 0 & -i & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & -i & 0 & 0 & i & 0 & 0\\ i & 0 & -i & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & i\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & -i & 0\\ -i & 0 & i & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & -i & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & i & 0 & 0 & 0 & 0\end{array}),\\ {{\mathscr{L}}}_{b}=(\begin{array}{ccccccccc}0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & i & 0 & 0 & -i & 0\\ 0 & 0 & 0 & 0 & -i & 0 & 0 & i & 0\\ 0 & 0 & 0 & 0 & 0 & i & 0 & 0 & 0\\ 0 & i & -i & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & i & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -i\\ 0 & -i & i & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & -i & 0 & 0\end{array}).\end{array}$$
(14)

where \({{\rm{\Delta }}}_{ij}={\omega }_{i}-{\omega }_{j}\,(i,j=1,2,\mathrm{3)}\) and \({{\rm{\Omega }}}_{i}(t)={{\rm{\Omega }}}_{i}{e}^{-i{{\rm{\Delta }}}_{3i}t}\,(i=1,\mathrm{2)}\). The non-Markovianity of the quantum dynamics of the three-level system is reflected in both \({{\mathscr{A}}}^{\mathrm{(1,0)}}(t)\) and \({{\mathscr{A}}}^{\mathrm{(0,1)}}(t)\), which are solved via the hierarchical equation

$$\begin{array}{rcl}\frac{\partial }{\partial t}{{\mathscr{A}}}^{(m,n)}(t) & = & -(m{\gamma }_{a}+n{\gamma }_{b}){{\mathscr{A}}}^{(m,n)}(t)+ {\mathcal H} (t){{\mathscr{A}}}^{(m,n)}(t)\\ & & +\,m\frac{{{\rm{\Gamma }}}_{a}{\gamma }_{a}}{2}{ {\mathcal L} }_{a}{{\mathscr{A}}}^{(m-\mathrm{1,}n)}(t)+{ {\mathcal L} }_{a}{{\mathscr{A}}}^{(m+\mathrm{1,}n)}(t)\\ & & +\,n\frac{{{\rm{\Gamma }}}_{b}{\gamma }_{b}}{2}{ {\mathcal L} }_{b}{{\mathscr{A}}}^{(m,n-\mathrm{1)}}(t)+{ {\mathcal L} }_{b}{{\mathscr{A}}}^{(m,n+\mathrm{1)}}(t),\end{array}$$
(15)

where \({{\mathscr{A}}}^{(m,n)}(t)=0\) if m or n < 0, \({{\mathscr{A}}}^{\mathrm{(0,0)}}(t)\equiv {\mathscr{A}}(t)\), and the initial condition is \({{\mathscr{A}}}^{(m,n)}\mathrm{(0)}=0\) for m or \(n\ne 0\).

With the reduced density operator \({\rho }_{{\rm{sys}}}(t)\) obtained by choosing the initial state of the total system as \(|{{\rm{\Psi }}}_{0}\rangle \equiv |D\mathrm{(0)}\rangle \otimes \mathrm{|0}\rangle \), the fidelity of the dark state of the three-level system can be written as

$$ {\mathcal F} (t)\equiv \sqrt{\langle D(t)|{\rho }_{{\rm{sys}}}(t)|D(t)\rangle }\mathrm{.}$$
(16)

This quantity can be used to characterize the leakage of the dark state to other levels.

Breakdown of the dark state

As shown in Eq. (5), for the interaction Hamiltonian \({H}_{{\rm{RWA}}}\) with only the rotating terms, the dark state persists when the three-level system couples to two zero-temperature baths. Thus, \( {\mathcal F} (t\mathrm{)=1}\) in this case. Below we demonstrate the dynamical evolution of the fidelity \( {\mathcal F} (t)\) of the dark state when the counter-rotating terms are included in the interaction Hamiltonian. For simplicity, we choose the same parameters for the two baths, i.e., \({\alpha }_{a}(t,s)={\alpha }_{b}(t,s)=\frac{{\rm{\Gamma }}\gamma }{2}{e}^{-\gamma |t-s|}\). We also choose the level of state |3〉 to be the zero point of the energy. The energy difference between |1〉 and |3〉 is taken as \({\omega }_{3}-{\omega }_{1}=\omega =1\). Other parameters with the frequency unit are expressed as the ratio to ω. To numerically calculate the quantum dynamics of the system, we need to truncate Eq. (15) at a given hierarchical order \({\mathscr{N}}\) so that \({{\mathscr{A}}}^{(m,n)}(t=\mathrm{0)}\) for \(m+n > {\mathscr{N}}\). As shown in Fig. 2(a), the numerical results at \({\mathscr{N}}=4\) (blue curve) are very close to the result at \({\mathscr{N}}=10\) (read dots), indicating that the results converge rapidly with the hierarchical order \({\mathscr{N}}\). In this case, we choose Γ = 1. Note that the coupling strength between the system and the bath can be characterized by Γ which corresponds to the decay rate of the system under Markovian approximation23. To distinguish different regimes of interactions between the system and the bath, we define Γ = 1 to be the deep-strong coupling regime, because the corresponding Markovian decay rate is now comparable to the system frequencies. Consequently, Γ = 0.01 and 0.1 are defined as the strong and ultrastrong coupling regimes, respectively. When Γ < 1, more accurate results are obtained at the same hierarchical order due to the even faster convergence of the numerical results in the weak coupling regime between the system and the baths. Thus, the truncation of Eq. (15) at \({\mathscr{N}}=4\) can already give reliable results for the model that we study here. However, we choose \({\mathscr{N}}=10\) in our following calculations to have the obtained results more accurate. Also, Fig. 2(b) and (c) show the difference between the results of \({\mathscr{N}}=10\) and \({\mathscr{N}}=20\). It is found that the difference is below 10−5, further indicating that the truncation is reliable.

Figure 2
figure 2

Fidelity evolution of the dark state. (a) Fidelity of the dark state calculated at different order \({\mathscr{N}}\) of the hierarchical equation in Eq. (15). The parameters of the two baths are Γ = 1 and γ = 0.6. (b) and (c) Fidelity difference \({\rm{\Delta }}F\equiv {F}_{10}-{F}_{20}\) at different values of the bath parameters, where \({F}_{10}\) (\({F}_{20}\)) is the fidelity of the dark state calculated using the hierarchical equation up to the order of \({\mathscr{N}}=10\) (20). (d) Fidelity of the dark state calculated at different values of Γ, where γ = 0.6 and \({\mathscr{N}}=10\). (e) Fidelity of the dark state calculated at different values of γ, where Γ = 0.1 and \({\mathscr{N}}=10\). In this figure and the following one, the parameters of the three-level system are \({\omega }_{1}=-1\), \({\omega }_{2}=-1.2\), \({\omega }_{3}=0\), \({{\rm{\Omega }}}_{1}=0.5\), and \({{\rm{\Omega }}}_{2}=0.2\).

Full size image

In Fig. 2(d), we show the fidelity evolution of the dark state by varying the coupling strength Γ between the system and the two baths. When Γ = 1, the fidelity of the dark state decays fast and then quickly reaches the stationary oscillations in this strong system-bath coupling regime. These stationary oscillations correspond to a dynamic equilibrium of the dark state under the combined actions of the drive fields and the baths. When decreasing the coupling strength Γ, the fidelity of the dark state decays slowly and takes a long time to reach the stationary oscillations. However, with a given correlation time (inverse of γ) of the two baths, the dark state has already exhibited appreciable leakage at Γ = 0.01 [see the blue curve in Fig. 2(d)]. In Fig. 2(e), we also show the effect of the environment correlation time (i.e., γ) on the fidelity evolution of the dark state. For the non-Markovian environment, which has a longer correlation time (i.e., smaller γ), the fidelity of the dark state decays slowly with the evolution t, in comparison with the Markovian environment with a shorter correlation time (i.e., larger γ). This reveals that the dark state leaks to other levels more slowly when it is coupled to a non-Markovian environment. The physical intuition for why the Markovian environment leads to a faster decay in fidelity may be related to the memory effect of the environment. Actually, owing to the memory effect of the environment, some information that is leaked into the bath can come back to the system in the non-Markovian case. Moreover, similar to Fig. 2(d), the fidelity of the dark state exhibits stationary oscillations at longer times even in the Markovian limit (large γ) of the baths. These oscillations are also due to the persistent drive applied to the system.

From our results above, we can conclude that when the coupling strength between the system and the bath becomes strong, the dark state is unstable under the influence of the counter-rotating terms even when the environments are at zero temperature. Thus, owing to the counter-rotating terms in the interaction Hamiltonian, the environment-induced effect on the dynamical evolution of the dark state cannot be eliminated by simply lowering the temperature of the environments.

Leakage reduction of the dark state

To reduce the leakage of the dark state, we introduce a leakage elimination operator24,25,26,27,28,29,30. When this leakage elimination operator is added, the total Hamiltonian of the open system becomes

$${H}_{{\rm{tot}}}={H}_{{\rm{sys}}}+{H}_{{\rm{bath}}}+{H}_{{\rm{int}}}+R(t),$$
(17)

where

$$R(t)=-c(t\mathrm{)(|1}\rangle \langle \mathrm{1|}+\mathrm{|2}\rangle \langle \mathrm{2|)}$$
(18)

is the leakage elimination operator that suppresses leakage from the system to the environment. In numerical calculations, we use the rectangular control pulse \(c(t)\) as an example, which has a period of τ [see Fig. 1(b)]. Within the time intervals \(l\tau \le t < l\tau +{\rm{\Delta }}\), where l is a positive integer to identify the pulse is in the lth period, the control pulse is switched on and has an intensity of h. For other times, the control pulse is turned off.

From Fig. 3(a), it can be seen that when the control pulse is applied, the fidelity of the dark state (gray and brown curves) decreases much slowly with the evolution time, in sharp contrast with the fidelity without the control pulse (red curve). Thus, the leakage elimination operator works quite effectively in reducing the leakage of the dark state. Also, the brown curve has a higher fidelity than the gray one, implying that a higher pulse intensity yields better control effect for suppressing the state leakage. Figure 3(b) presents the effect of the period of the control pulse on the leakage elimination. It is clear that the fidelity of dark state increases when decreasing the period of the control pulse. Therefore, both higher intensity and frequency of the control pulse can strengthen the effect of the leakage elimination on the dark state. In Fig. 3(c) and (d), we further show the effect of the inverse of the correlation time γ. Compared with Fig. 2(c), it can be seen that the leakage of the dark state is much eliminated for very non-Markovian baths with small values of γ. This is the distinct advantage of the non-Markovian baths.

Figure 3
figure 3

Fidelity evolution of the dark state in the presence of the control pulses. (a) Fidelity of the dark state when varying the strength of the pulse, where the period of the pulses is \(\tau =0.2\) and the duration of each pulse is \({\rm{\Delta }}=0.1\). The parameters of the two baths are Γ = 0.1 and γ = 0.6. (b) Fidelity of the dark state when varying the period of the pulses, where the strength of the pulse is \(h=\mathrm{0.6/}{\rm{\Delta }}\) and the duration is chosen to be \({\rm{\Delta }}\mathrm{=0.5}\tau \). The parameters of the two baths are also Γ = 0.1 and γ = 0.6. (c) Fidelity of the dark state when varying the inverse of the bath correlation time γ, where the strength of the pulse is h = 6, the period of the pulses is τ = 0.2, and the duration of each pulse is \({\rm{\Delta }}\mathrm{=0.1}\). The coupling strength of each bath is Γ = 0.1. (d) Fidelity of the dark state when varying the inverse of the bath correlation time γ, where the coupling strength of each bath to the system is Γ = 1. Other parameters are the same as in (c).

Full size image

Finally, we discuss more about the leakage elimination operator. Instead of using the eigenstates |1〉, |2〉 and |3〉, we can also use the dark state \(|D(t)\rangle \), bright state \(|B(t)\rangle \) and eigenstate |3〉 as the basis states of the Hilbert space of the three-level system, where the dark state is given in Eq. (2) and the bright state is

$$|B(t)\rangle =\frac{{{\rm{\Omega }}}_{1}}{{\rm{\Omega }}}{e}^{-i{\omega }_{1}t}|1\rangle +\frac{{{\rm{\Omega }}}_{2}}{{\rm{\Omega }}}{e}^{-i{\omega }_{2}t}|2\rangle \mathrm{.}$$
(19)

With this set of basis states, we can rewrite the interaction Hamiltonian in Eq. (3) as

$$\begin{array}{rcl}{H}_{{\rm{int}}} & = & \frac{{{\rm{\Omega }}}_{2}^{\ast }}{{\rm{\Omega }}}{e}^{i{\omega }_{1}t}|D(t)\rangle \langle 3|\sum _{k}({g}_{a,k}{a}_{k}+{g}_{a,k}^{\ast }{a}_{k}^{\dagger })\\ & & -\frac{{{\rm{\Omega }}}_{1}^{\ast }}{{\rm{\Omega }}}{e}^{i{\omega }_{2}t}|D(t)\rangle \langle 3|\sum _{k}({g}_{b,k}{b}_{k}+{g}_{a,k}^{\ast }{b}_{k}^{\dagger })\\ & & +\,\frac{{{\rm{\Omega }}}_{1}^{\ast }}{{\rm{\Omega }}}{e}^{i{\omega }_{1}t}|B(t)\rangle \langle 3|\sum _{k}({g}_{a,k}{a}_{k}+{g}_{a,k}^{\ast }{a}_{k}^{\dagger })\\ & & +\,\frac{{{\rm{\Omega }}}_{2}^{\ast }}{{\rm{\Omega }}}{e}^{i{\omega }_{2}t}|B(t)\rangle \langle 3|\sum _{k}({g}_{b,k}{b}_{k}+{g}_{a,k}^{\ast }{b}_{k}^{\dagger })+{\rm{H}}{\rm{.c}}{\rm{.}}\end{array}$$
(20)

From Eq. (20), it follows that there is no direct leakage from the dark state to the bright state, so the leakage of the dark state occurs only via the excited state |3〉. If the state |3〉 is kept unoccupied, the leakage of the dark state can be prevented. That is why we can use the leakage elimination operator in Eq. (18) to restore the quantum coherence of the dark state. As shown in refs 31–34, counterdiabatic driving usually introduces an additional term to cancel the effect of the nonadiabatic transition, while the leakage elimination operator can effectively enlarge the energy difference of the system. According to the adiabatic condition [see, e.g., Eq. (5) in ref. 32], this leakage elimination operator can have the similar effect as the counterdiabatic driving.

Discussion

We have studied the non-Markovian quantum dynamics of the dark state in a \({\rm{\Lambda }}\)-type three-level system coupled to two bosonic baths and revealed the effect of the counter-rotating terms in the system-bath coupling on the dark state. Due to these counter-rotating terms, the dark state leaks to other states even at zero temperature, in sharp contrast to the dark state within the RWA. Thus, whether the dark state is really dark or not depends on the validity of the RWA in the considered system. Actually, our numerical results have shown that the counter-rotating terms cannot be neglected for a strong system-bath interaction, because these terms yield appreciable leakage of the dark state in this strong coupling regime. To restore the quantum coherence of the dark state, we propose to apply a leakage elimination operator to the system. Our numerical results indicate that the leakage of the dark state can indeed be greatly suppressed with the help of this leakage elimination operator.

Our study reveals a possible mechanism for the dark state to leak and a way to fight against it. This may improve the precision of the experiments related to the dark state. While the studies of quantum dynamics beyond the RWA mainly focus on the two-level systems, our work provides insights into the dynamics of the three-level system beyond the RWA.

Methods

From the total Hamiltonian \({H}_{{\rm{tot}}}\) given in Sec. II, it follows that the Heisenberg equations for the field operators \({a}_{k}\) and \({b}_{k}\) are given by

$$\begin{array}{rcl}\frac{d}{dt}{a}_{k}(t) & = & -i{\omega }_{k}{a}_{k}(t)-i{g}_{a,k}^{\ast }(|1\rangle \langle 3|(t)+|3\rangle \langle 1|(t)),\\ \frac{d}{dt}{b}_{k}(t) & = & -i{\omega }_{k}{b}_{k}(t)-i{g}_{b,k}^{\ast }(|2\rangle \langle 3|(t)+|3\rangle \langle 2|(t)),\end{array}$$
(21)

were \({a}_{k}(t)\equiv {U}^{\dagger }(t){a}_{k}U(t)\) and \({b}_{k}(t)\equiv {U}^{\dagger }(t){b}_{k}U(t)\), with \(U(t)\) given in Eq. (11). Equation (21) can be formally solved as

$$\begin{array}{rcl}{a}_{k}(t) & = & {e}^{-i{\omega }_{k}t}{a}_{k}-i{g}_{a,k}^{\ast }{\int }_{0}^{t}ds{e}^{-i{\omega }_{k}(t-s)}[|1\rangle \langle 3|(s)+|3\rangle \langle 1|(s)],\\ {b}_{k}(t) & = & {e}^{-i{\omega }_{k}t}{b}_{k}-i{g}_{b,k}^{\ast }{\int }_{0}^{t}ds{e}^{-i{\omega }_{k}(t-s)}[|2\rangle \langle 3|(s)+|3\rangle \langle 2|(s)]\mathrm{.}\end{array}$$
(22)

Similar to Eq. (21), we can also derive the Heisenberg equations for the system operators \(|i\rangle \langle j|\), with \(i,j=1,2,3\). Then, substituting Eq. (22) into the Heisenberg equations for the system operators, we have

$$\begin{array}{rcl}\frac{d}{dt}|1\rangle \langle 1|(t) & = & i{{\rm{\Omega }}}_{1}{e}^{-i({\omega }_{3}-{\omega }_{1})t}\mathrm{|3}\rangle \langle \mathrm{1|(}t)-i{{\rm{\Omega }}}_{1}^{\ast }{e}^{i({\omega }_{3}-{\omega }_{1})t}\mathrm{|1}\rangle \langle \mathrm{3|(}t)+i\mathrm{(|3}\rangle \langle \mathrm{1|(}t)\\ & & -\,\mathrm{|3}\rangle \langle \mathrm{1|(}t))({\xi }_{a}(t)+{\xi }_{a}^{\dagger }(t)),\\ \frac{d}{dt}|2\rangle \langle 2|(t) & = & i{{\rm{\Omega }}}_{2}{e}^{-i({\omega }_{3}-{\omega }_{2})t}\mathrm{|3}\rangle \langle \mathrm{2|(}t)-i{{\rm{\Omega }}}_{2}^{\ast }{e}^{i({\omega }_{3}-{\omega }_{2})t}\mathrm{|2}\rangle \langle \mathrm{3|(}t)+i\mathrm{(|3}\rangle \langle \mathrm{2|(}t)\\ & & -\,\mathrm{|2}\rangle \langle \mathrm{3|(}t))({\xi }_{b}(t)+{\xi }_{b}^{\dagger }(t)),\\ \frac{d}{dt}|3\rangle \langle 3|(t) & = & i{{\rm{\Omega }}}_{1}^{\ast }{e}^{i({\omega }_{3}-{\omega }_{1})t}\mathrm{|1}\rangle \langle \mathrm{3|(}t)-i{{\rm{\Omega }}}_{1}{e}^{-i({\omega }_{3}-{\omega }_{1})t}\mathrm{|3}\rangle \langle \mathrm{1|(}t)+i{{\rm{\Omega }}}_{2}^{\ast }{e}^{i({\omega }_{3}-{\omega }_{2})t}\mathrm{|2}\rangle \langle \mathrm{3|(}t)\\ & & -\,i{{\rm{\Omega }}}_{2}{e}^{-i({\omega }_{3}-{\omega }_{2})t}\mathrm{|3}\rangle \langle \mathrm{2|(}t)+i\mathrm{(|1}\rangle \langle \mathrm{3|(}t)-\mathrm{|3}\rangle \langle \mathrm{1|(}t))({\xi }_{a}(t)+{\xi }_{a}^{\dagger }(t))\\ & & +\,i\mathrm{(|2}\rangle \langle \mathrm{3|(}t)-\mathrm{|3}\rangle \langle \mathrm{2|(}t))({\xi }_{b}(t)+{\xi }_{b}^{\dagger }(t)),\\ \frac{d}{dt}|3\rangle \langle 1|(t) & = & i({\omega }_{3}-{\omega }_{1}\mathrm{)|3}\rangle \langle \mathrm{1|(}t)+i{{\rm{\Omega }}}_{1}^{\ast }{e}^{i({\omega }_{3}-{\omega }_{1})t}\mathrm{(|1}\rangle \langle \mathrm{1|(}t)-\mathrm{|3}\rangle \langle \mathrm{3|(}t))\\ & & +\,i{{\rm{\Omega }}}_{2}^{\ast }{e}^{i({\omega }_{3}-{\omega }_{2})t}\mathrm{|2}\rangle \langle \mathrm{1|(}t)+i\mathrm{(|1}\rangle \langle \mathrm{1|(}t)-\mathrm{|3}\rangle \langle \mathrm{3|(}t))({\xi }_{a}(t)\\ & & +\,{\xi }_{a}^{\dagger }(t))+i\mathrm{|2}\rangle \langle \mathrm{1|(}t)({\xi }_{b}(t)+{\xi }_{b}^{\dagger }(t)),\\ \frac{d}{dt}|3\rangle \langle 2|(t) & = & i({\omega }_{3}-{\omega }_{2}\mathrm{)|3}\rangle \langle \mathrm{2|(}t)+i{{\rm{\Omega }}}_{2}^{\ast }{e}^{i({\omega }_{3}-{\omega }_{2})t}\mathrm{(|2}\rangle \langle \mathrm{2|(}t)-\mathrm{|3}\rangle \langle \mathrm{3|(}t))\\ & & +\,i{{\rm{\Omega }}}_{1}^{\ast }{e}^{i({\omega }_{3}-{\omega }_{1})t}\mathrm{|1}\rangle \langle \mathrm{2|(}t)+i\mathrm{(|2}\rangle \langle \mathrm{2|(}t)-\mathrm{|3}\rangle \langle \mathrm{3|(}t))({\xi }_{b}(t)\\ & & +\,{\xi }_{b}^{\dagger }(t))+i\mathrm{|1}\rangle \langle \mathrm{2|(}t)({\xi }_{a}(t)+{\xi }_{a}^{\dagger }(t)),\\ \frac{d}{dt}|2\rangle \langle 1|(t) & = & i({\omega }_{2}-{\omega }_{1}\mathrm{)|2}\rangle \langle \mathrm{1|(}t)+i{{\rm{\Omega }}}_{2}{e}^{-i({\omega }_{3}-{\omega }_{2})}\mathrm{|3}\rangle \langle \mathrm{1|(}t)-i{{\rm{\Omega }}}_{1}^{\ast }{e}^{i({\omega }_{3}-{\omega }_{1})t}\mathrm{|2}\rangle \langle \mathrm{3|(}t)\\ & & +\,i\mathrm{|3}\rangle \langle \mathrm{1|(}t)({\xi }_{b}(t)+{\xi }_{b}^{\dagger }(t))-i\mathrm{|2}\rangle \langle \mathrm{3|(}t)({\xi }_{a}(t)+{\xi }_{a}^{\dagger }(t)),\end{array}$$
(23)

with.., \(\mathrm{|2}\rangle \langle \mathrm{3|(}t)\), and \(\mathrm{|2}\rangle \langle \mathrm{1|(}t)\) equal to the conjugates of \(\mathrm{|3}\rangle \langle \mathrm{1|(}t)\), \(\mathrm{|3}\rangle \langle \mathrm{2|(}t)\), and \(\mathrm{|1}\rangle \langle \mathrm{2|(}t)\), respectively. Here \(|i\rangle \langle j|(t)\equiv {U}^{\dagger }(t)|i\rangle \langle j|U(t)\), and the noise operators \({\xi }_{a}(t)\) and \({\xi }_{b}(t)\) are defined as \({\xi }_{a}(t)\equiv {\sum }_{k}{g}_{a,k}{e}^{-i{\omega }_{k}t}{a}_{k}\) and \({\xi }_{b}(t)\equiv {\sum }_{k}{g}_{b,k}{e}^{-i{\omega }_{k}t}{b}_{k}\). Also, we have considered the resonant case with \({\omega }_{a}={\omega }_{3}-{\omega }_{1}\) and \({\omega }_{b}={\omega }_{3}-{\omega }_{2}\). Note that Eq. (23) is just the quantum Langvin equation of the system.

Let us introduce the Bargmann coherent states for the two bosonic baths,

$$|z\rangle \equiv \mathop{\otimes }\limits_{k}|{z}_{k}\rangle ={e}^{\sum _{k}{z}_{ak}{a}_{k}^{\dagger }+\sum _{k}{z}_{bk}{b}_{k}^{\dagger }}|0\rangle ,$$
(24)

which satisfies

$${a}_{k}|z\rangle ={z}_{ak}|z\rangle ,\,{a}_{k}^{\dagger }|z\rangle =\frac{\partial }{\partial {z}_{ak}}|z\rangle ,\,{b}_{k}|z\rangle ={z}_{bk}|z\rangle ,\,{b}_{k}^{\dagger }|z\rangle =\frac{\partial }{\partial {z}_{bk}}|z\rangle \mathrm{.}$$
(25)

As in ref. 35 for the spin-boson model, we project Eq. (23) onto the Bargmann coherent states and take the expectation value for each operator. Then, we convert the quantum Langevin equation in Eq. (23) to

$$\frac{\partial }{\partial t}{\mathscr{A}}(t,z)= {\mathcal H} (t){\mathscr{A}}(t,z)+[{z}_{at}+{\int }_{0}^{t}ds{\alpha }_{a}(t,s)\frac{\delta }{\delta {z}_{as}}]{ {\mathcal L} }_{a}{\mathscr{A}}(t,z)+[{z}_{bt}+{\int }_{0}^{t}ds{\alpha }_{b}(t,s)\frac{\delta }{\delta {z}_{bs}}]{ {\mathcal L} }_{b}{\mathscr{A}}(t,z),$$
(26)

where \(\,{\mathscr{A}}(t,z)\equiv {({{\mathscr{A}}}_{11}(t,z),{{\mathscr{A}}}_{22}(t,z),{{\mathscr{A}}}_{33}(t,z),{{\mathscr{A}}}_{31}(t,z),{{\mathscr{A}}}_{32}(t,z),{{\mathscr{A}}}_{21}(t,z),{{\mathscr{A}}}_{13}(t,z),{{\mathscr{A}}}_{23}(t,z),{{\mathscr{A}}}_{12}(t,z))}^{T}\), with T denoting the transpose of a matrix and \({{\mathscr{A}}}_{ij}(t,z)\equiv \langle {{\rm{\Psi }}}_{0}|i\rangle \langle j|(t)|z\rangle \langle z|{{\rm{\Psi }}}_{0}\rangle \). The matrix \( {\mathcal H} (t)\) is given in Eq. (13), and the matrices \({ {\mathcal L} }_{a}\) and \({ {\mathcal L} }_{b}\) are given in Eq. (14). Also, in deriving Eq. (26), we have used the functional chain rule

$$\frac{\partial }{\partial {z}_{k}}=\int ds\frac{\partial {z}_{s}}{\partial {z}_{k}}\frac{\delta }{\delta {z}_{s}},$$
(27)

and defined the noise functions,

$${z}_{at}=\sum _{k}{g}_{a,k}{e}^{-i{\omega }_{k}t}{z}_{ak},\,{z}_{bt}=\sum _{k}{g}_{b,k}{e}^{-i{\omega }_{k}t}{z}_{bk}\mathrm{.}$$
(28)

Below we solve Eq. (26) using a hierarchical-equation appraoch36. We extend the limits of integrals in Eq. (26) to infinity and define \({D}_{a}\equiv {\int }_{-\infty }^{+\infty }{\alpha }_{a}(t,s)\frac{\delta }{\delta {z}_{as}}\) and \({D}_{b}\equiv {\int }_{-\infty }^{+\infty }{\alpha }_{b}(t,s)\frac{\delta }{\delta {z}_{bs}}\). According to Eq. (26), \({\mathscr{A}}(t,z)\) only contains the noise for \(0\le s < t\). Then, for other values of s, we have \(\frac{\delta }{\delta {z}_{is}}{\mathscr{A}}(t,z)=0\), where \(i=a,b\). Thus, the extension of the integral limit is exact. Now Eq. (26) can be expressed as

$$\frac{\partial }{\partial t}{\mathscr{A}}(t,z)= {\mathcal H} (t){\mathscr{A}}(t,z)+({z}_{at}+{D}_{a}){ {\mathcal L} }_{a}{\mathscr{A}}(t,z)+({z}_{bt}+{D}_{b}){ {\mathcal L} }_{b}{\mathscr{A}}(t,z\mathrm{).}$$
(29)

We further define \({{\mathscr{A}}}^{(m,n)}(t,z)\equiv {D}_{a}^{m}{D}_{b}^{n}{\mathscr{A}}(t,z)\) and consider the Ornstein-Uhlenbeck correlation functions \({\alpha }_{i}(t,s)=\frac{{{\rm{\Gamma }}}_{i}{\gamma }_{i}}{2}{e}^{-{\gamma }_{i}|t-s|}\), with \(i=a,b\). Then, it follows that

$$\begin{array}{rcl}\frac{\partial }{\partial t}{{\mathscr{A}}}^{(m,n)}(t,z) & = & \frac{\partial }{\partial t}{[{\int }_{-\infty }^{\infty }ds{\alpha }_{a}(t,s)\frac{\delta }{\delta {z}_{as}}]}^{m}{[{\int }_{-\infty }^{\infty }ds{\alpha }_{b}(t,s)\frac{\delta }{\delta {z}_{bs}}]}^{n}{\mathscr{A}}(t,z)\\ & = & -(m{\gamma }_{a}+n{\gamma }_{b}){{\mathscr{A}}}^{(m,n)}(t,z)\\ & & +\,{[{\int }_{-\infty }^{\infty }ds{\alpha }_{a}(t,s)\frac{\delta }{\delta {z}_{as}}]}^{m}{[{\int }_{-\infty }^{\infty }ds{\alpha }_{b}(t,s)\frac{\delta }{\delta {z}_{bs}}]}^{n}\frac{\partial }{\partial t}{\mathscr{A}}(t,z\mathrm{).}\\ & = & -(m{\gamma }_{a}+n{\gamma }_{b}){{\mathscr{A}}}^{(m,n)}(t,z)+{D}_{a}^{m}{D}_{b}^{n}\frac{\partial }{\partial t}{\mathscr{A}}(t,z\mathrm{).}\end{array}$$
(30)

Using Eq. (29), we obtain

$$\begin{array}{c}{D}_{a}^{m}{D}_{b}^{n}\frac{\partial }{\partial t}{\mathscr{A}}(t,z)={D}_{a}^{m}{D}_{b}^{n}[ {\mathcal H} (t){\mathscr{A}}(t,z)+{z}_{at}{ {\mathcal L} }_{a}{\mathscr{A}}(t,z)+{ {\mathcal L} }_{a}{D}_{at}{\mathscr{A}}(t,z)\\ \quad \quad \quad \quad \quad \quad \quad \quad +\,{z}_{bt}{ {\mathcal L} }_{b}{\mathscr{A}}(t,z)+{ {\mathcal L} }_{b}{D}_{bt}{\mathscr{A}}(t,z)]\\ \quad \quad \quad \quad \quad \quad \quad ={D}_{at}^{m}{D}_{bt}^{n} {\mathcal H} (t){\mathscr{A}}(t,z)+{D}_{at}^{m}{D}_{bt}^{n}{z}_{at}{ {\mathcal L} }_{a}{\mathscr{A}}(t,z)\\ \quad \quad \quad \quad \quad \quad \quad \quad +\,{D}_{at}^{m}{D}_{bt}^{n}{ {\mathcal L} }_{a}{D}_{at}{\mathscr{A}}(t,z)\\ \quad \quad \quad \quad \quad \quad \quad \quad +\,{D}_{at}^{m}{D}_{bt}^{n}{z}_{bt}{ {\mathcal L} }_{b}{\mathscr{A}}(t,z)+{D}_{at}^{m}{D}_{bt}^{n}{ {\mathcal L} }_{b}{D}_{bt}{\mathscr{A}}(t,z)\\ \quad \quad \quad \quad \quad \quad \quad = {\mathcal H} (t){{\mathscr{A}}}^{(m,n)}(t,z)+m\frac{{{\rm{\Gamma }}}_{a}{\gamma }_{a}}{2}{ {\mathcal L} }_{a}{D}_{at}^{m-1}{D}_{bt}^{n}{\mathscr{A}}(t,z)\\ \quad \quad \quad \quad \quad \quad \quad \quad +\,{z}_{at}{ {\mathcal L} }_{a}{{\mathscr{A}}}^{(m,n)}(t,z)+{ {\mathcal L} }_{a}{D}_{at}^{m+1}{D}_{bt}^{n}{\mathscr{A}}(t,z)\\ \quad \quad \quad \quad \quad \quad \quad \quad +\,n\frac{{{\rm{\Gamma }}}_{b}{\gamma }_{b}}{2}{ {\mathcal L} }_{b}{D}_{at}^{m}{D}_{bt}^{n-1}{\mathscr{A}}(t,z)+{z}_{bt}{ {\mathcal L} }_{b}{D}_{at}^{m}{D}_{bt}^{n}{\mathscr{A}}(t,z)\\ \quad \quad \quad \quad \quad \quad \quad \quad +\,{ {\mathcal L} }_{b}{D}_{at}^{m}{D}_{bt}^{n+1}{\mathscr{A}}(t,z)\\ \quad \quad \quad \quad \quad \quad \quad = {\mathcal H} (t){{\mathscr{A}}}^{(m,n)}(t,z)+m\frac{{{\rm{\Gamma }}}_{a}{\gamma }_{a}}{2}{ {\mathcal L} }_{a}{{\mathscr{A}}}^{(m-\mathrm{1,}n)}(t,z)\\ \quad \quad \quad \quad \quad \quad \quad \quad +\,{z}_{at}{ {\mathcal L} }_{a}{{\mathscr{A}}}^{(m,n)}(t,z)+{ {\mathcal L} }_{a}{{\mathscr{A}}}^{(m+\mathrm{1,}n)}(t,z)+n\frac{{{\rm{\Gamma }}}_{b}{\gamma }_{b}}{2}{ {\mathcal L} }_{b}{{\mathscr{A}}}^{(m,n-\mathrm{1)}}(t,z)\\ \quad \quad \quad \quad \quad \quad \quad \quad +\,{z}_{bt}{ {\mathcal L} }_{b}{{\mathscr{A}}}^{(m,n)}(t,z)+{ {\mathcal L} }_{b}{{\mathscr{A}}}^{(m,n+\mathrm{1)}}(t,z\mathrm{).}\end{array}$$
(31)

Substitution of Eq. (31) into Eq. (30) gives

$$\begin{array}{rcl}\frac{\partial }{\partial t}{{\mathscr{A}}}^{(m,n)}(t,z) & = & -(m{\gamma }_{a}+n{\gamma }_{b}){{\mathscr{A}}}^{(m,n)}(t,z)+ {\mathcal H} (t){{\mathscr{A}}}^{(m,n)}(t,z)\\ & & +\,m\frac{{{\rm{\Gamma }}}_{a}{\gamma }_{a}}{2}{ {\mathcal L} }_{a}{{\mathscr{A}}}^{(m-\mathrm{1,}n)}(t,z)+{z}_{at}{ {\mathcal L} }_{a}{{\mathscr{A}}}^{(m,n)}(t,z)\\ & & +\,{ {\mathcal L} }_{a}{{\mathscr{A}}}^{(m+\mathrm{1,}n)}(t,z)+n\frac{{{\rm{\Gamma }}}_{b}{\gamma }_{b}}{2}{ {\mathcal L} }_{b}{{\mathscr{A}}}^{(m,n-\mathrm{1)}}(t,z)\\ & & +\,{z}_{bt}{ {\mathcal L} }_{b}{{\mathscr{A}}}^{(m,n)}(t,z)+{ {\mathcal L} }_{b}{{\mathscr{A}}}^{(m,n+\mathrm{1)}}(t,z\mathrm{).}\end{array}$$
(32)

This is a stochastic differential equation, because it involves the noise functions \({z}_{at}\) and \({z}_{bt}\).

Let \({{\mathscr{A}}}^{(m,n)}(t)\equiv {\mathcal M} \{{{\mathscr{A}}}^{(m,n)}(t,z)\}\), where the statistical average is defined as

$$ {\mathcal M} \{\ldots \}\equiv \prod _{k}\int \frac{{d}^{2}{z}_{ak}{d}^{2}{z}_{bk}}{{\pi }^{2}}{e}^{-|{z}_{ak}{|}^{2}-|{z}_{bk}{|}^{2}}\{\ldots \}\mathrm{.}$$
(33)

It follows from Eq. (32) that

$$\begin{array}{c}\frac{\partial }{\partial t}{{\mathscr{A}}}^{(m,n)}(t)= {\mathcal M} \{\frac{\partial }{\partial t}{{\mathscr{A}}}^{(m,n)}(t,z)\}\\ \quad \quad \quad \quad \,\,\,\,=\, {\mathcal M} \{{z}_{at}{ {\mathcal L} }_{a}{{\mathscr{A}}}^{(m,n)}(t,z)\}+ {\mathcal M} \{{z}_{bt}{ {\mathcal L} }_{b}{{\mathscr{A}}}^{(m,n)}(t,z)\}\\ \quad \quad \quad \quad \,\,\,\,\,\,\,\,-\,(m{\gamma }_{a}+n{\gamma }_{b}){{\mathscr{A}}}^{(m,n)}(t)+ {\mathcal H} (t){{\mathscr{A}}}^{(m,n)}(t)+m\frac{{{\rm{\Gamma }}}_{a}{\gamma }_{a}}{2}{ {\mathcal L} }_{a}{{\mathscr{A}}}^{(m-\mathrm{1,}n)}(t)\\ \quad \quad \quad \quad \,\,\,\,\,\,\,\,+{ {\mathcal L} }_{a}{{\mathscr{A}}}^{(m+\mathrm{1,}n)}(t)+\,n\frac{{{\rm{\Gamma }}}_{b}{\gamma }_{b}}{2}{ {\mathcal L} }_{b}{{\mathscr{A}}}^{(m,n-\mathrm{1)}}(t)+{ {\mathcal L} }_{b}{{\mathscr{A}}}^{(m,n+\mathrm{1)}}(t\mathrm{).}\end{array}$$
(34)

For an arbitrary element \({{\mathscr{A}}}_{ij}^{(m,n)}(t,z)\) of \({{\mathscr{A}}}^{(m,n)}(t,z)\),

$$\begin{array}{rcl} {\mathcal M} \{{z}_{at}{{\mathscr{A}}}_{ij}^{(m,n)}(t,z)\} & = & {\mathcal M} \{{z}_{at}{D}_{at}^{m}{D}_{bt}^{n}\langle {{\rm{\Psi }}}_{0}|i\rangle \langle j|(t)|z\rangle \langle z|{{\rm{\Psi }}}_{0}\rangle \}\\ & = & {\mathcal M} \{\langle {{\rm{\Psi }}}_{0}|i\rangle \langle j|(t)[{\xi }_{a}^{\dagger }(t{)]}^{m}{[{\xi }_{b}^{\dagger }(t)]}^{n}{\xi }_{a}(t)|z\rangle \langle z|{{\rm{\Psi }}}_{0}\rangle \}\\ & = & \langle {{\rm{\Psi }}}_{0}|i\rangle \langle j|(t)[{\xi }_{a}^{\dagger }(t{)]}^{m}{[{\xi }_{b}^{\dagger }(t)]}^{n}{\xi }_{a}(t)|{{\rm{\Psi }}}_{0}\rangle \\ & = & \mathrm{0,}\end{array}$$

because \({\xi }_{a}(t)|{{\rm{\Psi }}}_{0}\rangle ={\sum }_{k}{g}_{a,k}{e}^{-i{\omega }_{k}t}{a}_{k}|D\mathrm{(0)}\rangle \otimes |0\rangle ={\sum }_{k}{g}_{a,k}{e}^{-i{\omega }_{k}t}{a}_{k}|D\mathrm{(0)}\rangle \otimes {|0\rangle }_{a}\otimes {|0\rangle }_{b}=0\). Thus, we have \( {\mathcal M} \{{z}_{at}{ {\mathcal L} }_{a}{{\mathscr{A}}}^{(m,n)}(t,z)\}=0\). Similarly, because \({\xi }_{b}(t)|{{\rm{\Psi }}}_{0}\rangle ={\sum }_{k}{g}_{b,k}{e}^{-i{\omega }_{k}t}{b}_{k}|D\mathrm{(0)}\rangle \otimes |0\rangle =0\), we can also derive that \( {\mathcal M} \{{z}_{bt}{ {\mathcal L} }_{b}{{\mathscr{A}}}^{(m,n)}(t,z)\}=0\). Therefore, Eq. (34) is simplified as

$$\begin{array}{c}\frac{\partial }{\partial t}{{\mathscr{A}}}^{(m,n)}(t)=-(m{\gamma }_{a}+n{\gamma }_{b}){{\mathscr{A}}}^{(m,n)}(t)+ {\mathcal H} (t){{\mathscr{A}}}^{(m,n)}(t)+m\frac{{{\rm{\Gamma }}}_{a}{\gamma }_{a}}{2}{ {\mathcal L} }_{a}{{\mathscr{A}}}^{(m-\mathrm{1,}n)}(t)\\ \quad \quad \quad \quad \quad \,\,\,\,+{ {\mathcal L} }_{a}{{\mathscr{A}}}^{(m+\mathrm{1,}n)}(t)+n\frac{{{\rm{\Gamma }}}_{b}{\gamma }_{b}}{2}{ {\mathcal L} }_{b}{{\mathscr{A}}}^{(m,n-\mathrm{1)}}(t)+{ {\mathcal L} }_{b}{{\mathscr{A}}}^{(m,n+\mathrm{1)}}(t),\end{array}$$
(35)

where \({{\mathscr{A}}}^{(m,n)}(t)=0\) if m or n < 0. This is the hierarchical equation given in Eq. (15).

Moreover, there is the relation that

$$\begin{array}{rcl}{{\mathscr{A}}}_{ij}^{(m,n)}(t) & = & {\mathcal M} \{{{\mathscr{A}}}_{ij}^{(m,n)}(t,z)\}\\ & = & {\mathcal M} \{{D}_{at}^{m}{D}_{bt}^{n}\langle {{\rm{\Psi }}}_{0}|i\rangle \langle j|(t)|z\rangle \langle z|{{\rm{\Psi }}}_{0}\rangle \}\\ & = & {\mathcal M} \{\langle {{\rm{\Psi }}}_{0}|i\rangle \langle j|(t)[{\xi }_{a}^{\dagger }(t{)]}^{m}{[{\xi }_{b}^{\dagger }(t)]}^{n}|z\rangle \langle z|{{\rm{\Psi }}}_{0}\rangle \}\\ & = & \langle {{\rm{\Psi }}}_{0}|i\rangle \langle j|(t)[{\xi }_{a}^{\dagger }(t{)]}^{m}{[{\xi }_{b}^{\dagger }(t)]}^{n}|{{\rm{\Psi }}}_{0}\rangle \mathrm{.}\end{array}$$

Because \({\xi }_{a}\mathrm{(0)|}{{\rm{\Psi }}}_{0}\rangle =0\), and \({\xi }_{b}\mathrm{(0)|}{{\rm{\Psi }}}_{0}\rangle =0\), then \({{\mathscr{A}}}_{ij}^{(m,n)}\mathrm{(0)}=\langle {{\rm{\Psi }}}_{0}|i\rangle \langle j|{[{\xi }_{a}^{\dagger }\mathrm{(0)]}}^{m}{[{\xi }_{b}^{\dagger }\mathrm{(0)]}}^{n}|{{\rm{\Psi }}}_{0}\rangle =0\), if m or \(n\ne 0\). Therefore, the initial condition of Eq. (35) is \({{\mathscr{A}}}^{(m,n)}\mathrm{(0)}=0\) for m or...

Note that

$${{\mathscr{A}}}^{\mathrm{(0,0)}}(t)= {\mathcal M} \{{\mathscr{A}}(t,z)\}=\prod _{k}\int \frac{{d}^{2}{z}_{ak}{d}^{2}{z}_{bk}}{{\pi }^{2}}{e}^{-|{z}_{ak}{|}^{2}-|{z}_{bk}{|}^{2}}{\mathscr{A}}(t,z)\equiv {\mathscr{A}}(t),$$
(36)

where \({\mathscr{A}}(t)\equiv ({{\mathscr{A}}}_{ij}(t))\), with \(i,j=1,2,3\), and \({{\mathscr{A}}}_{ij}(t)=\langle {{\rm{\Psi }}}_{0}|i\rangle \langle j|(t)|{{\rm{\Psi }}}_{0}\rangle \). From Eq. (35), we get

$$\frac{\partial }{\partial t}{\mathscr{A}}(t)= {\mathcal H} (t){\mathscr{A}}(t)+{ {\mathcal L} }_{a}{{\mathscr{A}}}^{\mathrm{(1,0)}}(t)+{ {\mathcal L} }_{b}{{\mathscr{A}}}^{\mathrm{(0,1)}}(t),$$
(37)

which is the non-Markovian quantum Bloch equation in Eq. (12).