Positivity-preserving non-Markovian Master Equation for Open Quantum System Dynamics: Stochastic Schrödinger Equation Approach

In the environmental interaction picture, the interaction Hamiltonian reads:

Using the identity resolution of Bargmann coherent states

the Schrödinger equation, regarding the evolution of the pure state of the composite system $|\Psi_{{\rm tot}}\rangle$, can be rewritten in the Bargmann space representation (setting $\hbar=1$),

Here, we define a complex Gaussian process $z_{t}^{*}\equiv-i\sum_{k}g_{k}z_{k}^{*}e^{i\omega_{k}t}$, which satisfies the following relations: ${\mathcal{M}}(z_{t})={\mathcal{M}}(z_{t}z_{s})=0$, and $\alpha(t,s)\equiv{\mathcal{M}}(z_{t}z_{s}^{*})=\sum_{k}|g_{k}|^{2}e^{-i\omega_{k}(t-s)}$, where $\alpha(t,s)$ is the correlation function of the complex Gaussian process $z_{t}^{*}$. Then, using the chain rule $\frac{\partial(\cdot)}{\partial z_{k}^{*}}=\int_{0}^{t}ds\frac{\partial z_{s}^{*}}{\partial z_{k}^{*}}\frac{\delta(\cdot)}{\delta z_{s}^{*}}$, Eq. (4) can lead to a formal linear non-Markovian QSD equation:

By taking the statistical mean over all trajectories, the reduced density matrix of the system can be recovered

The functional derivative term in Eq. (5) can be formally written using a to-be-determined operator $\hat{O}$, defined as

which can be solved through an operator evolution equation:

where $\hat{H}_{\rm eff}$ is the effective Hamiltonian:

and $\bar{O}(t)\equiv\int_{0}^{t}ds\,\alpha(t,s)\hat{O}(t,s)$. Therefore, the formal linear QSD equation reads:

Generally, the structure of the exact O-operator is complicated. Only a few models can be solved with the exact O-operator. A compromised solution to this difficulty is to replace the exact O-operator with an approximated one. One can drop certain terms of the O-operator to simplify the calculation, called a truncation operation. How to truncate the O-operator depends on the type of interaction and the size of the system. Without the loss of generality, the O-operator can be written as a sum of all component operators [49]:

And the approximated O-operator after the truncation is defined as:

Consequently, the approximated QSD equation after the truncation reads:

where the trajectory $|\psi_{z}^{N}\rangle$ is the associated approximated trajectory.

