Dark state with counter-rotating dissipative channels

TL;DR

This paper investigates how counter-rotating terms affect the stability of dark states in a three-level quantum system, revealing leakage at zero temperature and proposing a method to restore coherence.

Contribution

It demonstrates the impact of counter-rotating terms on dark state stability beyond the RWA and introduces a leakage elimination technique to recover quantum coherence.

Findings

Counter-rotating terms cause dark state leakage at zero temperature.

Leakage elimination operator can restore dark state coherence.

Dark state stability is compromised without RWA assumptions.

Abstract

Dark state as a consequence of interference between different quantum states has great importance in the fields of chip-scale atomic clock and quantum information. For the $\Lambda$-type three-level system, this dark state is generally regarded as being dissipation-free because it is a superposition of two lowest states without dipole transition between them. However, previous studies are based on the rotating-wave approximation (RWA) by neglecting the counter-rotating terms in the system-environment interaction. In this work, we study non-Markovian quantum dynamics of the dark state in a $\Lambda$-type three-level system coupled to two bosonic baths and reveal the effect of counter-rotating terms on the dark state. In contrast to the dark state within the RWA, leakage of the dark state occurs even at zero temperature, as a result of these counter-rotating terms. Also, we present a method to restore the quantum coherence of the dark state by applying a leakage elimination operator to the system.

Equations

$$H_{\mathrm{sys}}$$

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

$$H_{\mathrm{int}}$$

$$H_{\mathrm{RWA}}=\sum _k(g_{a,k}{a}_k|3⟩⟨1|+g_{b,k}{b}_k|3⟩⟨2|)+H.c..$$

$$\frac{\partial }{\partial t}|D(t)⟩\otimes |0⟩=-i(H_{\mathrm{sys}}+H_{\mathrm{RWA}}+H_{\mathrm{bath}})|D(t)⟩\otimes |0⟩,$$

$${\rho }_{\mathrm{sys}}(t)=\sum _{i,j}{\rho }_{ij}^{(\mathrm{sys})}(t)|i⟩⟨j|,$$

$${\rho }_{ij}^{(\mathrm{sys})}(t)=⟨i|{\mathrm{Tr}}_{\mathrm{env}}\rho (t)|j⟩.$$

$${\rho }_{ij}^{(\mathrm{sys})}(t)$$

$$|j⟩⟨i|(t)\equiv U^{†}(t)|j⟩⟨i|U(t),$$

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

$$\frac{\partial }{\partial t}\mathcal{A}(t)=\mathcal{H}(t)\mathcal{A}(t)+{\mathcal{L}}_a{\mathcal{A}}^{(1,0)}(t)+{\mathcal{L}}_b{\mathcal{A}}^{(0,1)}(t),$$

$$\mathcal{H}(t)=\left(\begin{array}[]{ccccccccc}0&0&0&i\Omega_{1}(t)&0&0&-i\Omega_{1}^{*}(t)&0&0\\ 0&0&0&0&i\Omega_{2}(t)&0&0&-i\Omega_{2}^{*}(t)&0\\ 0&0&0&-i\Omega_{1}(t)&-i\Omega_{2}(t)&0&i\Omega_{1}^{*}(t)&i\Omega_{2}^{*}(t)&0\\ i\Omega_{1}^{*}(t)&0&-i\Omega_{1}^{*}(t)&i\Delta_{31}&0&i\Omega_{2}^{*}(t)&0&0&0\\ 0&i\Omega_{2}^{*}(t)&-i\Omega_{2}^{*}(t)&0&i\Delta_{32}&0&0&0&i\Omega_{1}^{*}(t)\\ 0&0&0&i\Omega_{2}(t)&0&i\Delta_{21}&0&-i\Omega_{1}^{*}(t)&0\\ -i\Omega_{1}(t)&0&i\Omega_{1}(t)&0&0&0&-i\Delta_{31}&0&-i\Omega_{2}(t)\\ 0&-i\Omega_{2}(t)&i\Omega_{2}(t)&0&0&-i\Omega_{1}(t)&0&-i\Delta_{32}&0\\ 0&0&0&0&i\Omega_{1}(t)&0&-i\Omega_{2}^{*}(t)&0&-i\Delta_{21}\\ \end{array}\right),$$

