Non-Markovian qubit dynamics in nonequilibrium environments

We consider a single qubit system $S$ with the states $|0\rangle$ and $|1\rangle$ interacting with a nonequilibrium environment $E$. The environmental initial state $\rho_{E}(0)$ is in nonequilibrium corresponding to the nonstationary statistical property of the environment. In the interaction picture with respect to the single qubit system Hamiltonian $H_{S}=(\hbar/2)\omega_{S}\sigma_{z}$ and the environment Hamiltonian $H_{E}$, the Hamiltonian of the total system $S+E$ can be written as Norris et al. (2016)

where $\omega_{S}$ is the frequency difference of the qubit system, $\sigma_{z}$ denotes the Pauli matrix in the single qubit basis $\mathcal{B}_{S}=\{|0\rangle,|1\rangle\}$, and the environmental noise $\epsilon(t)$ is subject to a stochastic process in the presence of a classical environment whereas it is a time dependent operator in the presence of a quantum environment.

The qubit system undergoes pure decoherence during its dynamical evolution due to the fact that $[H_{S},H_{SE}(t)]=0$, and the environmental decoherence effect on the qubit system is reflected in the dynamics of the coherence element. The dynamical evolution of the density matrix of the total system is described by the Liouville equation $d\rho_{SE}(t)/dt=-(i/\hbar)[H_{SE}(t),\rho_{SE}(t)]$. By taking a partial trace over the degrees of freedom of the environment, the reduced density matrix of the single qubit system can be written as

where $\langle\cdots\rangle=\mathrm{Tr}_{E}(\cdots\rho_{E})$ represents a trace over the environmental degrees of freedom with the environment state $\rho_{E}(t)$, $U_{SE}\bm{(}t;\epsilon(t)\bm{)}=\exp\big{[}-(i/\hbar)\int_{0}^{t}dt^{\prime}H_{SE}(t^{\prime})\big{]}$ denotes the unitary time evolution operator in terms of the $H_{SE}(t)$ in Eq. (1) and the single qubit Kraus operators $K_{S\mu}$ satisfy

with the decoherence function $D(t)=e^{-i\Phi(t)-\Gamma(t)}$ in terms of the time-dependent phase angle $\Phi(t)$ and decay factor $\Gamma(t)$ expressed respectively as

based on the statistical characteristics of the environmental noise $\epsilon(t)$. The $n$th-order cumulant $C_{n}(t_{1},\cdots,t_{n})$ depends on the environmental correlation function $\langle\epsilon(t_{1})\cdots\epsilon(t_{k})\rangle$ with $k\leq n$ and $\langle\cdots\rangle$ denoting an ensemble average with respect to the initial state of the environment $\rho_{E}(0)$. The odd-order and even-order cumulants of the environmental noise characterize the degrees of the asymmetry and broadness of the distribution, respectively. It is worth mentioning that the environmental noise $\epsilon(t)$ is stationary if and only if $[H_{E},\rho_{E}(0)]=0$ which corresponds to the environment in a certain equilibrium state initially, and in this case the qubit system undergoes no phase evolution since the odd-order cumulants of the environmental noise vanish Norris et al. (2016). For the general case $[H_{E},\rho_{E}(0)]\neq 0$ corresponding to the environment in an initial nonequilibrium state, the statistical property of the environment is nonstationary and the environmental odd-order correlation functions play an essential role in the unitary dynamical evolution of the qubit system.

Based on the Kraus operators expression for the single qubit reduced density matrix in Eq. (2), the diagonal elements are time independent and the off diagonal elements are time dependent

The dynamics of the single qubit system satisfies a time-local quantum master equation as

where $\phi(t)=d\Phi(t)/dt$ represents the frequency shift caused by the odd-order environmental correlations and $\gamma(t)=d\Gamma(t)/dt$ denotes the decoherence rate induced by the even-order environmental correlations. Due to the environmental decoherence, the diagonal elements of the qubit reduced density matrix do not evolve with time while the off-diagonal elements decay monotonously (Markovian behavior) or non-monotonously (non-Markovian behavior). The non-Markovianity quantifying the exchange of a flow of information between the single qubit system and environment can be, by taking the optimization over all pairs of initial states, expressed as Breuer et al. (2009); Laine et al. (2010)

