On the structure of the master equation for a two-level system coupled to a thermal bath

In order to derive a master equation, we need to compute ${\rm d}\mathcal{M}[P]/{\rm d}t$, where for simplicity in the notation we have renamed $P={{P}_{\xi ,\nu }}$. With this aim, we take into account that ${\rm d}\mathcal{M}[P]/{\rm d}t=\mathcal{M}[{\rm d}P/{\rm d}t]$ and consider the SLN equation (2) for P. In order to have a more compact notation, we also re-express the equation (2) in the interaction picture with respect to the system and in a Liouvillian form

$$\begin{eqnarray}&&\frac{{\rm d}P_{t}^{v}}{{\rm d}t}=\mathcal{L}(t)P_{t}^{v},\end{eqnarray} \tag{ 12 }$$

where $\mathcal{L}(t)$ is a super-operator acting over the stochastic projector P_t flattened as a vector, denoted as P^v_t. In detail, $\mathcal{L}(t)P_{t}^{v}={{\left( \frac{{\rm i}}{\hbar }\xi (t)[X(t),{{P}_{\xi ,\nu }}]+\frac{{\rm i}}{2}\nu (t)\{X(t),{{P}_{\xi ,\nu }}\} \right)}^{v}}$, where ${{\left( \cdots \right)}^{v}}$ denotes the vector form, and $X(t)={{{\rm e}}^{{\rm i}{{H}_{S}}t}}\hat{x}{{{\rm e}}^{-{\rm i}{{H}_{S}}t}}$. The equation above can be re-written as

$$\begin{eqnarray}\begin{array}{rcl} \frac{{\rm d}P_{t}^{v}}{{\rm d}t} & = & \mathop{\sum }\limits_{j=0,1}\mathop{\sum }\limits_{\lambda }L_{j\lambda }^{0}(t){{z}_{j\lambda }}P_{t}^{v} \\ {} & {} & +\mathop{\sum }\limits_{\lambda }L_{j\lambda }^{*}(t)z_{j\lambda }^{*}P_{t}^{v}, \\ \end{array}\end{eqnarray} \tag{ 13 }$$

where $L_{j\lambda }^{0}(t)=\partial \mathcal{L}(t)/\partial {{z}_{j\lambda }}$, and $L_{j\lambda }^{*}(t)=\partial \mathcal{L}(t)/\partial z_{j\lambda }^{*}$ are new super-operators. Note that, because of the interaction picture, these quantities are rotated with respect to H_S. In addition, they no longer depend on ${{z}_{j\lambda }}$ or $z_{j\lambda }^{*}$. The analytical average of the last equation leads to the evolution of the reduced density matrix of the system as

$$\begin{eqnarray}\begin{array}{rcl} \frac{{\rm d}\rho _{s}^{v}(t)}{{\rm d}t}=\mathcal{M}[\frac{{\rm d}P_{t}^{v}}{{\rm d}t}] & = & \mathop{\sum }\limits_{j=0,1}\mathop{\sum }\limits_{\lambda }L_{j\lambda }^{0}(t)\mathcal{M}[{{z}_{j\lambda }}P_{t}^{v}] \\ {} & {} & +\mathop{\sum }\limits_{\lambda }L_{j\lambda }^{*}(t)\mathcal{M}[z_{j\lambda }^{*}P_{t}^{v}]. \\ \end{array}\end{eqnarray} \tag{ 14 }$$

Unfortunately this is an open equation, in the sense that the derivative of ${{\rho }_{s}}(t)$ does not depend only on ${{\rho }_{s}}$ as expected for a proper master equation, but it also depends on the unknown quantities $\mathcal{M}[{{z}_{j\lambda }}P_{t}^{v}]$ and $\mathcal{M}[z_{j\lambda }^{*}P_{t}^{v}]$. However, as noted above, we now know that $\mathcal{M}[\cdots ]=\int {\rm d}\mu ({{z}_{0}})\int {\rm d}\mu ({{z}_{1}})\cdots$, with ${\rm d}\mu ({{z}_{j}})$ given by (9). Then, we can re-express

$$\begin{eqnarray}\begin{array}{rcl} \mathcal{M}[{{z}_{j\lambda }}P_{t}^{v}] & = & \mathcal{M}\left[ \frac{\partial P_{t}^{v}}{\partial z_{j\lambda }^{*}} \right]=A_{\lambda }^{*}, \\ \mathcal{M}[z_{j\lambda }^{*}P_{t}^{v}] & = & \mathcal{M}\left[ \frac{\partial P_{t}^{v}}{\partial {{z}_{j\lambda }}} \right]=A_{\lambda }^{0}, \\ \end{array}\end{eqnarray} \tag{ 15 }$$

where the right-hand side of the expressions have been obtained by considering an integration by parts of the left-hand side. As one can see from these expressions, the evolution of P^v_t is intimately related to the evolution of $\partial P_{t}^{v}/\partial z_{j\lambda }^{*}$ and $\partial P_{t}^{v}/\partial {{z}_{j\lambda }}$. Hence, the evolution equation for $\rho _{s}^{v}$ is written in the form

$$\begin{eqnarray}&&\frac{{\rm d}\rho _{s}^{v}(t)}{{\rm d}t}=\mathop{\sum }\limits_{j=0,1}\mathop{\sum }\limits_{\lambda }L_{j\lambda }^{0}(t)A_{\lambda }^{*}+\mathop{\sum }\limits_{\lambda }L_{j\lambda }^{*}(t)A_{\lambda }^{0}.\end{eqnarray} \tag{ 16 }$$

This is still an open equation that depends on two unknown quantities $A_{\lambda }^{*}(t)$ and $A_{\lambda }^{0}(t)$. To proceed further, we need to re-express these quantities in terms of $\rho _{s}^{v}(t)$. The first step is to calculate their evolution equation, taking into account the consistency conditions $({\rm d}/{\rm d}t)(\partial P/\partial {{z}_{j\lambda }})=(\partial /\partial {{z}_{j\lambda }})({\rm d}P/{\rm d}t)$, and $({\rm d}/{\rm d}t)(\partial P/\partial z_{j\lambda }^{*})=(\partial /\partial z_{j\lambda }^{*})({\rm d}P/{\rm d}t)$ (see appendix for details), together with equation (13)

$$\begin{eqnarray}\begin{array}{rcl} \frac{\partial }{\partial {{z}_{j\lambda }}}\frac{{\rm d}P_{t}^{v}}{{\rm d}t} & = & L_{j\lambda }^{0}(t)P_{t}^{v}+\mathop{\sum }\limits_{{{j}^{\prime} }=0,1}\mathop{\sum }\limits_{\lambda ^{\prime} }L_{{{j}^{\prime} }\lambda ^{\prime} }^{0}(t){{z}_{{{j}^{\prime} }\lambda ^{\prime} }}\frac{\partial P_{t}^{v}}{\partial {{z}_{j\lambda }}} \\ {} & {} & +\mathop{\sum }\limits_{\lambda ^{\prime} }L_{{{j}^{\prime} }\lambda ^{\prime} }^{*}(t)z_{{{j}^{\prime} }\lambda ^{\prime} }^{*}\frac{\partial P_{t}^{v}}{\partial {{z}_{j\lambda }}} \\ \frac{\partial }{\partial z_{j\lambda }^{*}}\frac{{\rm d}P_{t}^{v}}{{\rm d}t} & = & L_{j\lambda }^{*}(t)P_{t}^{v}+\mathop{\sum }\limits_{j=0,1}\mathop{\sum }\limits_{\lambda ^{\prime} }L_{j\lambda ^{\prime} }^{0}(t){{z}_{{{j}^{\prime} }\lambda ^{\prime} }}\frac{\partial P_{t}^{v}}{\partial z_{j\lambda }^{*}} \\ {} & {} & +\mathop{\sum }\limits_{\lambda ^{\prime} }L_{{{j}^{\prime} }\lambda ^{\prime} }^{*}(t)z_{{{j}^{\prime} }\lambda ^{\prime} }^{*}\frac{\partial P_{t}^{v}}{\partial z_{j\lambda }^{*}}. \\ \end{array}\end{eqnarray} \tag{ 17 }$$