One of the advantages of the non-Markovian QSD approach is that any reduced density operator recovered from quantum trajectories $\hat{\rho}_{r}={\mathcal{M}}(|\psi_{z}\rangle\langle\psi_{z}|)$ is always positivity preserved, even if quantum trajectories are numerically generated by the approximated QSD equation (13), $\hat{\rho}_{r}^{N}={\mathcal{M}}(|\psi^{N}_{z}\rangle\langle\psi^{N}_{z}|)$. (For the single trajectory, we know that $|\psi_{z}^{N}\rangle\langle\psi_{z}^{N}|$ must be positive semidefinite. The ensemble average ${\mathcal{M}}(|\psi_{z}^{N}\rangle\langle\psi_{z}^{N}|)$ can be considered as a convex combination of $|\psi_{z}^{N}\rangle\langle\psi_{z}^{N}|$, therefore, ${\mathcal{M}}(|\psi_{z}^{N}\rangle\langle\psi_{z}^{N}|)$ is also positive semidefinite. Principally, this is how we derive the positivity-preserving ME from the approximated QSD equation.

For a given exact QSD equation, the associated ME reads:

where $\hat{P}$ denotes the stochastic operator $\hat{P}\equiv|\psi_{z}\rangle\langle\psi_{z}|$. Using the Novikov theorem (see Appendix A), the two terms, ${\mathcal{M}}(z_{t}^{*}\hat{P})$ and ${\mathcal{M}}(z_{t}\hat{P})$, in the above equation, can be estimated,

As a result, the formal ME is obtained,

The above-derived ME is positivity-preserving since the reduced density matrix $\hat{\rho}_{r}$ is rigorously equivalent to the exact stochastic quantum trajectory governed by the QSD equation (5).

Next, we will demonstrate why the ME can not guarantee positivity if all the four exact O-operators in Eq. (16) are replaced by the approximated one $\hat{O}^{N}$. Following the similar method of obtaining Eq. (16), the approximated ME for the perturbative reduced density matrix $\hat{\rho}_{r}^{\prime}$ reads:

However, it is worth pointing out that the reduced density matrix $\hat{\rho}_{r}^{\prime}$ can violate positivity because the approximated ME can not be unraveled by the QSD equation (13). To clarify the difference between Eq. (17) and the approximated ME which is rigorously equivalent to the approximated QSD equation (13), we recover the ME starting from the identity $\hat{\rho}_{r}^{N}={\mathcal{M}}(|\psi_{z}^{N}\rangle\langle\psi_{z}^{N}|)$ and the QSD equation (13). The approximated ME reads

where $\hat{P}^{N}\equiv|\psi_{z}^{N}\rangle\langle\psi_{z}^{N}|$. After applying the Novikov theorem to simplify the term of ${\mathcal{M}}(z_{t}^{*}\hat{P}^{N})$, it is easy to verify that in the general case

This is why the above-mentioned approximated ME (17) cannot be unraveled by the QSD equation (13). To solve this problem, we need to know the exact value of $\delta_{z_{s}}\hat{P}^{N}$, therefore a new O-operator has to be introduced

where the new operator $\hat{O}_{d}(t,s,z^{*})$ is determined by an evolution equation,

together with the initial condition,

The subtle difference between $\hat{O}^{N}$ and $\hat{O}_{d}$ is just the reason of positivity violation in the ME (17). Consequently, the result of applying the Novikov theorem is revised to

where $\bar{O}_{d}^{\dagger}(t,z^{*})\equiv\int_{0}^{t}ds\alpha^{*}(t,s)\hat{O}_{d}^{\dagger}(t,s,z^{*})$. By substituting Eq. (22) into Eq. (18), we obtain the formal positivity-preserving approximated ME,

In Ref. [50, 26], it has been proved that the O-operator of a dissipative three-level system contains noise up to the first order. We use a noise-free operator $\hat{O}^{\{0\}}$ to replace the exact $\hat{O}$ in the QSD equation:

where $\bar{O}^{\{0\}}\equiv\int_{0}^{t}ds\alpha(t,s)\hat{O}^{\{0\}}(t,s)$, and the Lindblad operator $\hat{L}=\hat{J}_{-}$. The operator $\hat{O}^{\{0\}}$ is governed by its evolution equation

with its initial condition

According to Eq. (26), the equation is valid only when the operator $\hat{O}^{\{0\}}$ has the form of

where $f_{1}(t,s)$ and $g_{1}(t,s)$ are two to-be-determined evolution coefficients. By substituting the ansatz (28) into Eq. (26), the coefficients $f_{1}$ and $g_{1}$ can be determined by the following differential equations

associated with the initial conditions,

where $F_{1}(t)\equiv\int_{0}^{t}ds\alpha(t,s)f_{1}(t,s)$, and $G_{1}(t)\equiv\int_{0}^{t}ds\alpha(t,s)g_{1}(t,s)$. Subsequently, the time-evolution equation of the operator $\hat{O}_{d}(t,s,z^{*})$ reads

Note that the last functional derivative term in Eq. (20) has been eliminated because $\hat{O}^{\{0\}}$ is noise-free. Similarly, the ansatz of the operator $\hat{O}_{d}$ reads

By substituting the ansatz (32) into Eq. (31), the new set of coefficients, $f_{2}(t,s),g_{2}(t,s)$ and $p_{2}(t,s,s^{\prime})$, are determined by

with the initial conditions

By comparing Eq. (29) and Eq. (33), it is clear that $f_{1}=f_{2}$ and $g_{1}=g_{2}$, since they have the same evolution equations and the same initial conditions. As a result, $f_{1}$ and $g_{1}$, in the rest of the paper, will be replaced by $f_{2}$ and $g_{2}$, respectively.

After obtaining operators $\hat{O}_{d}$ and $\hat{O}^{\{0\}}$, the formal ME (23) for the dissipative three-level model reads

where $F_{2}(t)$, $G_{2}(t)$, and $P_{2}(t,s^{\prime})$ are defined as

Applying the Novikov theorem and the termination condition in Ref. [22], the term $\hat{J}_{-}^{2}\int_{0}^{t}dsP_{2}(t,s){\mathcal{M}}(z_{s}^{*}\hat{P}^{\{0\}})\hat{J}_{+}$ in the above ME can be further simplified to

Subsequently, the approximated positivity-preserving ME reads,

where the coefficient $P_{f^{*}}(t)\equiv\int_{0}^{t}\int_{0}^{t}dsds^{\prime}P_{2}(t,s)\alpha^{*}(s,s^{\prime})f^{*}_{2}(t,s^{\prime})$. Note that $\alpha(t,s)$ is the time correlation function, corresponding to a variety of spectra, white or colored. For simplicity, we assume the environment is described by an Ornstein-Uhlenbeck process. So, its correlation function is $\alpha(t,s)=a\gamma e^{-\gamma|t-s|}e^{-i\Omega(t-s)}$, where $1/\gamma$ is the scale of memory time and $\Omega$ is the central frequency of the environment. As a result, the coefficients’ evolution equations can be simplified from integrodifferential equations to differential equations that

with the initial conditions,

where $\tilde{P}_{2}(t)\equiv\int_{0}^{t}ds\alpha(t,s)P_{2}(t,s)$.

In this section, we compare the simulation results of the population of states of the three-level system using four different methods: (1) the exact ME, $\hat{\rho}_{r}={\mathcal{M}}(|\psi_{z}\rangle\langle\psi_{z}|)$ in Eq. (16); (2) the approximated positivity-preserving ME, $\hat{\rho}_{r}^{\{0\}}={\mathcal{M}}(|\psi_{z}^{\{0\}}\rangle\langle\psi_{z}^{\{0\}}|)$ in Eq. (23); (3) the approximated non-positivity-preserving ME, $\hat{\rho}_{r}^{\prime}$ in Eq. (17); (4) the approximated QSD, $|\psi_{z}^{\{0\}}\rangle$ in Eq. (25).

First of all, we plot the time evolution of the population of the dissipative three-level system, $\rho_{00}$, $\rho_{11}$, and $\rho_{22}$, generated by approximated QSD equation approach and the approximated positivity-preserving ME approach. The initial state of the system is prepared at an excited state: $|\psi_{z}(t=0)\rangle=|2\rangle$. We choose the frequency, $\omega=1$, and the central frequency of the environment $\Omega=0$. In a strong non-Markovian regime, $\gamma=0.05$, the simulation results from two methods, as shown in Fig. 2, overlap each other. Since the reduced density matrix generated from the approximated QSD approach is naturally positivity-preserving, the matched dynamics prove that our approximated ME gives rise to the same degree of accuracy as the approximated QSD equation. It indicates that the approximated ME can guarantee positivity.

Next, we plot the time evolution of the population of the ground state using three different ME approaches in Fig. 3. Using the same parameters as Fig. 2, we observe that the exact ME and our approximated positivity-preserving ME both preserve the positivity. Meanwhile, the simulation result of the non-positivity-preserving ME leads to failure due to the appearance of negative probabilities in some time intervals. Furthermore, the magnitude of the negative probability increases with time up to infinity. Consequently, the probabilities of the other two levels also increase to infinity simultaneously. If simply replaces the exact O-operator in the exact ME with the truncated operator $\hat{O}^{\{0\}}$, then the Eq. (17) can be explicitly written as

It is clear that the above approximated ME does not preserve positivity in some parameter regions, and may lead to meaningless physics.

To further demonstrate the importance of our method in studying non-Markovian dynamics, we plot Fig. 4, the time evolution of the population of the ground state in a moderate non-Markovian regime. When $\gamma=0.2$, a shorter memory time compared with the parameter $\gamma=0.05$ used in Fig. 3, the dynamics simulated by the approximated non-positivity-preserving ME $\rho_{r,00}^{\prime}$ do not contain any negative probabilities. However, its distance from the results of the exact ME approach is significantly larger compared with the results of our positivity-preserving ME. Comparing Figs. 3 and 4, we show that the reduced density matrix $\rho_{r}^{\{0\}}$ can guarantee positivity preservation from the Markovian to the strong non-Markovian regime. In contrast, the reduced density matrix $\rho_{r}^{\prime}$ cannot offer such confidence. Moreover, for different models and interested parameter spaces, our method is flexible for different approximations. It provides a robust method to obtain positivity-preserving MEs for the analysis of non-Markovian dynamics.