where $\mathcal{D}(\rho^{1},\rho^{2})=\frac{1}{2}\mathrm{tr}|\rho^{1}-\rho^{2}|$ denotes the trace distance between the states $\rho^{1}$ and $\rho^{2}$ and the optimal pair of initial states is chosen as the maximally coherent states $|\psi_{S}^{\pm}(0)\rangle=(|0\rangle\pm|1\rangle)/\sqrt{2}$. The single qubit dynamics displays non-Markovian behavior if the decoherence rate $\gamma(t)$ takes negative values in some time intervals.

The environmental influences on open quantum systems have been well modeled by RTN in a wide variety of dynamical processes in physics, chemistry and biology, such as, fluorescence process of single molecules Zheng and Brown (2003); Brokmann et al. (2003), decoherence and disentanglement processes in the presence of $1/f^{\alpha}$ noise Burkard (2009); Lo Franco et al. (2014); Castelano et al. (2016) and optical trapping and frequency modulation processes in quantum optics Ashkin et al. (1987); Silveri et al. (2017). Unlike the Gaussian noise of which all cumulants $C_{n}(t_{1},\cdots,t_{n})$ beyond $n=2$ are zero, the generalized RTN displays distinct non-Gaussianity. Its cumulants higher than second-order do not vanish and contribute to the dynamical decoherence of the qubit system. We here consider that the environmental noise $\epsilon(t)$ exhibits with generalized RTN property which is non-Markovian and nonstationary. The noise process transits randomly with an average transition rate $\eta$ between values $\pm 1$ with the amplitude $\chi$ which describes the environmental coupling. Its non-Markovian and nonstationary features are respectively characterized by a generalized master equation for the conditional probability with a memory kernel and by an initially nonstationary distribution. For the generalized RTN process, the environmental non-Markovian feature is described by a generalized master equation for the time evolution of the conditional probability Fuliński (1994)

where $K(t-\tau)$ is the memory kernel of the environmental noise and the conditional probability $P(\varepsilon,t|\varepsilon^{\prime},t^{\prime})$ and transition matrix $T$ are respectively expressed as

We consider the case that the environmental memory kernel is of an exponential form as $K(t-\tau)=\kappa e^{-\kappa(t-\tau)}$ with $\kappa$ denoting the memory decay rate. The smaller the decay rate $\kappa$ is, the stronger the environmental non-Markovian feature is. For the case $\kappa\rightarrow\infty$, namely, the memoryless case $K(t-\tau)=\delta(t-\tau)$, the environmental noise only exhibits Markovian feature.

By means of the Laplace transformation $\tilde{P}(\varepsilon,s|\varepsilon^{\prime},t^{\prime})=\int_{0}^{\infty}P(\varepsilon,t|\varepsilon^{\prime},t^{\prime})e^{-st}dt$, the conditional probability in Eq. (8) can be analytically expressed as

where $I_{2}$ is the $2\times 2$ identity matrix, $\widetilde{K}(s)=\kappa/(s+\kappa)$ represents the Laplace transform of the environmental memory kernel, and $\widetilde{\Omega}(s)=1/[s+2\lambda\widetilde{K}(s)]$. The environmental nonstationary feature arises from the initial distribution

where $a_{0}$ is the nonstationary parameter and $-1\leq a_{0}\leq 1$. For the case $a_{0}=0$, the environment is in equilibrium and the environmental noise only displays stationary feature.

Based on the non-Markovian and nonstationary features, the statistical characteristics of the environmental noise $\epsilon(t)$ are described by the first and second-order moments

where $\mathscr{L}^{-1}$ denotes the inverse Laplace transform. According to the Bayes’ theorem, the environmental higher-order moments satisfy the factorization

for the order of the time instants $t_{1}>\cdot\cdot\cdot>t_{n}$ $(n\geq 3)$.