If we now compute the average $\mathcal{M}[\cdots ]$ of the latter equation, we will find

$$\begin{eqnarray}\begin{array}{rcl} \frac{{\rm d}A_{\lambda }^{0}}{{\rm d}t} & = & L_{j\lambda }^{0}(t)\rho _{s}^{v}(t)+\mathop{\sum }\limits_{{{j}^{\prime} }=0,1}\mathop{\sum }\limits_{\lambda ^{\prime} }L_{{{j}^{\prime} }\lambda ^{\prime} }^{0}(t)B_{\lambda ^{\prime} \lambda }^{*0} \\ {} & {} & +\mathop{\sum }\limits_{{{j}^{\prime} }=0,1}\mathop{\sum }\limits_{\lambda ^{\prime} }L_{{{j}^{\prime} }\lambda ^{\prime} }^{*}(t)B_{\lambda ^{\prime} \lambda }^{00} \\ \frac{{\rm d}A_{\lambda }^{*}}{{\rm d}t} & = & L_{j\lambda }^{*}(t)\rho _{s}^{v}(t)+\mathop{\sum }\limits_{{{j}^{\prime} }=0,1}\mathop{\sum }\limits_{\lambda ^{\prime} }L_{{{j}^{\prime} }\lambda ^{\prime} }^{0}(t)B_{\lambda ^{\prime} \lambda }^{**} \\ {} & {} & +\mathop{\sum }\limits_{{{j}^{\prime} }=0,1}\mathop{\sum }\limits_{\lambda ^{\prime} }L_{{{j}^{\prime} }\lambda ^{\prime} }^{*}(t)B_{\lambda ^{\prime} \lambda }^{0*}. \\ \end{array}\end{eqnarray} \tag{ 18 }$$

which depends on terms of the form

$$\begin{eqnarray}\begin{array}{rcl} \mathcal{M}[{{z}_{{{j}^{\prime} }\lambda ^{\prime} }}\frac{\partial P_{t}^{v}}{\partial {{z}_{j\lambda }}}] & = & \mathcal{M}\left[ \frac{\partial }{\partial z_{{{j}^{\prime} }\lambda ^{\prime} }^{*}}\frac{\partial P_{t}^{v}}{\partial {{z}_{j\lambda }}} \right]=B_{\lambda ^{\prime} \lambda }^{*0}, \\ \mathcal{M}[{{z}_{{{j}^{\prime} }\lambda ^{\prime} }}\frac{\partial P_{t}^{v}}{\partial z_{j\lambda }^{*}}] & = & \mathcal{M}\left[ \frac{\partial }{\partial z_{{{j}^{\prime} }\lambda ^{\prime} }^{*}}\frac{\partial P_{t}^{v}}{\partial z_{j\lambda }^{*}} \right]=B_{\lambda ^{\prime} \lambda }^{**}, \\ \mathcal{M}[z_{{{j}^{\prime} }\lambda ^{\prime} }^{*}\frac{\partial P_{t}^{v}}{\partial {{z}_{j\lambda }}}] & = & \mathcal{M}\left[ \frac{\partial }{\partial {{z}_{{{j}^{\prime} }\lambda ^{\prime} }}}\frac{\partial P_{t}^{v}}{\partial {{z}_{j\lambda }}} \right]=B_{\lambda ^{\prime} \lambda }^{00}, \\ \mathcal{M}[z_{{{j}^{\prime} }\lambda ^{\prime} }^{*}\frac{\partial P_{t}^{v}}{\partial z_{j\lambda }^{*}}] & = & \mathcal{M}\left[ \frac{\partial }{\partial {{z}_{{{j}^{\prime} }\lambda ^{\prime} }}}\frac{\partial P_{t}^{v}}{\partial z_{j\lambda }^{*}} \right]=B_{\lambda ^{\prime} \lambda }^{0*}. \\ \end{array}\end{eqnarray} \tag{ 19 }$$

For the sake of simplicity, we note that, in the quantities $A_{\lambda }^{\beta }$ and $B_{\lambda {{\lambda }^{\prime} }}^{\beta {{\beta }^{\prime} }}$ ($\beta ,\beta ^{\prime} =*,0$), we skip the indexes j and ${{j}^{\prime} }j$ respectively. To reach an even higher order, we need to compute the evolution of the four quantities $(\partial /\partial z_{j\lambda ^{\prime} }^{*})(\partial P_{t}^{v}/\partial {{z}_{j\lambda }})$, $(\partial /\partial z_{j\lambda ^{\prime} }^{*})(\partial P_{t}^{v}/\partial z_{j\lambda }^{*})$, $(\partial /\partial {{z}_{j\lambda ^{\prime} }})(\partial P_{t}^{v}/\partial {{z}_{j\lambda }})$, and $(\partial /\partial {{z}_{j\lambda ^{\prime} }})(\partial P_{t}^{v}/\partial z_{j\lambda }^{*})$. At this point, we consider the above introduced consistency condition, and the evolution equation (17), to get the following