$$\mathcal{L}_{a}=\left(\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}\right),\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \mathcal{L}_{b}=\left(\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}\right).$$

$$\frac{\partial }{\partial t}{\mathcal{A}}^{(m,n)}(t)$$

$$\mathcal{F}(t)\equiv \sqrt{⟨D(t)|{\rho }_{\mathrm{sys}}(t)|D(t)⟩}.$$

$$H_{\mathrm{tot}}=H_{\mathrm{sys}}+H_{\mathrm{bath}}+H_{\mathrm{int}}+R(t),$$

$$R(t)=-c(t)(|1⟩⟨1|+|2⟩⟨2|)$$

$$|B(t)⟩=\frac{{\Omega }_1}{\Omega }{e}^{-i{\omega }_{1}t}|1⟩+\frac{{\Omega }_2}{\Omega }{e}^{-i{\omega }_{2}t}|2⟩.$$

$$H_{\mathrm{int}}$$

$$\frac{d}{dt}{a}_k(t)$$

$$\frac{d}{dt}{b}_k(t)$$

$$a_k(t)$$

$$b_k(t)$$

$$\frac{d}{dt}|1⟩⟨1|(t)$$

$$\frac{d}{dt}|2⟩⟨2|(t)$$

$$\frac{d}{dt}|3⟩⟨3|(t)$$

$$\frac{d}{dt}|3⟩⟨1|(t)$$

$$\frac{d}{dt}|3⟩⟨2|(t)$$

$$\frac{d}{dt}|2⟩⟨1|(t)$$

$$|z⟩\equiv \underset{k}{⨂}|z_k⟩=e^{{\sum }_k{z}_{ak}{a}_k^{†}+{\sum }_k{z}_{bk}{b}_k^{†}}|0⟩,$$

$$a_k|z⟩=z_{ak}|z⟩,\text{leavevmode}\text{leavevmode}\text{leavevmode}{a}_k^{†}|z⟩=\frac{\partial }{\partial z_{ak}}|z⟩,\text{leavevmode}\text{leavevmode}\text{leavevmode}{b}_k|z⟩=z_{bk}|z⟩,\text{leavevmode}\text{leavevmode}\text{leavevmode}{b}_k^{†}|z⟩=\frac{\partial }{\partial z_{bk}}|z⟩.$$

$$\frac{\partial }{\partial t}\mathcal{A}(t,z)$$

$$\frac{\partial }{\partial z_k}=\int ds\frac{\partial z_s}{\partial z_k}\frac{\delta }{\delta z_s},$$

Full text

Introduction