To obtain the quantum master equation (6) of the qubit system, we need to derive the exact expressions of the phase angle $\Phi(t)$ and decay factor $\Gamma(t)$ based on Eq. (4). However, it is hardly to calculate them for general environmental noise due to the fact that the number of order of the cumulants approaches infinity. It is difficult to derive the exact quantum master equation (first-order differential equation) for reduced density matrix $\rho_{S}(t)$ and only a few physical models can be exactly solved. For the general case, we need to truncate the environmental correlations effects to some finite order within some approximations. The time-ordered factorization rule for the higher-order moments in Eq. (13) results in all the higher-order moments closely associated with the second-order moments, which makes it possible to derive the closed higher-order differential equation for reduced density matrix $\rho_{S}(t)$ to contain all the infinite environmental correlations effects. As an effective methodology, we derive the closed higher-order differential equation for the reduced density matrix $\rho_{S}(t)$. We can solve the closed higher-order differential equation to obtain the exact expression of the qubit reduced density matrix and to derive the exact solutions of the frequency shift $\phi(t)$ and decoherence rate $\gamma(t)$.

In terms of the statistical characteristics of the environmental noise $\epsilon(t)$ in Eqs. (12) and (13), the second-order moment of the environmental noise obeys the closed differential relation

Based on the closed differential relation for the moments in Eq. (14) and the relation between cumulants and moments of the environmental noise

with ${n\choose m}=n!/[(n-m)!m!]$, the dynamical evolution of the reduced density matrix of the single qubit system is governed by a closed third-order differential equation

where the initial conditions satisfy

We have derived the closed third-order differential equation for the single qubit reduced density matrix in the presence of a nonequilibrium environment exhibiting non-Markovian and nonstationary RTN properties. We stress that the reduced density matrix in the closed differential equation (16) contains all the environmental correlations and it is exact without any approximation. The exact solution of the reduced density matrix in Eq. (16) can be analytically expressed as

where $\tilde{D}(s)$ is Laplace transform of the decoherence function $D(t)$ given by

The exact expression of the decoherence function $D(t)$ can be obtained based on the approach we derived in Refs. Cai and Zheng, 2018; Cai, 2020. We here focus on the general case that the denominator of $\tilde{D}(s)$ has three different roots $r_{n}(n=1,2,3)$ and the decoherence function $D(t)$ can be exactly expressed as

where the coefficients satisfy $c_{n}^{r}=\lim\limits_{s\rightarrow r_{n}}\big{[}(s-r_{n})\tilde{D}_{r}(s)\big{]}$ and $c_{n}^{i}=\lim\limits_{s\rightarrow r_{n}}\big{[}(s-r_{n})\tilde{D}_{i}(s)\big{]}$. Correspondingly, the frequency shift and decoherence rate can be respectively written as

In the following, we study the disentanglement dynamics of a two qubit system $T$ consisting of two noninteracting identical single qubits independently coupled to its local nonequilibrium environment. With respect to the two qubit system Hamiltonian $H_{T}=(\hbar/2)\omega_{S}\sigma_{z}\otimes I_{2}+(\hbar/2)\omega_{S}I_{2}\otimes\sigma_{z}$ and the environment Hamiltonian $H_{E}$, the Hamiltonian of the total system $T+E$ in the interaction picture can be expressed as

We construct the reduced density matrix of the two qubit system in the standard product basis $\mathcal{B}_{T}=\{|1\rangle=|11\rangle,|2\rangle=|10\rangle,|3\rangle=|01\rangle,|4\rangle=|00\rangle\}$. Based on the two qubit basis and by taking a partial trace over the environmental degrees of freedom, the reduced density matrix of the two qubit system can be expressed as

where $U_{TE}\bm{(}t;\epsilon(t)\bm{)}=\exp\big{[}-(i/\hbar)\int_{0}^{t}dt^{\prime}H_{TE}(t^{\prime})\big{]}$ denotes the unitary time evolution operator in terms of the Hamiltonian $H_{TE}(t)$ in Eq. (22) and the two qubit Kraus operators $K_{T\mu}(t)=K_{S\nu}(t)\otimes K_{S\upsilon}(t)(\nu,\upsilon=1,2)$ are the tensor products of the single qubit Kraus operators

According to the two qubit Kraus operators expression for the reduced density matrix in Eq. (23), the diagonal elements are time independent

and time-dependent off diagonal elements can be written as

By taking the optimization over all pairs of initial states, the non-Markovianity quantifying the flow of information exchange between the two qubit system and environment can be expressed as Breuer et al. (2009); Laine et al. (2010)