$$\begin{eqnarray}\begin{array}{rcl} \frac{{{\partial }^{2}}}{\partial {{z}_{{{j}^{\prime }}\lambda ^{\prime} }}\partial {{z}_{j\lambda }}}\frac{{\rm d}P}{{\rm d}t} & = & L_{j\lambda }^{0}(t)\frac{\partial P_{t}^{v}}{\partial {{z}_{{{j}^{\prime }}\lambda ^{\prime} }}}+L_{{{j}^{\prime }}\lambda ^{\prime} }^{0}(t)\frac{\partial P_{t}^{v}}{\partial {{z}_{j\lambda }}} \\ {} & {} & +\mathop{\sum }\limits_{{{j}^{\prime }}=0,1}\mathop{\sum }\limits_{\lambda ^{\prime} ^{\prime} }\left( L_{{{j}^{\prime }}\lambda ^{\prime} ^{\prime} }^{0}(t){{z}_{{{j}^{\prime }}\lambda ^{\prime} ^{\prime} }}+L_{{{j}^{\prime }}\lambda ^{\prime} ^{\prime} }^{*}(t)z_{{{j}^{\prime }}\lambda ^{\prime} }^{*} \right)\frac{{{\partial }^{2}}P_{t}^{v}}{\partial {{z}_{{{j}^{\prime }}\lambda ^{\prime} }}\partial {{z}_{j\lambda }}} \\ \frac{{{\partial }^{2}}}{\partial z_{{{j}^{\prime }}\lambda ^{\prime} }^{*}\partial {{z}_{j\lambda }}}\frac{{\rm d}P}{{\rm d}t} & = & L_{j\lambda }^{0}(t)\frac{\partial P_{t}^{v}}{\partial z_{{{j}^{\prime }}\lambda ^{\prime} }^{*}}+L_{{{j}^{\prime }}\lambda ^{\prime} }^{*}(t)\frac{\partial P_{t}^{v}}{\partial {{z}_{j\lambda }}} \\ {} & {} & +\mathop{\sum }\limits_{{{j}^{\prime }}=0,1}\mathop{\sum }\limits_{\lambda ^{\prime} ^{\prime} }\left( L_{{{j}^{\prime }}\lambda ^{\prime} ^{\prime} }^{0}(t){{z}_{{{j}^{\prime }}\lambda ^{\prime} ^{\prime} }}+L_{{{j}^{\prime }}\lambda ^{\prime} ^{\prime} }^{*}(t)z_{{{j}^{\prime }}\lambda ^{\prime} ^{\prime} }^{*} \right)\frac{{{\partial }^{2}}P_{t}^{v}}{\partial z_{{{j}^{\prime }}\lambda ^{\prime} }^{*}\partial {{z}_{j\lambda }}} \\ \frac{{{\partial }^{2}}}{\partial {{z}_{{{j}^{\prime }}\lambda ^{\prime} }}\partial z_{j\lambda }^{*}}\frac{{\rm d}P_{t}^{v}}{{\rm d}t} & = & L_{j\lambda }^{*}(t)\frac{\partial P_{t}^{v}}{\partial {{z}_{{{j}^{\prime }}\lambda ^{\prime} }}}+L_{{{j}^{\prime }}\lambda ^{\prime} }^{0}(t)\frac{\partial P_{t}^{v}}{\partial z_{j\lambda }^{*}} \\ {} & {} & +\mathop{\sum }\limits_{{{j}^{\prime }}=0,1}\mathop{\sum }\limits_{\lambda ^{\prime} ^{\prime} }\left( L_{{{j}^{\prime }}\lambda ^{\prime} ^{\prime} }^{0}(t){{z}_{{{j}^{\prime }}\lambda ^{\prime} ^{\prime} }}+L_{{{j}^{\prime }}\lambda ^{\prime} ^{\prime} }^{*}(t)z_{{{j}^{\prime }}\lambda ^{\prime} ^{\prime} }^{*} \right)\frac{{{\partial }^{2}}P_{t}^{v}}{\partial {{z}_{{{j}^{\prime }}\lambda ^{\prime} }}\partial z_{j\lambda }^{*}} \\ \frac{{{\partial }^{2}}}{\partial z_{{{j}^{\prime }}\lambda ^{\prime} }^{*}\partial z_{j\lambda }^{*}}\frac{{\rm d}P_{t}^{v}}{{\rm d}t} & = & L_{j\lambda }^{*}(t)\frac{\partial P_{t}^{v}}{\partial z_{{{j}^{\prime }}\lambda ^{\prime} }^{*}}+L_{{{j}^{\prime }}\lambda ^{\prime} }^{*}(t)\frac{\partial P_{t}^{v}}{\partial z_{j\lambda }^{*}} \\ {} & {} & +\mathop{\sum }\limits_{{{j}^{\prime }}=0,1}\mathop{\sum }\limits_{\lambda ^{\prime} ^{\prime} }\left( L_{{{j}^{\prime }}\lambda ^{\prime} ^{\prime} }^{0}(t){{z}_{{{j}^{\prime }}\lambda ^{\prime} ^{\prime} }}+L_{{{j}^{\prime }}\lambda ^{\prime} ^{\prime} }^{*}(t)z_{{{j}^{\prime }}\lambda ^{\prime} ^{\prime} }^{*} \right)\frac{{{\partial }^{2}}P_{t}^{v}}{\partial z_{{{j}^{\prime }}\lambda ^{\prime} }^{*}\partial z_{j\lambda }^{*}}. \\ \end{array}\end{eqnarray} \tag{ 20 }$$

Computing the stochastic average of the latter equations gives rise to the following set of equations

$$\begin{eqnarray}\begin{array}{rcl} \frac{{\rm d}B_{\lambda ^{\prime} \lambda }^{00}}{{\rm d}t} & = & L_{j\lambda }^{0}(t)A_{\lambda ^{\prime} }^{0}+L_{{{j}^{\prime }}\lambda ^{\prime} }^{*}(t)A_{\lambda }^{0} \\ {} & {} & +\mathop{\sum }\limits_{{{j}^{\prime }}=0,1}\mathop{\sum }\limits_{\lambda ^{\prime} }\left( L_{{{j}^{\prime }}\lambda ^{\prime} }^{0}(t)C_{\lambda ^{\prime} ^{\prime} \lambda ^{\prime} \lambda }^{*00}+L_{{{j}^{\prime }}\lambda ^{\prime} }^{*}(t)C_{\lambda ^{\prime} ^{\prime} \lambda ^{\prime} \lambda }^{000} \right) \\ \frac{{\rm d}B_{\lambda ^{\prime} \lambda }^{*0}}{{\rm d}t} & = & L_{j\lambda }^{0}(t)A_{\lambda ^{\prime} }^{*}+L_{{{j}^{\prime }}\lambda ^{\prime} }^{0}(t)A_{\lambda }^{0} \\ {} & {} & +\mathop{\sum }\limits_{{{j}^{\prime }}=0,1}\mathop{\sum }\limits_{\lambda ^{\prime} }\left( L_{{{j}^{\prime }}\lambda ^{\prime} }^{0}(t)C_{\lambda ^{\prime} ^{\prime} \lambda ^{\prime} \lambda }^{**0}+L_{{{j}^{\prime }}\lambda ^{\prime} }^{*}(t)C_{\lambda ^{\prime} ^{\prime} \lambda ^{\prime} \lambda }^{0*0} \right) \\ \frac{{\rm d}B_{\lambda ^{\prime} \lambda }^{**}}{{\rm d}t} & = & L_{j\lambda }^{*}(t)A_{\lambda ^{\prime} }^{*}+L_{{{j}^{\prime }}\lambda ^{\prime} }^{*}(t)A_{\lambda }^{*} \\ {} & {} & +\mathop{\sum }\limits_{{{j}^{\prime }}=0,1}\mathop{\sum }\limits_{\lambda ^{\prime} }\left( L_{{{j}^{\prime }}\lambda ^{\prime} }^{0}(t)C_{\lambda ^{\prime} ^{\prime} \lambda ^{\prime} \lambda }^{***}+L_{{{j}^{\prime }}\lambda ^{\prime} }^{*}(t)C_{\lambda ^{\prime} ^{\prime} \lambda ^{\prime} \lambda }^{0**} \right) \\ \frac{{\rm d}B_{\lambda ^{\prime} \lambda }^{0*}}{{\rm d}t} & = & L_{j\lambda }^{*}(t)A_{\lambda ^{\prime} }^{0}+L_{{{j}^{\prime }}\lambda ^{\prime} }^{0}(t)A_{\lambda }^{*} \\ {} & {} & +\mathop{\sum }\limits_{{{j}^{\prime }}=0,1}\mathop{\sum }\limits_{\lambda ^{\prime} }\left( L_{{{j}^{\prime }}\lambda ^{\prime} }^{0}(t)C_{\lambda ^{\prime} ^{\prime} \lambda ^{\prime} \lambda }^{*0*}+L_{{{j}^{\prime }}\lambda ^{\prime} }^{*}(t)C_{\lambda ^{\prime} ^{\prime} \lambda ^{\prime} \lambda }^{00*} \right), \\ \end{array}\end{eqnarray} \tag{ 21 }$$

where the $C_{\lambda ^{\prime\prime} \lambda ^{\prime} \lambda }^{\alpha \beta \gamma }$, with $\alpha ,\beta ,\gamma =0,*$ are given by partial derivatives of the projector over $z_{{{j}^{{\prime\prime} }}{{\lambda }^{{\prime\prime} }}}^{\alpha }$, $z_{{{j}^{\prime} }\lambda ^{\prime} }^{\beta }$ and $z_{j\lambda }^{\gamma }$.

The equations (16), (18) and (21) form an open system of mutually dependent equations. Because of this structure, the evolution of $\rho _{S}^{v}(t)$ can be written in terms of a complex series of terms that only depend on $\rho _{S}^{v}$ at past times. The first two terms of the series can be obtained as follows:

• (i)  
  First integrate analytically (18), obtaining a formal solution for $A_{\lambda }^{0}(t)$, and $A_{\lambda }^{*}(t)$, which comes in terms of different Bs

  $$\begin{eqnarray}\begin{array}{rcl} A_{\lambda }^{0}(t) & = & A_{\lambda }^{0}(0)+\int _{0}^{t}{\rm d}s[L_{j\lambda }^{0}(s)\rho _{s}^{v}(s) \\ {} & {} & \left. +\mathop{\sum }\limits_{{{j}^{\prime} }=0,1}\mathop{\sum }\limits_{\lambda ^{\prime} }\left( L_{{{j}^{\prime} }\lambda ^{\prime} }^{0}(s)B_{\lambda ^{\prime} \lambda }^{*0}(s)+L_{{{j}^{\prime} }\lambda ^{\prime} }^{*}(s)B_{\lambda ^{\prime} \lambda }^{00} \right) \right] \\ A_{\lambda }^{*}(t) & = & A_{\lambda }^{*}(0)+\int _{0}^{t}{\rm d}s[L_{j\lambda }^{*}(s)\rho _{s}^{v}(s) \\ {} & {} & \left. +\mathop{\sum }\limits_{{{j}^{\prime} }=0,1}\mathop{\sum }\limits_{\lambda ^{\prime} }\left( L_{{{j}^{\prime} }\lambda ^{\prime} }^{0}(s)B_{\lambda ^{\prime} \lambda }^{**}(s)+L_{{{j}^{\prime} }\lambda ^{\prime} }^{*}(s)B_{\lambda ^{\prime} \lambda }^{0*}(s) \right) \right]. \\ \end{array}\end{eqnarray} \tag{ 22 }$$

  Note that all initial states of A, B, C, etc are zero;
• (ii)  
  In a similar way, compute the different Bs by analytically integrating (21). From this, an equation here labelled as (21b) is obtained, that expresses the Bs in terms of As and Cs;
• (iii)  
  Insert the values of A given by (22) into (21b), which gives rise to an equation here labelled as (21c).
• (iv)  
  Since (21c) will come again in terms of Bs, replace on its right-hand side the expression of B given by (21b), obtaining (21d).
• (v)  
  In the resulting equation (21d) replace again the values of A given by (22), as well as the values of B given by (21b).

Equations (21a, b, c) are not written here for simplicity in the presentation, but are straight forwardly derivable following the above steps. The process shall be repeated until an equation for the As is obtained, which comes only in terms of $\rho _{s}^{v}$ at different times, and higher order terms $C_{\lambda ^{\prime\prime} \lambda ^{\prime} \lambda }^{\alpha \beta \gamma }$. Neglecting these higher order terms, the following master equation is obtained

$$\begin{eqnarray}\begin{array}{rcl} \frac{{\rm d}\rho _{s}^{v}(t)}{{\rm d}t} & = & \mathop{\sum }\limits_{j=0,1}\int _{0}^{t}{\rm d}s\{\alpha _{jj}^{0*}(t-s)\mathcal{L}_{j}^{0}(t)\mathcal{L}_{j}^{*}(s) \\ {} & {} & \left. +\int _{0}^{t}{\rm d}s\alpha _{jj}^{*0}(t-s)\mathcal{L}_{j}^{*}(t)\mathcal{L}_{j}^{0}(s) \right\}{{\rho }^{v}}(s) \\ {} & {} & +\mathop{\sum }\limits_{\frac{\xi ,\beta ,\gamma ,\eta }{\xi \ne \beta ,\gamma \ne \eta }}\int _{0}^{t}{\rm d}s\int _{0}^{s}{\rm d}{{s}^{\prime} }\int _{0}^{{{s}^{\prime} }}{\rm d}{{s}^{{\prime\prime} }}\alpha _{jj}^{\xi \beta }(t-{{s}^{\prime} })\alpha _{jj}^{\gamma \eta }(s-{{s}^{{\prime\prime} }}) \\ {} & {} & \times \mathcal{U}_{j{{j}^{\prime} }j{{j}^{\prime} }}^{\xi ,\beta ,\gamma ,\eta }(t,s,{{s}^{\prime} },{{s}^{{\prime\prime} }})\rho _{s}^{v}(s) \\ {} & {} & +\mathop{\sum }\limits_{\frac{\xi ,\beta ,\gamma ,\eta }{\xi \ne \beta ,\gamma \ne \eta }}\int _{0}^{t}{\rm d}s\int _{0}^{s}{\rm d}{{s}^{\prime} }\int _{0}^{{{s}^{\prime} }}{\rm d}{{s}^{{\prime\prime} }}\alpha _{jj}^{\xi \beta }(t-{{s}^{{\prime\prime} }})\alpha _{jj}^{\gamma \eta }(s-{{s}^{\prime} }) \\ {} & {} & \times \mathcal{U}_{j{{j}^{\prime} }{{j}^{\prime} }j}^{\xi ,\gamma ,\eta ,\beta }(t,s,{{s}^{\prime} },{{s}^{{\prime\prime} }})\rho _{s}^{v}(s)+\cdots , \\ \end{array}\end{eqnarray} \tag{ 23 }$$

where $\xi ,\beta ,\eta ,\nu =0,*$, and $j,i,n,m=0,1$, and ${{\rho }^{v}}(t)$ is the reduced density matrix of the system flattened as a vector. Note that because of the restrictions of the sums, we are left with terms proportional to $\alpha _{jj}^{0*}\alpha _{jj}^{0*}$, $\alpha _{jj}^{0*}\alpha _{jj}^{*0}$, $\alpha _{jj}^{*0}\alpha _{jj}^{*0}$, or $\alpha _{jj}^{*0}\alpha _{jj}^{0*}$. We have also defined $\mathcal{U}_{jinm}^{\alpha ,\beta ,\gamma ,\eta }(t,s,{{s}^{\prime} },{{s}^{{\prime\prime} }})=\mathcal{L}_{j}^{\alpha }(t)\mathcal{L}_{i}^{\beta }(s)\mathcal{L}_{n}^{\gamma }({{s}^{\prime} })\mathcal{L}_{m}^{\eta }({{s}^{{\prime\prime} }})$.The super-operators appearing in this expression are $\mathcal{L}_{0}^{0}(t)=\mathcal{L}_{0}^{*}(t)={{\mathcal{L}}_{M}}(t)$, $\mathcal{L}_{1}^{0}(t)=2{{\mathcal{L}}_{M}}(t)$, and $\mathcal{L}_{1}^{*}(t)=\mathcal{L}_{M}^{c}(t)$, where

$$\begin{eqnarray}\begin{array}{rcl} {{\mathcal{L}}_{M}}(t){{A}^{v}} & = & {\rm i}{{(\left[ X(t),A(t) \right])}^{v}}; \\ \mathcal{L}_{M}^{c}(t){{A}^{v}} & = & -\frac{1}{2}{{(\{X(t),P(t)\})}^{v}}. \\ \end{array}\end{eqnarray} \tag{ 24 }$$

As before, the quantity A^v in the lhs is a vector with size ${{2}^{2N}}$ flattened from matrix A with size ${{2}^{N}}\times {{2}^{N}}$, where N is the OQS dimension. Finally, in this expression we have also settled

$$\begin{eqnarray}\begin{array}{rcl} \alpha _{00}^{0*}(t,s) & = & \frac{1}{2\hbar }\mathop{\sum }\limits_{\lambda }g_{\lambda }^{2}{{{\rm e}}^{-{\rm i}{{\omega }_{\lambda }}(t-s)}}={{\alpha }_{R}}(t-s) \\ \alpha _{11}^{0*}(t,s) & = & \frac{{\rm i}}{2\hbar }\mathop{\sum }\limits_{\lambda }f_{\lambda }^{2}{{{\rm e}}^{-{\rm i}{{\omega }_{\lambda }}(t-s)}} \\ {} & = & -2{{\alpha }_{I}}(t-s)\theta (t-s). \\ \end{array}\end{eqnarray} \tag{ 25 }$$