Electromagnetically induced transparency discovered in quantum optics has long been an important effect in physics (see, e.g., [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 dots [2, 3], nitrogen-vacancy center [4] and silicon-vacancy center in diamond [5, 6]. In fact, dark state can have different applications in physics. The atomic clocks based on coherent population trapping [7, 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 state [11, 12, 13]. Operations on the quantum states like squeezing [14] or decay suppression [15, 16] can also be conducted with the help of dark state. Moreover, dark state can have important applications to the slow light [17] and photocell [18].

Here we consider the $\Lambda$-type three-level system where dipole transition between the lowest two states is forbidden. 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 emission [1]. Studies in two-level systems have shown that the counter-rotating terms can change the ground state [19, 20, 21], but there were few studies regarding the influence of the counter-rotating terms on the dark state [22]. 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 $\Lambda$-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 $\Lambda$-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

1 The model Hamiltonian

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

[TABLE]

where $\omega_{1}$, $\omega_{2}$, and $\omega_{3}$ are the three energy levels of the system and the frequencies of the two pumping fields are $\omega_{a}$ and $\omega_{b}$, respectively, with $\Omega_{1}$ and $\Omega_{2}$ characterizing their coupling strengths to the three-level system. To focus on the effect of the counter-rotating dissipative channels, we take $\Omega_{1}$ and $\Omega_{2}$ to be independent of time. Also, the effects like dephasing channels or non-adiabatic transitions are not included here because they were well studied [13]. Here we consider the resonant pumping case with $\omega_{a}=\omega_{3}-\omega_{1}$ and $\omega_{b}=\omega_{3}-\omega_{2}$. This $\Lambda$-type three-level system has a dark state [1],

[TABLE]

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

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^{{\dagger}}_{k}a_{k}+\sum_{k}\omega_{b,k}b^{{\dagger}}_{k}b_{k}$. The interaction between the system and the two baths is

[TABLE]

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 $|1\rangle$ and $|2\rangle$ is forbidden in the considered $\Lambda$-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

[TABLE]

Let the two bosonic baths be in the vacuum state $|0\rangle\equiv|0\rangle_{a}\otimes|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|0\rangle=0$. Thus,

[TABLE]

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 $|3\rangle$, so $|D(t)\rangle\otimes|0\rangle$ is not a solution of the Schrödinger equation, $\frac{\partial}{\partial t}|\Psi\rangle=-iH_{\rm tot}|\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_{\rm 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.

2 Non-Markovian quantum Bloch equation

The reduced density operator of the system can be written as

[TABLE]

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

[TABLE]

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

[TABLE]

where ${\rm Tr}$ denotes the trace over the degrees of freedom of the total system. Thus, $\rho^{\rm(sys)}_{ij}(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

[TABLE]

with

[TABLE]

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,

[TABLE]

with $T$ being the time-ordering operator. Here we study the case with the initial state of the total system given $|\Psi_{0}\rangle\equiv|D(0)\rangle\otimes|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{\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 $\gamma_{i}$, i.e., the inverse of the correlation times of the two baths. The coupling strength between the system and the $i$th bath is given by $\Gamma_{i}$ which corresponds to the decay rate under Markovian approximation [23]. Using the Heisenberg equation, we can derive the following non-Markovian quantum Bloch equation for the expectation values of the system operators (see method):

[TABLE]

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

[TABLE]

[TABLE]

where $\Delta_{ij}=\omega_{i}-\omega_{j}\ (i,j=1,2,3)$ and $\Omega_{i}(t)=\Omega_{i}e^{-i\Delta_{3i}t}\ (i=1,2)$. The non-Markovianity of the quantum dynamics of the three-level system is reflected in both $\mathcal{A}^{(1,0)}(t)$ and $\mathcal{A}^{(0,1)}(t)$, which are solved via the hierarchical equation

[TABLE]

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

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

[TABLE]

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

3 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)=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{\Gamma\gamma}{2}e^{-\gamma|t-s|}$. We also choose the level of state $|3\rangle$ to be the zero point of the energy. The energy difference between $|1\rangle$ and $|3\rangle$ is taken as $\omega_{3}-\omega_{1}=\omega=1$. Other parameters with the frequency unit are expressed as the ratio to $\omega$. To numerically calculate the quantum dynamics of the system, we need to truncate Eq. (LABEL:be2) at a given hierarchical order $\mathcal{N}$ so that $\mathcal{A}^{(m,n)}(t=0)$ for $m+n>\mathcal{N}$. As shown in Fig. 2(a), the numerical results at $\mathcal{N}=4$ (blue curve) are very close to the result at $\mathcal{N}=10$ (read dots), indicating that the results converge rapidly with the hierarchical order $\mathcal{N}$. Note that in this case, we choose $\Gamma=1$, corresponding to a strong coupling between the system and the two baths (the corresponding Markovian decay rate is now comparable to the system frequencies). When $\Gamma<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. (LABEL:be2) at $\mathcal{N}=4$ can already give reliable results for the model that we study here. However, we choose $\mathcal{N}=10$ in our following calculations to have the obtained results more accurate. Also, Fig. 2(b) shows the difference between the results of $\mathcal{N}=10$ and $\mathcal{N}=20$. It is found that the difference is below $10^{-5}$, further indicating that the truncation is reliable.

In Fig. 2(c), we show the fidelity evolution of the dark state by varying the coupling strength $\Gamma$ between the system and the two baths. When $\Gamma=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. However, when decreasing the coupling strength $\Gamma$, the fidelity of the dark state decays slowly and takes a long time to reach the stationary oscillations. Thus, with a given correlation time (inverse of $\gamma$) of the two baths, the dark state has a slow leakage with time when the system-bath coupling is weak. In Fig. 2(d), we also show the effect of the environment correlation time (i.e., $\gamma$) on the fidelity evolution of the dark state. For the non-Markovian environment, which has a longer correlation time (i.e., smaller $\gamma$), 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 $\gamma$). This reveals that the dark state leaks to other levels more slowly when it is coupled to a non-Markovian environment. Moreover, similar to Fig. 2(c), the fidelity of the dark state exhibits stationary oscillations at longer times even in the Markovian limit (large $\gamma$) of the baths. These oscillations are also due to the persistent drive applied to the system.

From our results above, we can conclude that 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.

4 Leakage reduction of the dark state

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

[TABLE]

where

[TABLE]

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 $\tau$ [see Fig. 1(b)]. Within the time intervals $l\tau\leq t<l\tau+\Delta$, where $l$ is a positive integer to identify the pulse is in the $l$th 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 (black and blue 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 black curve has a higher fidelity than the blue 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 Figs. 3(c) and 3(d), we further show the effect of the inverse of the correlation time $\gamma$. 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 $\gamma$. This is the distinct advantage of the non-Markovian baths.

Finally, we discuss more about the leakage elimination operator. Instead of using the eigenstates $|1\rangle$, $|2\rangle$ and $|3\rangle$, we can also use the dark state $|D(t)\rangle$, bright state $|B(t)\rangle$ and eigenstate $|3\rangle$ 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

[TABLE]

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

[TABLE]

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\rangle$. If the state $|3\rangle$ 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. The leakage elimination operator in Eq. (18) also have the effect of counterdiabatic driving [31, 32, 33, 34]. Because it can effectively enlarge the energy difference of the system, the adiabatic condition can be better fulfilled. Here we have proved the equivalence of the leakage elimination operator for $|3\rangle$ to the one for $|D(t)\rangle$ in the case of dissipative channels, but the equivalence about the counterdiabatic effects is not clear.

Discussion

We have studied the non-Markovian quantum dynamics of the dark state in a $\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. 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. Also, people can compare the leakage indicated in our article with the leakage induced by other mechanisms (e.g., dephasing channel, finite temperature effect or non-adiabatic transition) to see whether the effect of counter-rotating terms can be neglected in their cases.

Methods

From the total Hamiltonian $H_{\mathrm{tot}}$ given in Sec. II, it follows that the Heisenberg equations for the field operators $a_{k}$ and $b_{k}$ are given by

[TABLE]

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

[TABLE]

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. (LABEL:sbo) into the Heisenberg equations for the system operators, we have

[TABLE]

with $|1\rangle\langle 3|(t)$, $|2\rangle\langle 3|(t)$, and $|2\rangle\langle 1|(t)$ equal to the conjugates of $|3\rangle\langle 1|(t)$, $|3\rangle\langle 2|(t)$, and $|1\rangle\langle 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,

[TABLE]

which satisfies

[TABLE]

As in [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

[TABLE]

where $\mathcal{A}(t,z)\equiv(\mathcal{A}_{11}(t,z),\mathcal{A}_{22}(t,z),\mathcal{A}_{33}(t,z),\mathcal{A}_{31}(t,z),\mathcal{A}_{32}(t,z),\mathcal{A}_{21}(t,z),\mathcal{A}_{13}(t,z),\mathcal{A}_{23}(t,z),\mathcal{A}_{12}(t,z))^{T}$, with $T$ denoting the transpose of a matrix and $\mathcal{A}_{ij}(t,z)\equiv\langle\Psi_{0}|i\rangle\langle{j}|(t)|z\rangle\langle{z}|\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

[TABLE]

and defined the noise functions,

[TABLE]

Below we solve Eq. (26) using a hierarchical-equation appraoch [36]. 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), $\mathcal{A}(t,z)$ only contains the noise for $0\leq s<t$. Then, for other values of $s$, we have $\frac{\delta}{\delta z_{is}}\mathcal{A}(t,z)=0$, where $i=a,b$. Thus, the extension of the integral limit is exact. Now Eq. (26) can be expressed as

[TABLE]

We further define $\mathcal{A}^{(m,n)}(t,z)\equiv D_{a}^{m}D_{b}^{n}\mathcal{A}(t,z)$ and consider the Ornstein-Uhlenbeck correlation functions $\alpha_{i}(t,s)=\frac{\Gamma_{i}\gamma_{i}}{2}e^{-\gamma_{i}|t-s|}$, with $i=a,b$. Then, it follows that

[TABLE]

Using Eq. (29), we obtain

[TABLE]

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

[TABLE]

This is a stochastic differential equation, because it involves the noise functions $z_{at}$ and $z_{bt}$.

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

[TABLE]

It follows from Eq. (32) that

[TABLE]

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

[TABLE]

because $\xi_{a}(t)|\Psi_{0}\rangle=\sum_{k}g_{a,k}e^{-i\omega_{k}t}a_{k}|D(0)\rangle\otimes|0\rangle=\sum_{k}g_{a,k}e^{-i\omega_{k}t}a_{k}|D(0)\rangle\otimes|0\rangle_{a}\otimes|0\rangle_{b}=0$. Thus, we have $\mathcal{M}\left\{z_{at}\mathcal{L}_{a}\mathcal{A}^{(m,n)}(t,z)\right\}=0$. Similarly, because $\xi_{b}(t)|\Psi_{0}\rangle=\sum_{k}g_{b,k}e^{-i\omega_{k}t}b_{k}|D(0)\rangle\otimes|0\rangle=0$, we can also derive that $\mathcal{M}\left\{z_{bt}\mathcal{L}_{b}\mathcal{A}^{(m,n)}(t,z)\right\}=0$. Therefore, Eq. (34) is simplified as

[TABLE]

where $\mathcal{A}^{(m,n)}(t)=0$ if $m$ or $n<0$. This is the hierarchical equation given in Eq. (LABEL:be2).

Moreover, there is the relation that

[TABLE]

Because $\xi_{a}(0)|\Psi_{0}\rangle=0$, and $\xi_{b}(0)|\Psi_{0}\rangle=0$, then $\mathcal{A}^{(m,n)}_{ij}(0)=\langle\Psi_{0}|i\rangle\langle{j}|[\xi_{a}^{\dagger}(0)]^{m}[\xi_{b}^{\dagger}(0)]^{n}|\Psi_{0}\rangle=0$, if $m$ or $n\neq 0$. Therefore, the initial condition of Eq. (35) is $\mathcal{A}^{(m,n)}(0)=0$ for $m$ or $n\neq 0$.

Note that

[TABLE]

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

[TABLE]

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

Acknowledgements

This work is supported by the National Key Research and Development Program of China (Grant No. 2016YFA0301200), and the NSAF (Grant Nos. U1330201 and U1530401). L.-A.W. is supported by Spanish MINECO/FEDER Grant FIS2015-69983-P, Basque Government Grant IT986-16 and UPV/EHU UFI 11/55. T.Y. is supported by the DOD/AF/AFOSR Grant No. FA9550-12-1-0001, and he thanks CSRC for hospitality during his visit.

Author contributions statement

Z.Y.Z. and M.C. performed the derivations and numerical calculations under the guidance of J.Q.Y.. Also, L.A.W. and T.Y. participated in the discussions. All authors contributed to the interpretation of the work and the writing of the manuscript.

Additional information

Competing financial interests: The authors declare no competing financial interests.