where the optimal pair of initial states can be chosen as the maximally entangled states super-decoherent Bell states $|\psi_{T}^{\pm}(0)\rangle=(|00\rangle\pm|11\rangle)/\sqrt{2}$ or sub-decoherent Bell states $|\varphi_{T}^{\pm}(0)\rangle=(|01\rangle\pm|10\rangle)/\sqrt{2}$ Addis et al. (2013, 2014). It is obvious that similar to the single qubit case, the two qubit dynamics also displays non-Markovian behavior if the decoherence rate $\gamma(t)$ takes negative values sometimes which depends on both the environmental coupling $\chi$ and memory decay rate $\kappa$. Because of the entanglement between the two identical single qubits for the optimal pair of maximally entangled states, the non-Markovianity of the two qubit system is less than twice that of the single qubit system, namely, $\mathcal{N}_{T}<2\mathcal{N}_{S}$.

Due to the environmental effects on its evolution, the two qubit system undergoes dynamical disentanglement. For a two qubit system, all the entanglement measures are compatible. In order to study the nonstationary and non-Markovian effects of the nonequilibrium environment on the disentanglement dynamics of the two qubit system, we here use the concurrence to quantify the entanglement defined as Wootters (1998); Yu and Eberly (2004)

where $\lambda_{i}(t)$ are the eigenvalues of the matrix $\zeta(t)=\rho_{T}(t)(\sigma_{y}\otimes\sigma_{y})\rho_{T}^{*}(t)(\sigma_{y}\otimes\sigma_{y})$ arranged in decreasing order with $\rho_{T}^{*}(t)$ denoting the complex conjugation of the two qubit reduced density matrix $\rho_{T}(t)$ in the two qubit basis $\mathcal{B}_{T}$. The concurrence $\mathscr{C}(t)$ varies from the maximum 1 for a maximally entangled state to the minimum 0 for a completely disentangled state. For pure quantum state, the entanglement corresponds to nonlocal correlations whereas it is not the general case for mixed states due to the fact that the environmental noise gives rise to the decay of nonlocal correlations. The nonlocality can be identified by the violation of the Bell inequalities in the presence of entanglement ($\mathscr{C}(t)>0$). The Clauser-Horne-Shimony-Holt (CHSH) form of Bell function has been widely used to determine whether there are nonlocal correlations of the entangled system. The maximum Bell function $\mathscr{B}(t)$ for an entangled two qubit system can be, based on the Horodecki criterion, expressed as Horodecki et al. (1995)

where the subscripts $j,k=1,2,3$ and $\mu_{j}(t)$ and $\mu_{k}(t)$ are functions in terms of the elements of the two qubit reduced density matrix. If the Bell function $\mathscr{B}(t)$ is larger than the classical threshold $\mathscr{B}_{\mathrm{th}}=2$, the quantum correlations of the entangled two qubit system cannot be reproduced by any classical local model.

It is known that the Bell states and Werner mixed states of a two qubit system play an essential role in quantum computation and quantum information Nielsen and Chuang (2000). The two qubit reduced density matrix expressed in Eq. (23) for initial composite Bell states and extended Werner states has an $X$ structure both initially and during the dynamical evolution. The concurrence $\mathscr{C}(t)$ for an initial $X$ structure reduced density matrix of a two qubit system can be computed in a particular form as Yu and Eberly (2007)

where

The time dependent maximum CHSH-Bell function $\mathscr{B}(t)$ for an X structure two qubit density matrix can be expressed analytically as Derkacz and Jakóbczyk (2005); Mazzola et al. (2010a)

where $\mathscr{B}_{1}(t)=2\sqrt{\mu_{1}(t)+\mu_{2}(t)}$ and $\mathscr{B}_{2}(t)=2\sqrt{\mu_{1}(t)+\mu_{3}(t)}$ with

We first focus on the initial states of the system in the composite Bell states of the form Mazzola et al. (2010b)

where the initial state parameter $c$ is real and satisfies $-1\leq c\leq 1$. It has, by studying the quantum mutual information, quantum discord and classical correlations of the dynamics, shown that for the initial states in Eq. (34) there is a sudden transition from classical to quantum decoherence for the two qubit system coupled to a nonequilibrium environment exhibiting with generalized RTN property, and the nonequilibrium feature of the environment can delay the critical time of the transition of decoherence from classical to quantum Basit et al. (2020). The concurrence at time $t$ in Eq. (28) for the two qubit system prepared in the initial states of Eq. (34) can be reduced to