Note also that $\alpha _{00}^{*0}(t,s)={{(\alpha _{00}^{0*}(t,s))}^{*}}$, and $\alpha _{11}^{*0}(t,s)={{(\alpha _{11}^{0*}(t,s))}^{*}}$, and that $\alpha _{jj}^{00}(t,s)=\alpha _{jj}^{**}(t,s)=0$ for any j.

The number of members of each class, $A_{\lambda }^{\alpha }$ (referred here as first order), $B_{\lambda ^{\prime} \lambda }^{\alpha \beta }$ (second order), $C_{\lambda ^{\prime\prime} \lambda ^{\prime} \lambda }^{\alpha \beta \gamma }$ (third order), etc, are $P\times {{2}^{k}}$, where k is the order number, and P is the number of complex noises, in our case P = 2 (see figure 1). Thus, the master equation (23) is written in terms of an infinite series of terms, wherein the first two terms are in the form of an integral over ${{\rho }_{s}}$ involving one correlation function (5), and the following eight terms are in the form of a triple integral over ${{\rho }_{s}}$ involving two correlation functions. The next terms will come in the form of five integrals involving three correlation functions, and so on and so forth. The derivation of the equation does not rely on any approximation, but the resulting equation unveils a very complex structure with an infinite hierarchy of terms. Hence, for practical uses, the hierarchy should be truncated, which would of course lead to a certain error and to an approximated result.

Zoom In Zoom Out Reset image size

For this reason, it is important to take into account the convergence properties of the series. If the system parameters are such that the relative weight of the terms involving a single integral is smaller than that of the terms involving a triple integral, then the case is not tractable because the equation cannot be truncated in a controlled way. In addition, it is worth noticing at this point that the above-described hierarchical structure does not appear because of the use of any expansion approximation, such as the weak coupling. Rather, it is a structure intrinsically attached to the OQS problem. Hence, for systems or parameter ranges where the hierarchy structure does not converge (i.e. successive terms are not smaller than the previous ones), one might conclude that these particular instances shall not be correctly dealt with master equations. Indeed, in these situations the dynamics of the open quantum system appear to alter the environment so dramatically that obtaining a closed equation within the Hilbert space of the system is not feasible in practice. In contrast, for these problems where the subsequent terms are less important, one can safely truncate the equation. This will be seen more clearly in the next section, particularly when analysing the case where the environment is characterized with an exponential correlation function.

Weak coupling and the Markov limit

For consistency, let us check that equation (23) coincides with the master equation up to the second order in the coupling parameter g derived from other approaches. Indeed, each correlation function $\alpha _{jj}^{\alpha \beta }$ will be of the order g², so that if we just keep second order terms, the equation takes the form

$$\begin{eqnarray}\begin{array}{rcl} \frac{{\rm d}\rho _{s}^{v}(t)}{{\rm d}t} & = & \mathop{\sum }\limits_{j=0,1}\int _{0}^{t}{\rm d}s\{\hat{\alpha }_{jj}^{0*}(t-s)\mathcal{L}_{j}^{0}(t)\mathcal{L}_{j}^{*}(s) \\ {} & {} & \left. +\int _{0}^{t}{\rm d}s\hat{\alpha }_{jj}^{*0}(t-s)\mathcal{L}_{j}^{*}(t)\mathcal{L}_{j}^{0}(s) \right\}\rho _{s}^{v}(s). \\ \end{array}\end{eqnarray} \tag{ 26 }$$

Replacing the values $\mathcal{L}_{0}^{0}(t)=\mathcal{L}_{0}^{*}(t)={{\mathcal{L}}_{M}}(t)$, $\mathcal{L}_{1}^{0}(t)=2{{\mathcal{L}}_{M}}(t)$, and $\mathcal{L}_{1}^{*}(t)=\mathcal{L}_{M}^{c}(t)$, we find

$$\begin{eqnarray}\begin{array}{rcl} \frac{{\rm d}\rho _{s}^{v}(t)}{{\rm d}t} & = & -\int _{0}^{t}{\rm d}s\left\{ \operatorname{Re}\left[ \hat{\alpha }_{00}^{0*}(t-s) \right]{{\mathcal{L}}_{M}}(t){{\mathcal{L}}_{M}}(s) \right. \\ {} & {} & +{\rm i}\hat{\alpha }_{11}^{0*}(t-s){{\mathcal{L}}_{M}}(t)\mathcal{L}_{M}^{c}(s) \\ {} & {} & +{\rm i}\hat{\alpha }_{11}^{*0}(t-s)\mathcal{L}_{M}^{c}(t){{\mathcal{L}}_{M}}(s)\}\rho _{s}^{v}(s), \\ \end{array}\end{eqnarray} \tag{ 27 }$$

with ${{\mathcal{L}}_{M}}(t)$ and $\mathcal{L}_{M}^{c}(t)$ defined in (24). Re-writing the resulting equation in its matrix form, it is identical to the time-convoluted master equation [5, 22] up to second order in g

$$\begin{eqnarray}\begin{array}{rcl} \frac{{\rm d}{{\rho }_{s}}(t)}{{\rm d}t} & = & -\int _{0}^{t}{\rm d}s{{\alpha }_{T}}(t-s)X(t)X(s){{\rho }_{s}}(s) \\ {} & {} & -\int _{0}^{t}{\rm d}s\alpha _{T}^{*}(t-s){{\rho }_{s}}(s)X(s)X(t)+\int _{0}^{t}{\rm d}\tau {{\alpha }_{T}}(t-s)X(s) \\ {} & {} & \times {{\rho }_{s}}(s)X(t)+\int _{0}^{t}{\rm d}s\alpha _{T}^{*}(t-s)X(t){{\rho }_{s}}(s)X(s)+\mathcal{O}({{g}^{4}}). \\ \end{array}\end{eqnarray} \tag{ 28 }$$

As noted above, the equation is convoluted because the reduced density matrix is within the integral, and not local in time as in the so-called convolution-less master equation [5]. However, taking into account that in the interaction picture with respect to the system ${{\rho }_{s}}(s)={{\rho }_{s}}(t)+\mathcal{O}({{g}^{2}})$, and that the ${{\rho }_{s}}(s)$ already appears in second order terms, one can write the latter equation as

$$\begin{eqnarray}&&\begin{array}{l} \frac{{\rm d}{{\rho }_{s}}(t)}{{\rm d}t}=-\int _{0}^{t}{\rm d}s{{\alpha }_{T}}(t-s)X(t)X(s){{\rho }_{s}}(t) \\ -\int _{0}^{t}{\rm d}s\alpha _{T}^{*}(t-s){{\rho }_{s}}(t)X(s)X(t) \\ +\int _{0}^{t}{\rm d}s{{\alpha }_{T}}(t-s)X(s){{\rho }_{s}}(t)X(t). \\ \end{array}\end{eqnarray} \tag{ 29 }$$

This is already a time local second order master equation, identical to the one derived following the TCL projection operator technique [5], and the second order stochastic Schrödinger equations [23, 24].

If one considers the Markov limit of the master equation (23), by considering $\alpha (t)=\Gamma \delta (t)$, the terms involving three integrals (and two correlation functions) vanish, while the terms involving a single integral give rise to the well-known Lindblad equation [25]

$$\begin{eqnarray}&&\frac{{\rm d}{{\rho }_{s}}(t)}{{\rm d}t}=-\Gamma \left( X(t)X(t){{\rho }_{s}}(t)+{{\rho }_{s}}(t)X(t)X(t)) \right.\left. +2X(t){{\rho }_{s}}(t)X(t) \right).\end{eqnarray} \tag{ 30 }$$

An alternative to deriving the dynamics of the reduced density operator of the system is to re-write (2) as

$$\begin{eqnarray}&&\frac{{\rm d}P_{t}^{v}}{{\rm d}t}={{\mathcal{L}}_{S}}P_{t}^{v}+\mathcal{E}P_{t}^{v},\end{eqnarray} \tag{ 31 }$$