The initial concurrence of the two qubit system prepared in the composite Bell states in Eq. (34) can be expressed as $\mathscr{C}(0)=2|c|$ since the initial value of the decoherence function satisfies $D(0)=1$. Therefore, the entanglement of the two qubit system exists except for the special case $c=0$. For the case $-1\leq c<0$, the concurrence at time $t$ exists if the threshold value of the decoherence function satisfies $|D(t)|>|D_{\mathrm{th}}|=\sqrt{(1+c)/(1-c)}$ whereas if it exists at time $t$ for the case $0<c\leq 1$, the threshold value of the decoherence function satisfies $|D(t)|>|D_{\mathrm{th}}|=\sqrt{(1-c)/(1+c)}$. The time dependent maximum CHSH-Bell function $\mathscr{B}(t)$ for the initial states of the two qubit system of Eq. (34) can be reduced to

The presence of entanglement $\mathscr{C}(t)>0$ is a necessary condition to achieve nonlocality. The close relation between $\mathscr{B}(t)$ and $\mathscr{C}(t)$ for the two qubit system prepared in the initial composite Bell states of Eq. (34) can be expressed as

For the case $-1\leq c<0$, the classical threshold $\mathscr{C}_{\mathrm{th}}$ which corresponds to the Bell function $\mathscr{B}(t)\geq\mathscr{B}_{\mathrm{th}}=2$ only exists for $-1/2\leq c<0$ and can be expressed as $\mathscr{C}_{\mathrm{th}}=(1-c)\sqrt{1/4-c^{2}}-(1+c)$ whereas for $-1\leq c<-1/2$ the maximum CHSH-Bell function $\mathscr{B}(t)$ is always larger than the threshold $\mathscr{B}_{\mathrm{th}}=2$. Similarly, for the case $0<c\leq 1$, the threshold $\mathscr{C}_{\mathrm{th}}$ for the Bell function larger than the threshold $\mathscr{B}_{\mathrm{th}}=2$ exists for $0<c\leq 1/2$ and can be expressed as $\mathscr{C}_{\mathrm{th}}=(1+c)\sqrt{1/4-c^{2}}-(1-c)$ while the maximum CHSH-Bell function $\mathscr{B}(t)$ is always larger than the threshold $\mathscr{B}_{\mathrm{th}}=2$ for $1/2<c\leq 1$. Therefore, the two qubit system always displays quantum nonlocality at time $t$ for the case $1/2<|c|\leq 1$ whereas for the case $0<|c|\leq 1/2$ the threshold value of the decoherence function should satisfies $|D(t)|>|D_{\mathrm{th}}|=\sqrt[4]{1/2-c^{2}}$ to ensure that the concurrence $\mathscr{C}(t)$ is larger than the classical threshold $\mathscr{C}_{\mathrm{th}}$.

We now focus on the case that the two qubit system is prepared in the extended Werner states initially as Bellomo et al. (2008a); Lo Franco et al. (2014)

where $0\leq r\leq 1$ denotes the purity parameter of the initial states, $I_{4}$ is the $4\times 4$ identity matrix and $|\psi_{T}(0)\rangle=\alpha|00\rangle+\beta|11\rangle$ and $|\varphi_{T}(0)\rangle=\alpha|01\rangle+\beta|10\rangle$ represent the initial entangled states with the complex coefficients $\alpha$, $\beta$ and $|\alpha|^{2}+|\beta|^{2}=1$. For the case $\alpha=\pm\beta=1/\sqrt{2}$, the extended Werner states corresponds to a subclass of Bell-diagonal states, namely, the Werner states Verstraete and Wolf (2002); Horodecki et al. (2009). The concurrence for the two qubit system prepared in the extended Werner states initially of Eq. (38) can be reduced to

The entanglement of the two qubit system exists if the initial value of concurrence $\mathscr{C}(0)$ in the extended Werner states is larger than zero, correspondingly

except for the special cases $\alpha=0$ and $\alpha=1$. The concurrence exists at time $t$ if the threshold value of the decoherence function satisfies $|D(t)|>|D_{\mathrm{th}}|=\sqrt{(1-r)/|\alpha\beta|}/2$. The time dependent maximum CHSH-Bell function $\mathscr{B}(t)$ for the initial extended Werner states of Eq. (38) can be reduced to