where we have defined $\mathcal{E}={\rm i}{{L}^{-}}{{z}_{0t}}+{\rm i}{{L}^{-}}{{z}_{1t}}-{{L}^{+}}z_{1t}^{*}$, with the super-operators defined as ${{L}^{-}}{{\rho }^{v}}={{(\{\hat{x},\rho \})}^{v}}$, and ${{L}^{+}}{{\rho }^{v}}={{(\{\hat{x},\rho \})}^{v}}$. Importantly, z_0t is now a real Gaussian coloured noise, such that $\mathcal{M}[{{z}_{0t}}{{z}_{0s}}]={{\alpha }_{R}}(t-s)$. In addition, $\mathcal{M}[{{z}_{1t}}z_{1s}^{*}]=-{{\alpha }_{I}}(t-s)\theta (t-s)$, which we simplify ${{\alpha }_{I}}(t-s)$, while $\mathcal{M}[z_{1t}^{*}{{z}_{1s}}]=-{{\alpha }_{I}}(t-s)\theta (s-t)$, here denoted as $\alpha _{I}^{*}(t-s)$. Notice that ${{\alpha }_{I}}(0)=\alpha _{I}^{*}(0)$. In order to compute the reduced density matrix, we consider the noise average of (31) as

$$\begin{eqnarray}&&\frac{{\rm d}\rho _{s}^{v}}{{\rm d}t}={{\mathcal{L}}_{S}}\rho _{s}^{v}+\rho 1,\end{eqnarray} \tag{ 32 }$$

where we have defined ${{\rho }_{1}}=\mathcal{M}[{{P}_{1}}]=\mathcal{M}[\mathcal{E}P_{t}^{v}]$, or more in detail

$$\begin{eqnarray}\begin{array}{rcl} {{\rho }_{1}} & = & {\rm i}{{L}^{-}}\int _{0}^{t}{\rm d}s{{\alpha }_{R}}(t-s)\mathcal{M}\left[ \frac{\delta P_{t}^{v}}{\delta {{z}_{0s}}} \right] \\ {} & {} & +{\rm i}{{L}^{-}}\int _{0}^{t}{\rm d}s\alpha _{I}^{*}(t-s)\mathcal{M}\left[ \frac{\delta P_{t}^{v}}{\delta z_{1s}^{*}} \right] \\ {} & {} & -{{L}^{+}}\int _{0}^{t}{\rm d}s{{\alpha }_{I}}(t-s)\mathcal{M}\left[ \frac{\delta P_{t}^{v}}{\delta {{z}_{1s}}} \right]. \\ \end{array}\end{eqnarray} \tag{ 33 }$$

Here, we have performed the integration by using parts of the Gaussian integrals to find that in general

$$\begin{eqnarray}\begin{array}{rcl} \mathcal{M}[{{z}_{0t}}P_{t}^{v}] & = & \int _{0}^{t}{\rm d}s{{\alpha }_{R}}(t-s)\mathcal{M}\left[ \frac{\delta P_{t}^{v}}{\delta {{z}_{0s}}} \right], \\ \mathcal{M}[{{z}_{1t}}P_{t}^{v}] & = & \int _{0}^{t}{\rm d}s\alpha _{I}^{*}(t-s)\mathcal{M}\left[ \frac{\delta P_{t}^{v}}{\delta z_{1s}^{*}} \right], \\ \mathcal{M}[z_{1t}^{*}P_{t}^{v}] & = & \int _{0}^{t}{\rm d}s{{\alpha }_{I}}(t-s)\mathcal{M}\left[ \frac{\delta P_{t}^{v}}{\delta {{z}_{1s}}} \right], \\ \end{array}\end{eqnarray} \tag{ 34 }$$

a result that is known as the Novikov theorem (see for instance [28]). As before, we proceed further by computing the time evolution of ${{\rho }_{1}}$, which can be written as a sum of the time evolution of three different types of terms, namely: $\mathcal{M}[{{z}_{0t}}P_{t}^{v}]$, $\mathcal{M}[{{z}_{1t}}P_{t}^{v}]$, and $\mathcal{M}[z_{1t}^{*}P_{t}^{v}]$. In the above expression, the integral can be extended to infinity because the stochastic projector is only non-zero in the interval $[0,t]$. Considering this, the time derivative of the first term of (34), $\mathcal{M}[{{z}_{0t}}P_{t}^{v}]$, can be written as

$$\begin{eqnarray}\begin{array}{rcl} \frac{{\rm d}\mathcal{M}[{{z}_{0t}}P_{t}^{v}]}{{\rm d}t} & = & \int _{0}^{\infty }{\rm d}s\frac{{\rm d}{{\alpha }_{R}}(t-s)}{{\rm d}t}\mathcal{M}\left[ \frac{\delta P_{t}^{v}}{\delta {{z}_{0s}}} \right] \\ {} & {} & +\int _{0}^{\infty }{\rm d}s{{\alpha }_{R}}(t-s)\mathcal{M}\left[ \frac{\delta }{\partial {{z}_{0s}}}({{\mathcal{L}}_{S}}P_{t}^{v}+\mathcal{E}P_{t}^{v}) \right], \\ \end{array}\end{eqnarray} \tag{ 35 }$$

where we have considered equation (31) to describe $\mathcal{M}\left[ (\delta /\delta {{z}_{0s}})({\rm d}P_{t}^{v}/{\rm d}t) \right]$. The last term in (35) can be rewritten as

$$\begin{eqnarray}&&\qquad \quad \mathcal{M}\left[ \frac{\delta }{\delta {{z}_{0s}}}({{\mathcal{L}}_{S}}P_{t}^{v}+\mathcal{E}P_{t}^{v}) \right]=\delta (t-s){{L}^{-}}\mathcal{M}[P_{t}^{v}]+\mathcal{M}\left[ \left( {{\mathcal{L}}_{S}}+\mathcal{E} \right)\frac{\delta P_{t}^{v}}{\delta {{z}_{0s}}} \right].\qquad \quad \end{eqnarray} \tag{ 36 }$$

From this, we can express (35) as

$$\begin{eqnarray}\begin{array}{rcl} \frac{{\rm d}\mathcal{M}[{{z}_{0t}}P_{t}^{v}]}{{\rm d}t} & = & \int _{0}^{\infty }{\rm d}s\frac{{\rm d}{{\alpha }_{R}}(t-s)}{{\rm d}t}\mathcal{M}\left[ \frac{\delta P_{t}^{v}}{\delta {{z}_{0s}}} \right]+{{\alpha }_{R}}(0){{L}^{-}}\rho _{s}^{v} \\ {} & {} & +\int _{0}^{\infty }{\rm d}s{{\alpha }_{R}}(t-s)\mathcal{M}\left[ \left( {{\mathcal{L}}_{S}}+\mathcal{E}) \right)\frac{\delta P_{t}^{v}}{\delta {{z}_{0s}}} \right], \\ \end{array}\end{eqnarray} \tag{ 37 }$$

where we have taken into account that $(\delta {{z}_{0t}}/\delta {{z}_{0s}})=\delta (t-s)$, and also that $\rho _{s}^{v}=\mathcal{M}[P_{t}^{v}]$. The time derivative of the two other terms, $\mathcal{M}[{{z}_{1t}}P_{t}^{v}]$ and $\mathcal{M}[z_{1t}^{*}P_{t}^{v}]$, contributing to ${{\rho }_{1}}$ can also be computed in a similar way, so that we can express

$$\begin{eqnarray}&&\frac{{\rm d}{{\rho }_{1}}}{{\rm d}t}={{\mathcal{L}}_{S}}{{\rho }_{1}}+{{\hat{\rho }}_{1}}+{{V}_{0}}\rho _{s}^{v}+{{\rho }_{2}},\end{eqnarray} \tag{ 38 }$$

where we have defined ${{V}_{0}}=-{{\alpha }_{R}}(0){{L}^{-}}{{L}^{-}}-{\rm i}{{\alpha }_{I}}(0)({{L}^{-}}{{L}^{+}}+{{L}^{+}}{{L}^{-}})$, and

$$\begin{eqnarray}&&\begin{array}{l} {{{\hat{\rho }}}_{1}}={\rm i}{{L}^{-}}\int _{0}^{\infty }{\rm d}s\frac{{\rm d}{{\alpha }_{R}}(t-s)}{{\rm d}t}\mathcal{M}\left[ \frac{\delta P_{t}^{v}}{\delta {{z}_{0s}}} \right] \\ \qquad +\;{\rm i}{{L}^{-}}\int _{0}^{\infty }{\rm d}s\frac{{\rm d}\alpha _{I}^{*}(t-s)}{{\rm d}t}\mathcal{M}\left[ \frac{\delta P_{t}^{v}}{\delta z_{1s}^{*}} \right] \\ \qquad -\;{{L}^{+}}\int _{0}^{\infty }{\rm d}s\frac{{\rm d}{{\alpha }_{I}}(t-s)}{{\rm d}t}\mathcal{M}\left[ \frac{\delta P_{t}^{v}}{\delta {{z}_{1s}}} \right]. \\ \end{array}\end{eqnarray} \tag{ 39 }$$

Also, we have defined the operator

$$\begin{eqnarray}&&\begin{array}{l} {{\rho }_{2}}=\mathcal{M}\left[ \mathcal{E}\left\{ {\rm i}{{L}^{-}}\int _{0}^{\infty }{\rm d}s{{\alpha }_{R}}(t-s)\frac{\delta P_{t}^{v}}{\delta {{z}_{0s}}}+{\rm i}{{L}^{-}} \right. \right. \\ \qquad \int _{0}^{\infty }{\rm d}s\left. \left. \alpha _{I}^{*}(t-s)\frac{\delta P_{t}^{v}}{\delta z_{1s}^{*}}-{{L}^{+}}\int _{0}^{\infty }{\rm d}s{{\alpha }_{I}}(t-s)\frac{\delta P_{t}^{v}}{\delta {{z}_{1s}}} \right\} \right]. \\ \end{array}\end{eqnarray} \tag{ 40 }$$

To rewrite ${{\rho }_{2}}$ in a different way, we shall consider the fact that ${{L}^{-}}{{L}^{+}}\equiv {{L}^{+}}{{L}^{-}}$, and that following integration by parts we have terms of the form

$$\begin{eqnarray}&&\int _{0}^{t}{\rm d}s{{\alpha }_{R}}(t-s)\mathcal{M}[{{z}_{0t}}\frac{\delta P}{\delta {{z}_{0s}}}]=\int _{0}^{t}{\rm d}{{s}^{\prime} }\int _{0}^{t}{\rm d}s{{\alpha }_{R}}(t-{{s}^{\prime} }){{\alpha }_{R}}(t-s)\mathcal{M}[\frac{{{\delta }^{2}}P}{\delta {{z}_{0{{s}^{\prime} }}}\delta {{z}_{0s}}}].\end{eqnarray} \tag{ 41 }$$

Note that the terms in ${{\rho }_{2}}$ that depend on ${{\alpha }_{I}}(t-s)$ and $\alpha _{I}^{*}(t-s)$ can be written in a similar way. Considering this, and also the fact that in all these terms the integrals can be extended to infinity, we can re-write ${{\rho }_{2}}$ as

$$\begin{eqnarray}&&\begin{array}{l} {{\rho }_{2}}={\rm i}{{L}^{-}}\mathcal{M}[\int _{0}^{\infty }{\rm d}s{{\alpha }_{R}}(t-s)\int _{0}^{\infty }{\rm d}{{s}^{\prime }}\left\{ {\rm i}{{L}^{-}}{{\alpha }_{R}}(t-{{s}^{\prime }})\frac{{{\delta }^{2}}P_{t}^{v}}{\delta {{z}_{0{{s}^{\prime }}}}\delta {{z}_{0s}}}+{\rm i}{{L}^{-}}\alpha _{I}^{*}(t-{{s}^{\prime }})\frac{{{\delta }^{2}}P_{t}^{v}}{\delta z_{1{{s}^{\prime }}}^{*}\delta {{z}_{0s}}} \right. \\ \left. \left. -{{L}^{+}}{{\alpha }_{I}}(t-{{s}^{\prime }})\frac{{{\delta }^{2}}P_{t}^{v}}{\delta {{z}_{1{{s}^{\prime }}}}\delta {{z}_{0s}}} \right\} \right]+{\rm i}{{L}^{-}}\mathcal{M}[\int _{0}^{\infty }{\rm d}s\alpha _{I}^{*}(t-s)\int _{0}^{\infty }{\rm d}{{s}^{\prime }} \\ \times \left\{ {\rm i}{{L}^{-}}{{\alpha }_{R}}(t-{{s}^{\prime }})\frac{{{\delta }^{2}}P_{t}^{v}}{\delta {{z}_{0{{s}^{\prime }}}}\delta z_{1s}^{*}}+{\rm i}{{L}^{-}}\alpha _{I}^{*}(t-{{s}^{\prime }})\frac{{{\delta }^{2}}P_{t}^{v}}{\delta z_{1{{s}^{\prime }}}^{*}\delta z_{1s}^{*}} \right.\left. \left. -{{L}^{+}}{{\alpha }_{I}}(t-{{s}^{\prime }})\frac{{{\delta }^{2}}P_{t}^{v}}{\delta {{z}_{1{{s}^{\prime }}}}\delta z_{1s}^{*}} \right\} \right] \\ -{{L}^{+}}\mathcal{M}[\int _{0}^{\infty }{\rm d}s{{\alpha }_{I}}(t-s)\int _{0}^{\infty }{\rm d}{{s}^{\prime }}\left\{ {\rm i}{{L}^{-}}{{\alpha }_{R}}(t-{{s}^{\prime }})\frac{{{\delta }^{2}}P_{t}^{v}}{\delta {{z}_{0{{s}^{\prime }}}}\delta {{z}_{1s}}}+{\rm i}{{L}^{-}}\alpha _{I}^{*}(t-{{s}^{\prime }})\frac{{{\delta }^{2}}P_{t}^{v}}{\delta z_{1{{s}^{\prime }}}^{*}\delta {{z}_{1s}}} \right. \\ \left. \left. -{{L}^{+}}{{\alpha }_{I}}(t-{{s}^{\prime }})\frac{{{\delta }^{2}}P_{t}^{v}}{\delta {{z}_{1{{s}^{\prime }}}}\delta {{z}_{1s}}} \right\} \right]. \\ \end{array}\end{eqnarray} \tag{ 42 }$$

After some straightforward calculations, it is found that the evolution of ${{\rho }_{2}}$ can be written similarly to (38)

$$\begin{eqnarray}&&\frac{{\rm d}{{\rho }_{2}}}{{\rm d}t}={{\mathcal{L}}_{S}}{{\rho }_{2}}+{{\hat{\rho }}_{2}}+V{{\rho }_{1}}+{{\rho }_{3}},\end{eqnarray} \tag{ 43 }$$

where we have defined a super-operator $V={\rm i}{{L}^{-}}{{\alpha }_{R}}(0)+{\rm i}{{L}^{-}}{{\alpha }_{I}}(0)-{{L}^{+}}{{\alpha }_{I}}(0)$. In this expression ${{\hat{\rho }}_{2}}$ is defined as the time derivative of only the correlation functions in ${{\rho }_{2}}$, according to (42). For instance, the first three terms of ${{\hat{\rho }}_{2}}$ can be written as

$$\begin{eqnarray}&&\begin{array}{l} {\rm i}{{L}^{-}}\mathcal{M}[\int _{0}^{\infty }{\rm d}s\int _{0}^{\infty }{\rm d}{{s}^{\prime} }\{{\rm i}{{L}^{-}}\frac{{\rm d}}{{\rm d}t}\left( {{\alpha }_{R}}(t-s){{\alpha }_{R}}(t-{{s}^{\prime} }) \right) \\ \quad \times \frac{{{\delta }^{2}}P_{t}^{v}}{\delta {{z}_{0{{s}^{\prime} }}}\delta {{z}_{0s}}}+{\rm i}{{L}^{-}}\frac{{\rm d}}{{\rm d}t}\left( {{\alpha }_{R}}(t-s)\alpha _{I}^{*}(t-{{s}^{\prime} }) \right)\frac{{{\delta }^{2}}P_{t}^{v}}{\delta z_{1{{s}^{\prime} }}^{*}\delta {{z}_{0s}}} \\ \quad \left. \left. -{{L}^{+}}\frac{{\rm d}}{{\rm d}t}\left( {{\alpha }_{R}}(t-s){{\alpha }_{I}}(t-{{s}^{\prime} }) \right)\frac{{{\delta }^{2}}P_{t}^{v}}{\delta {{z}_{1{{s}^{\prime} }}}\delta {{z}_{0s}}} \right\} \right]. \\ \end{array}\end{eqnarray} \tag{ 44 }$$

Even though the hierarchy appears to become more and more complex, at any order $n\gt 1$, the terms can always be conveniently recast in a similar way, that leads to a general equation of the form

$$\begin{eqnarray}&&\frac{{\rm d}{{\rho }_{n}}}{{\rm d}t}={{\mathcal{L}}_{S}}{{\rho }_{n}}+{{\hat{\rho }}_{n}}+V{{\rho }_{n-1}}+{{\rho }_{n+1}}.\end{eqnarray} \tag{ 45 }$$

The difficulty of solving this hierarchy numerically relies on the complexity of the terms ${{\hat{\rho }}_{n}}$, that in principle cannot be written in terms of the hierarchy terms ${{\rho }_{n}}$. One possibility to tackle this problem is to compute these terms with a parallel hierarchy, but this second hierarchy turns out to be coupled to yet another one and so on, and the resulting structure becomes too cumbersome to be dealt with reasonably. A way out of this problem is to consider a simple correlation function of the form $\alpha (t)=({{g}_{R}}+{\rm i}{{g}_{I}}){{{\rm e}}^{-\gamma t}}$. Then, ${{\hat{\rho }}_{n}}=-n\gamma {{\rho }_{n}}$, and the general term of the hierarchy (45) can be written as

$$\begin{eqnarray}&&\frac{{\rm d}{{\rho }_{n}}}{{\rm d}t}={{\mathcal{L}}_{S}}{{\rho }_{n}}-n\gamma {{\rho }_{n}}+V{{\rho }_{n-1}}+{{\rho }_{n+1}}.\end{eqnarray} \tag{ 46 }$$

For $n\gt 1$, while for n = 1 we have $({\rm d}{{\rho }_{1}}/{\rm d}t)={{\mathcal{L}}_{S}}{{\rho }_{1}}-\gamma {{\rho }_{1}}+{{V}_{0}}\rho _{s}^{v}+{{\rho }_{2}}$. Interestingly, the latter expression is very similar to the hierarchy obtained by Tanimura in [17, 26] (see also [18] for a review) from a path integral representation. As in the Tanimura derivation, each matrix in the hierarchy (46) depends on the previous and the next matrix. However, the present result differs from the Tanimura one; however, in our case each term ${{\rho }_{n}}$ is coupled to the term ${{\rho }_{n+1}}$ without any preceding super-operator, while in the Tanimura case the link with ${{\rho }_{n+1}}$ appears multiplied by a super-operator. Interestingly, Tanimura shows in [17] that his expansion is non-perturbative, i.e. its different terms do not correspond to the terms found in a perturbative expansion.

In addition, following Tanimura, the hierarchy can be truncated with a terminator of the form ${{\rho }_{N}}(t)\approx {\rm exp} [{{\mathcal{L}}_{S}}t]{{\rho }_{N}}(0)$ for $N\gamma$ sufficiently large [18]. A similar result can be obtained in the present case. To show this, let us consider the Nth order term of the hierarchy, which can formally be integrated as

$$\begin{eqnarray*}&&\rho _{N}^{I}(t)=\rho _{N}^{I}(0)+\int _{0}^{t}{\rm d}s{{{\rm e}}^{-N\gamma (t-s)}}({{\rho }_{N+1}}(s)+V{{\rho }_{N-1}}(s)).\end{eqnarray*}$$

Here we have considered interaction picture, such that for any n, $\rho _{n}^{I}(t)={\rm exp} [{{\mathcal{L}}_{S}}t]{{\rho }_{n}}(t)$. For $N\gamma \gg 1$, we can approximate ${\rm exp} [-N\gamma (t-s)]\approx (1/(N\gamma ))\delta (t-s)$, which replaced in the above expression leads to $\rho _{N}^{I}(t)\approx \rho _{N}^{I}(0)$ for large $N\gamma$. In the non-interaction picture, this can be written as ${{\rho }_{N}}(t)\approx {\rm exp} [{{\mathcal{L}}_{S}}t]{{\rho }_{N}}(0)$, which brings the same result as in the Tanimura case.

As already discussed in the previous section, the above argumentation shows that only in certain situations, i.e. when $N\gamma \gg 1$, the hierarchy presented in this section can be truncated, and the corresponding reduced density operator obtained.

Numerical example

In the following, we consider the evolution of a two level system with Hamiltonian ${{H}_{S}}=({{\omega }_{s}}/2){{\sigma }_{z}}$, and coupling to an environment with a coupling operator $\hat{x}={{\sigma }_{x}}.$ The environment is assumed to be characterized with a correlation function of the form $\alpha (t)={{{\rm e}}^{-\gamma t}}$. With this example, it will be shown: first, that the derived hierarchy (46) corresponds to the result obtained with the SLN equation (2); second, that for non-Markovian situations ($g\gamma \gt 1$), several members of the hierarchy are necessary to correctly reproduce the non-exponential decay of the system.

Figure 2 displays a comparison between the result obtained from a reduced density matrix using the hierarchy (solid black curve), and the one obtained with the SLN equation (2), in this last case when considering a different number of trajectories M in the sum (11). Indeed, when M = 1000 the curve given by the SLN is almost indistinguishable from the one obtained with the hierarchy, at least up to times of the order of 15. For longer times the statistics of the SLN equation begin to be more demanding, and more trajectories would be needed to get to the correct averaged result.

Zoom In Zoom Out Reset image size

In figure 3 we show the evolution of $\langle {{\sigma }_{z}}(t)\rangle$ obtained with the hierarchical equations (46) when considering different hierarchical levels $N=1,2,3$ and 4. It can be extracted from the figure that in a strongly non-Markovian situation (upper panel, with $g\gamma \gt 1$), there is a large difference between the results predicted by the hierarchy truncated at N = 1, and higher level truncations. While for N = 1 an exponential decay is observed, higher hierarchy levels unveil strong departures from this exponential form, which are due to the non-Markovian effects. The lower panel reproduces a situation in which, despite of the presence of a relative large coupling, the correlation time of the environment $1/\gamma$ is small as compared to the system evolution time scales $\sim 1/g$, so that the curves for all hierarchy levels practically show an exponential decay. For a decaying rate $\gamma \lt 1$, i.e. longer environment correlation times than the ones considered in the plots, the hierarchy performance is less optimal as there are situations in which it does not converge; at least at the relatively small N here considered.

Zoom In Zoom Out Reset image size