format: text

text:
# A perturbative approach to non-Markovian stochastic Schrödinger equations

## I Introduction

A common problem in physics is to model open quantum systems. They consists of a small system immersed in a bath (environment). Due to the large Hilbert space of the bath it is convenient to describe the system by its reduced state. The reduced state is defined as

$$\rho _{\mathrm{red}}(t)=\mathrm{Tr}_{\mathrm{bath}}[|\mathrm{\Psi }(t)\mathrm{\Psi }(t)|].$$
(1)
where $`|\mathrm{\Psi }(t)`$ is the combined system state, found from the Schrödinger equation for the open quantum system.

It has been shown Nak58 ; Zwa60  by a projection-operator method that we can write a general master equation for the reduced state as

$$\dot{\rho }_{\mathrm{red}}(t)=\frac{i}{\mathrm{}}[\widehat{H}(t),\rho _{\mathrm{red}}(t)]+_0^t𝒦(t,s)[\widehat{L}]\rho _{\mathrm{red}}(s)𝑑s$$
(2)
where $`𝒦(t,s)[\widehat{L}]`$ is the ‘memory time’ superoperator. It (operators on) the system operator $`\widehat{L}`$ and represents how the bath affects the system. The problem with this equation is that in general $`𝒦(t,s)[\widehat{L}]`$ can not be explicitly evaluated.

The most notable approximation used is the Born-Markov one. This arises when the environmental influences on the system are instantaneous. Mathematical consistency requires that this results in a Lindblad master equation, of the form Lin76

$$\dot{\rho }_{\mathrm{red}}(t)=\frac{i}{\mathrm{}}[\widehat{H}(t),\rho _{\mathrm{red}}(t)]+\gamma 𝒟[\widehat{L}]\rho _{\mathrm{red}}(t),$$
(3)
where $`𝒟[\widehat{L}]`$ is the superoperator that represent the damping of the system into the bath. It has the form

$$𝒟[\widehat{L}]\rho _{\mathrm{red}}=\widehat{L}\rho _{\mathrm{red}}\widehat{L}^{}\frac{1}{2}\widehat{L}^{}\widehat{L}\rho _{\mathrm{red}}\frac{1}{2}\rho _{\mathrm{red}}\widehat{L}^{}\widehat{L}.$$
(4)
This equation can be solved deterministically GarZol00  or by the stochastic Schrödinger approach GarZol00 ; DalCasMol92 ; GarParZol92 ; Car93 .

For the non-Markovian situation there have been many attempts at finding solutions to Eq. (2). However, some have the problem that it is hard to ensure the positivity requirement on $`\rho _{\mathrm{red}}(t)`$ Red65 . A method that does ensure the positivity requirement on the reduced state is the non-Markovian stochastic Schrödinger equation (SSE) approach Dio96 ; DioStr97 ; DioGisStr98 ; StrDioGis99 ; Str01 ; Cre00 ; Bud00 ; GamWis02 . A non-Markovian SSE generates stochastic pure states $`|\psi _z(t)`$ that should satisfy

$$\rho _{\mathrm{red}}(t)=E[|\psi _z(t)\psi _z(t)|],$$
(5)
where $`z(t)`$ is some noise function which is non-white and $`E`$ denotes the ensemble average over $`z(t)`$. To solve these non-Markovian SSE one has to take into account the past behavior of the system and bath, giving rise to a functional derivative in the attempt to derive a SSE. This presents a problem as for most systems an exact solution to the functional derivative does not exists. Thus at present an exact non-Markovian SSE only exists for simple systems, which can be solve exactly via other methods, like the undriven two level atom (TLA). For this and more examples see Ref. DioGisStr98 ; GamWis02 .

Recently Yu, Diósi, Gisin and Strunz (YDGS) have developed explicitly a ‘post-Markovian’ perturbation method to first order that allows solutions for systems that are close to the Markovian limit YuDioGisStr99 ; YuDioGisStr00 . In this paper we present a perturbation method that can be carried to arbitrary order and so is not limited to the post Markovian regime. However we must place a requirement on the form of the memory functions. This requirement is that the memory function must take the form

$$\alpha (ts)=\underset{j=1}{\overset{J}{}}|G_j|^2e^{\kappa _j|ts|/2i(\omega _j\mathrm{\Omega })(ts)},$$
(6)
for some finite (and, in practice, relatively small) $`J`$. It should be noted also that we have not proven convergence of our perturbation theory and this theory is only valid for a zero-temperature bath.

The format of this paper is as follows. In Sec. II we present a general outline of the theory of non-Markovian SSE. This is basically a summary of the results of Refs. Dio96 ; DioStr97 ; DioGisStr98 ; StrDioGis99 ; GamWis02 . In Sec. III our perturbation method is presented. In Secs. IV we outline Imamo$`\overline{\mathrm{g}}`$lu enlarged system method Ima94 ; SteIma96 . In Sec. V we apply our perturbation method to a simple system, a driven TLA and compare our results for $`\rho _{\mathrm{red}}(t)`$ with the enlarged system methods. In Sec. VI we investigate YDGS post-Markovian perturbation method YuDioGisStr99 ; YuDioGisStr00 . Finally we conclude with a discussion of the potential applications of our results in Sec. VII.

## II Non-Markovian Stochastic Schrödinger Equations

In this section we will present an outline of the theory we presented in GamWis02 , which is an extension of Diósi, Gisin and Strunz (DGS) diffusive Non-Markovian SSEs Dio96 ; DioStr97 ; DioGisStr98 ; StrDioGis99  which allows for real-valued noise $`z(t)`$.

### II.1 Underlying Dynamics

The non-Markovian SSEs developed in references GamWis02 ; Dio96 ; DioStr97 ; DioGisStr98 ; StrDioGis99  are valid when the dynamics of the open quantum system can be described by the total Hamiltonian

$$\widehat{H}_{\mathrm{tot}}=\widehat{H}_{\mathrm{sys}}\widehat{1}+\widehat{1}\widehat{H}_{\mathrm{bath}}+\widehat{V}.$$
(7)
The system Hamiltonian is $`\widehat{H}_{\mathrm{sys}}=\widehat{H}_\mathrm{\Omega }+\widehat{H}`$. The bath is modelled by a collection of harmonic oscillators, so the Hamiltonian for the bath is

$$\widehat{H}_{\mathrm{bath}}=\mathrm{}\underset{k}{}\omega _k\widehat{a}_k^{}\widehat{a}_k,$$
(8)
where $`\widehat{a}_k`$ and $`\omega _k`$ are the lowering operator and angular frequency of the $`k^{\mathrm{th}}`$ mode respectively. This is the standard model for the electromagnetic field. The interaction Hamiltonian has the form

$$\widehat{V}=i\mathrm{}(\widehat{L}\widehat{b}^{}\widehat{L}^{}\widehat{b}),$$
(9)
where we have defined the bath lowering operators $`\widehat{b}`$ as $`\widehat{b}=_kg_k\widehat{a}_k`$. That is, the coupling amplitude of the $`k^{\mathrm{th}}`$ mode to the system is $`g_k`$.

For calculation purposes we define the non-Markovian SSE in an interaction picture. This allows us to move the fast dynamics placed on the state by the Hamiltonians $`\widehat{H}_\mathrm{\Omega }`$ and $`\widehat{H}_{\mathrm{bath}}`$ to the operators. The unitary evolution operator for this transformations is

$$U(t,0)=e^{\frac{i}{\mathrm{}}(\widehat{H}_\mathrm{\Omega }\widehat{1}+\widehat{1}\widehat{H}_{\mathrm{bath}})(t0)}.$$
(10)
Thus the combined state in the interaction picture is define as

$$|\mathrm{\Psi }(t)=U^{}(t,0)|\mathrm{\Psi }(t)_{\mathrm{Sch}},$$
(11)
and an arbitrary operator $`\widehat{A}`$ becomes

$$\widehat{A}_{\mathrm{int}}=U^{}(t,0)\widehat{A}U(t,0).$$
(12)
This allows us to write the Schrödinger equation as

$$d_t|\mathrm{\Psi }(t)=\frac{i}{\mathrm{}}[\widehat{H}_{\mathrm{int}}(t)+\widehat{V}_{\mathrm{int}}(t)]|\mathrm{\Psi }(t),$$
(13)
where the Hamiltonians are

$`\widehat{H}_{\mathrm{int}}(t)`$ $`=`$ $`U^{}(t,0)\widehat{H}U(t,0),`$ (14)
and

$`\widehat{V}_{\mathrm{int}}(t)`$ $`=`$ $`i\mathrm{}[\widehat{L}e^{i\mathrm{\Omega }t}\widehat{b}_{\mathrm{int}}^{}(t)\widehat{L}^{}e^{i\mathrm{\Omega }t}\widehat{b}_{\mathrm{int}}(t)],`$ (15)
where

$`\widehat{b}_{\mathrm{int}}(t)`$ $`=`$ $`{\displaystyle \underset{k}{}}g_k\widehat{a}_ke^{i\omega _kt}.`$ (16)
Here we have finally restricted the form of $`\widehat{H}_\mathrm{\Omega }`$ to be such that $`\widehat{L}`$ in the interaction picture simply rotates in the complex plane at frequency $`\mathrm{\Omega }`$. That is $`\widehat{L}_{\mathrm{int}}(t)=\widehat{L}e^{i\mathrm{\Omega }t}`$.

### II.2 Non-Markovian SSE-Defined

A non-Markovian SSE is a stochastic differential equation for the system state vector $`|\psi _z(t)`$ containing some non-white noise $`z(t)`$. It has the property that when $`|\psi _z(t)\psi _z(t)|`$ is averaged over all possible $`z(t)`$ one obtains $`\rho _{\mathrm{red}}`$(t). It should be noted that for a single $`\rho _{\mathrm{red}}(t)`$, $`z(t)`$ can take many different functional forms, and we label these different forms as stochastic unravelings GamWis02 .

In Ref. GamWis02  we showed that non-Markovian SSEs can be derived from quantum measurement theory (QMT), where the different unravelings correspond to different measurements on the bath. The two unravelings we considered were the ‘coherent’ or ‘DGS’ Dio96 ; DioStr97 ; DioGisStr98 ; StrDioGis99  unraveling and the ‘quadrature’ unraveling. A special case of our quadrature unraveling was published in Ref. BasGhi02 .

As in the Markov limit we can define (at least) two non-Markovian SSEs, for each unraveling: one for $`z(t)`$ chosen from an ostensible distribution (a guessed distribution) and the other for its actual distribution. The former gives a non-Markovian SSE linear in the unnormalised state $`|\stackrel{~}{\psi }_z(t)`$, while the latter gives a non-Markovian SSE non-linear in the normalized state $`|\psi _z(t)`$. In Ref. GamWis02  we came to the conclusion that the solution of the actual non-Markovian SSE at time $`t`$ gives the state the system will be in if a measurement of the bath is performed at that time. Unlike in the Markov case, linking of the states through time to make a trajectory turns out to be a convenient fiction. However, it has been suggested that such trajectories can be given an interpretation within a non-standard QMT Lou01 ; BarLou02 .

#### II.2.1 Coherent Unravelling-Outlined

The first unravelling we consider is the ‘coherent’ unravelling. This unravelling arises when the bath is projected into a coherent state. We define the coherent state as

$$|\{a_k\}=\underset{k}{}\frac{1}{\sqrt{\pi }}e^{|a_k|^2/2}\underset{n_k}{}\frac{a_k^{n_k}}{\sqrt{n_k!}}|n_k,$$
(17)
so that $`\widehat{1}=_k|a_ka_k|d^2a_k`$.

In a measurement we can define an operator for the measurement process, the noise operator. For this measurement it must have the coherent basis as its eigenstate, so the noise operator is

$$\widehat{z}(t)=\widehat{b}_{\mathrm{int}}(t)e^{i\mathrm{\Omega }t}=\underset{k}{}g_k\widehat{a}_ke^{i\mathrm{\Omega }_kt},$$
(18)
where $`\mathrm{\Omega }_k=\omega _k\mathrm{\Omega }`$. The noise function (eigenvalue of the noise operator) is

$$z(t)=\underset{k}{}g_ka_ke^{i\mathrm{\Omega }_kt}.$$
(19)
where $`a_k`$ are the results of the projection in the coherent basis.

If we assume an ostensible distribution for $`a_k`$ as being the overlap of the coherent state with vacuum state, that is, it has the form

$$\mathrm{\Lambda }(\{a_k\})=\{0_k\}|\{a_k\}\{a_k\}|\{0_k\}=\pi ^Ke^{_k|a_k|^2},$$
(20)
where $`K=_k`$. With this ostensible distribution the noise function has the following correlations

$`\stackrel{~}{E}[z(t)z^{}(s)]`$ $`=`$ $`\alpha (ts),`$ (21a)
$`\stackrel{~}{E}[z(t)z(s)]`$ $`=`$ $`0.`$ (21b)
where the tilde above the $`E`$ refers to a average over the ostensible distribution. In Eq. (21a) we have defined $`\alpha (ts)`$, this function we label the memory function. On a microscopic level it has the form

$$\alpha (ts)=\underset{k}{}|g_k|^2e^{i\mathrm{\Omega }_k(ts)}.$$
(22)
Using the above ostensible distribution we can define a linear conditional system state as

$$|\stackrel{~}{\psi }_z(t)=\frac{\{a_k\}|\mathrm{\Psi }(t)}{\sqrt{\mathrm{\Lambda }(\{a_k\})}}.$$
(23)
Taking the time derivative and using Eq. (13) we get a linear differential equation for $`|\stackrel{~}{\psi }_z(t)`$ of the from

$`_t|\stackrel{~}{\psi }_z(t)`$ $`=`$ $`[{\displaystyle \frac{i}{\mathrm{}}}\widehat{H}_{\mathrm{int}}(t)+z^{}(t)\widehat{L}\widehat{L}^{}{\displaystyle _0^t}\alpha (ts)`$ (24)
$`\times {\displaystyle \frac{\delta }{\delta z^{}(s)}}ds]|\stackrel{~}{\psi }_z(t),`$
where $`\delta /\delta z^{}(s)`$ represents a functional derivative. For a derivation of this equation see Ref. DioStr97 ; GamWis02 . The functional derivative in this equation stops us from calling this equation a non-Markovian SSE, as it means that $`_t|\stackrel{~}{\psi }_z(t)`$ does not depend upon the state $`|\stackrel{~}{\psi }_z(t)`$ at all times for a single function $`z(t)`$, but rather also upon states for other noise functions. That is, we cannot stochastically choose $`z(t)`$ in order to generate a trajectory independent of other trajectories. Instead, all possible trajectories would have to be calculated in parallel, which in calculation terms amounts to solving the complete Schrödinger equation Eq. (13). However, as explained in Refs. GamWis02 ; DioGisStr98 ; StrDioGis99  if we can make the following ansatz

$$\frac{\delta }{\delta z^{}(s)}|\stackrel{~}{\psi }_z(t)={}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s)|\stackrel{~}{\psi }_z(t),$$
(25)
then we can write a linear non-Markovian SSE as

$$_t|\stackrel{~}{\psi }_z(t)=\left[\frac{i}{\mathrm{}}\widehat{H}_{\mathrm{int}}(t)+z^{}(t)\widehat{L}\widehat{L}^{}{}_{}{}^{(0)}\widehat{F}_{z}^{}(t)\right]|\stackrel{~}{\psi }_z(t),$$
(26)
where the operator functional $`{}_{}{}^{(0)}\widehat{F}_{z}^{}(t)`$ is defined as

$`{}_{}{}^{(0)}\widehat{F}_{z}^{}(t)={\displaystyle _0^t}\alpha (ts){}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s)𝑑s.`$ (27)
The significance of the superscripts $`(0)`$ proceeding these operators will become apparent in Sec. III.

To derive the actual (non-linear) non-Markovian SSE we need to condition the state on a noise function that is equivalent to the actual probability distribution,

$$P(\{a_k\},t)=\mathrm{\Psi }(t)|\{a_k\}\{a_k\}|\mathrm{\Psi }(t).$$
(28)
For most systems $`|\mathrm{\Psi }(t)`$ is unknown. Nevertheless we can use a Girsanov transformation GamWis02 ; DioGisStr98  to relate the actual noise function to the ostensible noise function. In this case,

$$z(t)=z_\mathrm{\Lambda }(t)+_0^t\alpha (ts)\widehat{L}_s𝑑s,$$
(29)
where $`z_\mathrm{\Lambda }(t)`$ is equivalent to the noise function used in the ostensible case, satisfying the correlations defined in Eqs. (21a) and (21b). With the correct $`z(t)`$ the actual non-Markovian SSE for the normalised state is GamWis02 ; DioGisStr98

$`d_t|\psi _z(t)`$ $`=`$ $`[{\displaystyle \frac{i}{\mathrm{}}}\widehat{H}_{\mathrm{int}}(t)(\widehat{L}^{}\widehat{L}^{}_t){}_{}{}^{(0)}\widehat{F}_{z}^{}(t)`$ (30)
$`+(\widehat{L}^{}\widehat{L}^{}_t){}_{}{}^{(0)}\widehat{F}_{z}^{}(t)_t+z^{}(t)`$
$`\times (\widehat{L}\widehat{L}_t)]|\psi _z(t).`$
The notation $`\widehat{L}_t`$ is short hand for $`\psi _z(t)|\widehat{L}|\psi _z(t)`$. From Eq. (26) and Eq. (30) if the operator functional $`{}_{}{}^{(0)}\widehat{F}_{z}^{}(t)`$ is known for all time and for each noise function $`z(t)`$ we can solve the coherent non-Markovian SSE.

#### II.2.2 Quadrature Unravelling-Outlined

To obtain a non-Markovian SSE with real noise, it is natural to consider a quadrature noise operator,

$$\widehat{z}(t)=\widehat{b}_{\mathrm{int}}(t)e^{i\omega _0t}e^{i\varphi }+\widehat{b}_{\mathrm{int}}^{}(t)e^{i\omega _0t}e^{i\varphi },$$
(31)
where $`\widehat{b}_{\mathrm{int}}(t)`$ is defined in equation (16) and $`\varphi `$ is some arbitrary phase. The phase $`\varphi `$ defines the measured quadrature: an $`x`$-quadrature measurement occurs when $`\varphi `$ is set to zero, and the conjugate measurement of the $`y`$-quadrature occurs when $`\varphi =\pi /2`$. Unless otherwise stated we will set $`\varphi `$ to zero.

The measurement basis for the bath measurement is $`|\{q_k\}`$ and must satisfy

$$\widehat{z}(t)|\{q_k\}=z(t)|\{q_k\}.$$
(32)
The problem with this noise function is that in general it is hard (maybe impossible) to work out a time-independent eigenstate in the interaction picture. However, we can find this eigenstate if we make the assumptions that for every mode $`k`$ there exists another mode, which we can label $`k`$, such that $`\mathrm{\Omega }_k=\mathrm{\Omega }_k`$ and $`g_k=g_k^{}`$. These assumptions simply mean that the modes coupled to the system come in symmetric pairs about the frequency $`\mathrm{\Omega }`$. Without loss of generality we can take the $`g_k`$’s to be real, absorbing any phases in the definitions of the bath operators. With all of these assumptions we can rewrite equation (31) as

$$\widehat{z}(t)=\underset{k>0}{}2g_k[\widehat{X}_k^+\mathrm{cos}(\mathrm{\Omega }_kt)+\widehat{Y}_k^{}\mathrm{sin}(\mathrm{\Omega }_kt)].$$
(33)
Here we have introduced the two-mode quadrature operators

$`\widehat{X}_k^\pm `$ $`=`$ $`(\widehat{x}_k\pm \widehat{x}_k)/\sqrt{2},`$ (34a)
$`\widehat{Y}_k^\pm `$ $`=`$ $`(\widehat{y}_k\pm \widehat{y}_k)/\sqrt{2},`$ (34b)
where $`\widehat{x}_k`$ and $`\widehat{y}_k`$ are the quadratures of $`\widehat{a}_k`$:

$$\widehat{a}_k=(\widehat{x}_k+i\widehat{y}_k)/\sqrt{2}.$$
(35)
The measurement basis that satisfies Eq. (32), in the x-quadrature representation is

$$|\{q_k\}=\underset{k>0}{}\frac{dx^{}}{\sqrt{2\pi }}|\frac{X_k^+x^{}}{\sqrt{2}}_k|\frac{X_k^++x^{}}{\sqrt{2}}_ke^{iY_k^{}x^{}}.$$
(36)
With this basis and the above noise operator the noise function for the quadrature measurement is

$$z(t)=\underset{k>0}{}2g_k[X_k^+\mathrm{cos}(\mathrm{\Omega }_kt)+Y_k^{}\mathrm{sin}(\mathrm{\Omega }_kt)],$$
(37)
which by definition is real.

Furthermore under the above assumptions the memory function $`\alpha (ts)`$ in Eq. (22) reduces to

$$\beta (ts)=2\underset{k>0}{}|g_k|^2\mathrm{cos}[\mathrm{\Omega }_k(ts)].$$
(38)
As in the coherent case we define the ostensible distribution as the overlap between the vacuum state and $`|\{q_k\}`$, that is

$$\mathrm{\Lambda }(\{X_k,Y_k\})=\pi ^{K/2}e^{_{k>0}(X_{k}^{+}{}_{}{}^{2}+Y_{k}^{}{}_{}{}^{2})}.$$
(39)
With this distribution the correlation for the real-valued noise function is

$$\stackrel{~}{E}[z(t)z(s)]=\beta (ts),$$
(40)
where the tilde, like before, means an average over the ostensible distribution. For this ostensible distribution the differential equation for $`|\stackrel{~}{\psi }_z(t)`$ is

$`_t|\stackrel{~}{\psi }_z(t)`$ $`=`$ $`[{\displaystyle \frac{i}{\mathrm{}}}\widehat{H}_{\mathrm{int}}(t)+z(t)\widehat{L}\widehat{L}_x{\displaystyle _0^t}\beta (ts)`$ (41)
$`\times {\displaystyle \frac{\delta }{\delta z(s)}}ds]|\stackrel{~}{\psi }_z(t),`$
where $`\widehat{L}_x=\widehat{L}+\widehat{L}^{}`$. Making the Ansatz,

$$\frac{\delta }{\delta z(s)}|\stackrel{~}{\psi }_z(t)={}_{}{}^{(0)}\widehat{q}_{z}^{}(t,s)|\stackrel{~}{\psi }_z(t),$$
(42)
the linear non-Markovian SSE is

$$_t|\stackrel{~}{\psi }_z(t)=\left[\frac{i}{\mathrm{}}\widehat{H}_{\mathrm{int}}(t)+z(t)\widehat{L}\widehat{L}_x{}_{}{}^{(0)}\widehat{Q}_{z}^{}(t)\right]|\stackrel{~}{\psi }_z(t),$$
(43)
where

$`{}_{}{}^{(0)}\widehat{Q}_{z}^{}(t)={\displaystyle _0^t}\beta (ts){}_{}{}^{(0)}\widehat{q}_{z}^{}(t,s)𝑑s.`$ (44)
To derive the actual non-Markovian SSE we need to calculate the correct noise function. The Girsanov transformation giving the actual real-valued $`z(t)`$ is GamWis02

$`z(t)`$ $`=`$ $`z_\mathrm{\Lambda }(t)+{\displaystyle _0^t}\widehat{L}_x_s\beta (ts)𝑑s,`$ (45)
where $`z_\mathrm{\Lambda }(t)`$ satisfies the correlations defined in Eq. (40). The actual non-Markovian SSE for the quadrature unravelling is

$`d_t|\psi _z(t)`$ $`=`$ $`[{\displaystyle \frac{i}{\mathrm{}}}\widehat{H}_{\mathrm{int}}(t)(\widehat{L}_x\widehat{L}_x_t){}_{}{}^{(0)}\widehat{Q}_{z}^{}(t)`$ (46)
$`+(\widehat{L}_x\widehat{L}_x_t){}_{}{}^{(0)}\widehat{Q}_{z}^{}(t)_t+z(t)`$
$`\times (\widehat{L}\widehat{L}_t)]|\psi _z(t).`$
Thus, if $`{}_{}{}^{(0)}\widehat{Q}_{z}^{}(t)`$ is known for $`z(t)`$ and all time then we can solve the quadrature non-Markovian SSE.

## III Perturbation Method

To solve the non-Markovian SSE, and hence find $`\rho _{\mathrm{red}}(t)`$, for the coherent or quadrature unravelling we have to work out the operator functionals $`{}_{}{}^{(0)}\widehat{F}_{z}^{}(t)`$ and $`{}_{}{}^{(0)}\widehat{Q}_{z}^{}(t)`$ respectively. This has been done exactly only for systems for which an analytical solution for $`\rho _{\mathrm{red}}(t)`$ may be found by other means DioGisStr98 ; StrDioGis99 ; Cre00  or for systems with a small number of bath modes GamWis02 . In this section we to propose our perturbation technique for working out these functionals when exact solutions are not possible.

### III.1 Perturbation Approach for the Coherent Unravelling

The perturbation that we are going to propose is only valid for memory functions of the form

$$\alpha (ts)=\underset{j=1}{\overset{J}{}}\alpha ^{(j)}(ts),$$
(47a)
where
$$\alpha ^{(j)}(ts)=|G_j|^2e^{\kappa _j|ts|/2}e^{i\mathrm{\Omega }_j(ts)}.$$
(47b)
In principle this is always a valid decomposition for the memory function as in the $`J\mathrm{}`$ and $`\kappa _j0`$ limit this memory function approaches the microscopic memory function displayed in Eq. (22). In Ref. SteIma96  the authors suggest that in practice most environments can be simulated with $`J`$ being quite small.

With this expansion for the memory function the functional $`{}_{}{}^{(0)}\widehat{F}_{z}^{}(t)`$ can be written as

$${}_{}{}^{(0)}\widehat{F}_{z}^{}(t)=\underset{j}{}{}_{}{}^{(0)}\widehat{F}_{z}^{(j)}(t).$$
(48)
where

$`{}_{}{}^{(0)}\widehat{F}_{z}^{(j)}(t)={\displaystyle _0^t}\alpha ^{(j)}(ts){}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s)𝑑s.`$ (49)
To calculate these operator functionals we set up a set of coupled nonlinear differential equations for $`{}_{}{}^{(0)}\widehat{F}_{z}^{(j)}(t)`$. Taking the time derivative of Eq. (49) we get

$`_t{}_{}{}^{(0)}\widehat{F}_{z}^{(j)}(t)`$ $`=`$ $`\alpha ^{(j)}(0){}_{}{}^{(0)}\widehat{f}_{z}^{}(t,t)+{\displaystyle _0^t}[_t\alpha ^{(j)}(ts)]`$ (50)
$`\times {}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s)ds+{\displaystyle _0^t}\alpha ^{(j)}(ts)`$
$`\times _t{}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s)ds.`$
The first term is easily evaluated using

$$^{(0)}\widehat{f}_z(t,t)=\widehat{L},$$
(51)
as derived in Appendix A. The second term is where our earlier decomposition of $`\alpha (ts)`$ is used. We chose $`\alpha ^{(j)}(ts)`$ such that $`_t\alpha ^{(j)}(ts)\alpha ^{(j)}(ts)`$. This results in the second term equaling

$$(\frac{\kappa _j}{2}+i\mathrm{\Omega }_j){}_{}{}^{(0)}\widehat{F}_{z}^{(j)}(t).$$
(52)
The third term involves the partial derivative $`_t[{}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s)]`$. To find this we use the fact that

$$_t\frac{\delta }{\delta z^{}(s)}|\stackrel{~}{\psi }_z(t)=\frac{\delta }{\delta z^{}(s)}_t|\stackrel{~}{\psi }_z(t),$$
(53)
which is called the consistency condition in DioGisStr98 . This consistency condition is only valid for $`ts`$ this is because at time $`t=s`$ the functional derivative is not well defined. Using Eq. (25) we can write the left-handed side (LHS) of the consistency condition as

$`_t{\displaystyle \frac{\delta }{\delta z^{}(s)}}|\stackrel{~}{\psi }_z(t)`$ $`=`$ $`[_t{}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s)]|\stackrel{~}{\psi }_z(t)`$ (54)
$`+{}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s)_t|\stackrel{~}{\psi }_z(t).`$
Substituting Eq. (26) in for $`_t|\stackrel{~}{\psi }_z(t)`$ gives

$`_t{\displaystyle \frac{\delta }{\delta z^{}(s)}}|\stackrel{~}{\psi }_z(t)`$ $`=`$ $`[_t{}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s){\displaystyle \frac{i}{\mathrm{}}}{}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s)\widehat{H}_{\mathrm{int}}(t)`$ (55)
$`+z^{}(t){}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s)\widehat{L}{}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s)`$
$`\times \widehat{L}^{}{}_{}{}^{(0)}\widehat{F}_{z}^{}(t)]|\stackrel{~}{\psi }_z(t).`$
Using Eqs. (26) and (25) the right-handed side (RHS) of the consistency condition gives

$`{\displaystyle \frac{\delta }{\delta z^{}(s)}}_t|\stackrel{~}{\psi }_z(t)`$ $`=`$ $`[{\displaystyle \frac{i}{\mathrm{}}}\widehat{H}_{\mathrm{int}}(t){}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s)+z^{}(t)\widehat{L}`$ (56)
$`\times {}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s)\widehat{L}^{}{}_{}{}^{(0)}\widehat{F}_{z}^{}(t){}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s)`$
$`\widehat{L}^{}{\displaystyle \frac{\delta }{\delta z^{}(s)}}{}_{}{}^{(0)}\widehat{F}_{z}^{}(t)]|\stackrel{~}{\psi }_z(t).`$
Equating the LHS with the RHS gives

$`_t{}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s)`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}[\widehat{H}_{\mathrm{int}}(t),{}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s)]+z^{}(t)[\widehat{L},`$ (57)
$`{}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s)][\widehat{L}^{}{}_{}{}^{(0)}\widehat{F}_{z}^{}(t),{}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s)]`$
$`\widehat{L}^{}{\displaystyle \frac{\delta }{\delta z^{}(s)}}{}_{}{}^{(0)}\widehat{F}_{z}^{}(t).`$
Substituting this equation with Eqs. (51) and (52) into Eq. (50) we get

$`_t{}_{}{}^{(0)}\widehat{F}_{z}^{(j)}(t)`$ $`=`$ $`|G_j|^2\widehat{L}({\displaystyle \frac{\kappa _j}{2}}+i\mathrm{\Omega }_j){}_{}{}^{(0)}\widehat{F}_{z}^{(j)}(t)+z^{}(t)`$
$`\times [\widehat{L},{}_{}{}^{(0)}\widehat{F}_{z}^{(j)}(t)]{\displaystyle \frac{i}{\mathrm{}}}[\widehat{H}_{\mathrm{int}}(t),{}_{}{}^{(0)}\widehat{F}_{z}^{(j)}(t)]`$
$`[\widehat{L}^{}{}_{}{}^{(0)}\widehat{F}_{z}^{}(t),{}_{}{}^{(0)}\widehat{F}_{z}^{(j)}(t)]\widehat{L}^{}`$
$`\times {\displaystyle \underset{k}{}}{}_{}{}^{(1)}\widehat{F}_{z}^{(j,k)}(t),`$
where $`{}_{}{}^{(1)}\widehat{F}_{z}^{(j,k)}(t)`$ is our first order functional. It has the form

$`{}_{}{}^{(1)}\widehat{F}_{z}^{(j,k)}(t)={\displaystyle _0^t}\alpha ^{(j)}(ts){}_{}{}^{(1)}\widehat{f}_{z}^{(k)}(t,s)𝑑s,`$ (59)
where we have used the following Ansatz

$$\frac{\delta }{\delta z^{}(s)}{}_{}{}^{(0)}\widehat{F}_{z}^{(k)}(t)={}_{}{}^{(1)}\widehat{f}_{z}^{(k)}(t,s).$$
(60)
If we knew the form of $`{}_{}{}^{(1)}\widehat{F}_{z}^{(j,k)}(t)`$ then Eq. (III.1) could be solved numerically.

To find the form of $`{}_{}{}^{(1)}\widehat{F}_{z}^{(j,k)}(t)`$ we can take the time derivative of Eq. (59). Doing this we get

$`_t{}_{}{}^{(1)}\widehat{F}_{z}^{(j,k)}(t)`$ $`=`$ $`\alpha ^{(j)}(0){}_{}{}^{(1)}\widehat{f}_{z}^{(k)}(t,t)+{\displaystyle _0^t}[_t\alpha ^{(j)}(ts)]`$ (61)
$`\times {}_{}{}^{(1)}\widehat{f}_{z}^{(k)}(t,s)ds+{\displaystyle _0^t}\alpha ^{(j)}(ts)`$
$`\times _t{}_{}{}^{(1)}\widehat{f}_{z}^{(k)}(t,s)ds.`$
The first term is easy to work out. From Eq. (III.1) it follows that

$${}_{}{}^{(1)}\widehat{f}_{z}^{(k)}(t,t)=[\widehat{L},{}_{}{}^{(0)}\widehat{F}_{z}^{(k)}(t)].$$
(62)
The second term as before also simply evaluates to

$$\left(\frac{\kappa _j}{2}+i\mathrm{\Omega }_j\right){}_{}{}^{(1)}\widehat{F}_{z}^{(j,k)}(t).$$
(63)
The third term is worked out via a new consistency condition,

$$_t\frac{\delta }{\delta z^{}(s)}{}_{}{}^{(0)}\widehat{F}_{z}^{(k)}(t)=\frac{\delta }{\delta z^{}(s)}_t{}_{}{}^{(0)}\widehat{F}_{z}^{(k)}(t).$$
(64)
Substituting Eqs. (60) and (III.1) into this consistency condition gives

$`_t{}_{}{}^{(1)}\widehat{f}_{z}^{(k)}(t,s)`$ $`=`$ $`({\displaystyle \frac{\kappa _k}{2}}+i\mathrm{\Omega }_k){}_{}{}^{(1)}\widehat{f}_{z}^{(k)}(t,s){\displaystyle \frac{i}{\mathrm{}}}[\widehat{H}_{\mathrm{int}}(t),{}_{}{}^{(1)}\widehat{f}_{z}^{(k)}(t,s)]+z^{}(t)[\widehat{L},{}_{}{}^{(1)}\widehat{f}_{z}^{(k)}(t,s)]`$ (65)
$`[\widehat{L}^{}{\displaystyle \underset{l}{}}{}_{}{}^{(1)}\widehat{f}_{z}^{(l)}(t,s),{}_{}{}^{(0)}\widehat{F}_{z}^{(k)}(t)][\widehat{L}^{}{\displaystyle \underset{l}{}}{}_{}{}^{(0)}\widehat{F}_{z}^{(l)}(t),{}_{}{}^{(1)}\widehat{f}_{z}^{(k)}(t,s)]`$
$`\widehat{L}^{}{\displaystyle \underset{l}{}}{\displaystyle \frac{\delta }{\delta z^{}(s)}}{}_{}{}^{(1)}\widehat{F}_{z}^{(k,l)}(t).`$
Substituting all these terms into Eq. (61) gives

$`_t{}_{}{}^{(1)}\widehat{F}_{z}^{(j,k)}(t)`$ $`=`$ $`|G_j|^2[\widehat{L},{}_{}{}^{(0)}\widehat{F}_{z}^{(k)}(t)]\left({\displaystyle \frac{\kappa _j}{2}}+i\mathrm{\Omega }_j\right){}_{}{}^{(1)}\widehat{F}_{z}^{(j,k)}(t)\left({\displaystyle \frac{\kappa _k}{2}}+i\mathrm{\Omega }_k\right){}_{}{}^{(1)}\widehat{F}_{z}^{(j,k)}(t)`$ (66)
$`{\displaystyle \frac{i}{\mathrm{}}}[\widehat{H}_{\mathrm{int}}(t),{}_{}{}^{(1)}\widehat{F}_{z}^{(j,k)}(t)]+z^{}(t)[\widehat{L},{}_{}{}^{(1)}\widehat{F}_{z}^{(j,k)}(t)][\widehat{L}^{}{\displaystyle \underset{l}{}}{}_{}{}^{(1)}\widehat{F}_{z}^{(j,l)}(t),{}_{}{}^{(0)}\widehat{F}_{z}^{(k)}(t)]`$
$`[\widehat{L}^{}{\displaystyle \underset{l}{}}{}_{}{}^{(0)}\widehat{F}_{z}^{(l)}(t),{}_{}{}^{(1)}\widehat{F}_{z}^{(j,k)}(t,s)]\widehat{L}^{}{\displaystyle \underset{l}{}}{}_{}{}^{(2)}\widehat{F}_{z}^{(j,k,l)}(t).`$
Where the last term is the second order functional, which equals

$${}_{}{}^{(2)}\widehat{F}_{z}^{(j,k,l)}(t)=_0^t\alpha ^{(j)}(ts)\frac{\delta }{\delta z^{}(s)}{}_{}{}^{(1)}\widehat{F}_{z}^{(k,l)}(t)𝑑s.$$
(67)
Here we see that we can develop a general way for setting up an $`n^{\mathrm{th}}`$ order differential equations. The $`n^{\mathrm{th}}`$ order functional is

$${}_{}{}^{(n)}\widehat{F}_{z}^{(j,k,\mathrm{},l)}(t)=_0^t\alpha ^{(j)}(ts){}_{}{}^{(n)}\widehat{f}_{z}^{(k,\mathrm{},l)}(t,s)𝑑s,$$
(68)
where we have used the Ansatz

$$\frac{\delta }{\delta z^{}(s)}{}_{}{}^{(n1)}\widehat{F}_{z}^{(k,\mathrm{},l)}(t)={}_{}{}^{(n)}\widehat{f}_{z}^{(k,\mathrm{},l)}(t,s).$$
(69)
The differential equation for the $`n^{\mathrm{th}}`$ order functional is

$`_t{}_{}{}^{(n)}\widehat{F}_{z}^{(j,k,\mathrm{},l)}(t)`$ $`=`$ $`\alpha ^{(j)}(0){}_{}{}^{(n)}\widehat{f}_{z}^{(k,\mathrm{},l)}(t,t)`$
$`+{\displaystyle _0^t}[_t\alpha ^{(j)}(ts)]{}_{}{}^{(n)}\widehat{f}_{z}^{(k,\mathrm{},l)}(t,s)𝑑s`$
$`+{\displaystyle _0^t}\alpha ^{(j)}(ts)_t{}_{}{}^{(n)}\widehat{f}_{z}^{(k,\mathrm{},l)}(t,s)ds.`$
The first term can always be calculated by the $`(n1)^{\mathrm{th}}`$ differential equation. The second term is always simple to calculate as $`_t\alpha ^{(j)}(ts)\alpha ^{(j)}(ts)`$ and the third term is always calculable by the $`(n1)^{\mathrm{th}}`$ order consistency condition

$$_t\frac{\delta }{\delta z^{}(s)}{}_{}{}^{(n1)}\widehat{F}_{z}^{(k,\mathrm{},l)}(t)=\frac{\delta }{\delta z^{}(s)}_t{}_{}{}^{(n1)}\widehat{F}_{z}^{(k,\mathrm{},l)}(t).$$
(71)
The $`n^{\mathrm{th}}`$ order perturbation method propose is to terminate this series by setting $`{}_{}{}^{(n)}\widehat{F}_{z}^{(j,k,\mathrm{},l)}(t)`$ equal to an arbitrary operator. The simplest scheme would be to set this operator to zero, but to keep the theory consistent with the Markov limit for all orders, we set $`{}_{}{}^{(n)}\widehat{F}_{z}^{(j,k,\mathrm{},l)}(t)`$ in the following manner. The zeroth order perturbation aries when we use the approximation

$${}_{}{}^{(0)}\widehat{F}_{z}^{(j)}(t)_0^t\alpha ^{(j)}(ts)𝑑s\frac{\delta }{\delta z^{}(t)}=_0^t\alpha ^{(j)}(ts)𝑑s\widehat{L}.$$
(72)
Note that the approximation here is the replacement of $`\delta /\delta z^{}(s)`$ by $`\delta /\delta z^{}(t)`$. The first order perturbation arises when we use the approximation

$`{}_{}{}^{(1)}\widehat{F}_{z}^{(j,k)}(t)`$ $``$ $`{\displaystyle _0^t}\alpha ^{(j)}(ts)𝑑s{\displaystyle \frac{\delta {}_{}{}^{(0)}\widehat{F}_{z}^{(k)}(t)}{\delta z^{}(t)}}`$ (73)
$`=`$ $`{\displaystyle _0^t}\alpha ^{(j)}(ts)𝑑s[\widehat{L},{}_{}{}^{(0)}\widehat{F}_{z}^{(k)}(t)]`$
and $`{}_{}{}^{(0)}\widehat{F}_{z}^{(j)}(t)`$ is calculated via Eq. (III.1). The $`n^{\mathrm{th}}`$ order perturbation arises when we use the approximation

$`{}_{}{}^{(n)}\widehat{F}_{z}^{(j,k,\mathrm{},l)}(t)`$ $``$ $`{\displaystyle _0^t}\alpha ^{(j)}(ts)𝑑s{\displaystyle \frac{\delta {}_{}{}^{(n1)}\widehat{F}_{z}^{(k,\mathrm{},l)}(t)}{\delta z^{}(t)}}`$
$`=`$ $`{\displaystyle _0^t}\alpha ^{(j)}(ts)𝑑s[\widehat{L},{}_{}{}^{(n1)}\widehat{F}_{z}^{(k,\mathrm{},l)}(t)]`$
and $`{}_{}{}^{(0)}\widehat{F}_{z}^{(j)}(t),\mathrm{},{}_{}{}^{(n1)}\widehat{F}_{z}^{(j,\mathrm{},k)}(t)`$ are calculated via Eqs. (III.1), (66) and (III.1). The physical motivations for choosing this type of expansion are;
a) For most system the memory function will decay and thus the most dominant term in the functional derivative will be the value as $`st`$.
b) Only $`{}_{}{}^{(0)}\widehat{F}_{z}^{(j)}(t)`$ affects the system directly, so the further removed the approximation the more accurate we expect the approximation to be.
c) In the Markovian limit, only the zero order term is needed.

To summarize this perturbation method, for environments which can be modelled by Eq. (47), it is possible to obtain a perturbative solution for the coherent non-Markovian SSE. From these SSEs it is possible to generate a perturbative solution for $`\rho _{\mathrm{red}}(t)`$, which by definition will always be positive. The number of coupled complex differential equations that are required for this technique is

$$d^2(J^n+J^{n1}+\mathrm{}+J)+d+J=d^2J\frac{J^n1}{J1}+d+J$$
(75)
where $`d`$ is the system dimension, $`J`$ is the number of exponentials required to simulate the memory function and $`n`$ is the order of the perturbation. The first term represents the number of equations needed to simulate the functional derivative. The next term $`d`$ is for the $`d`$ complex amplitudes of the system. The final term $`J`$ is for the stochastic equations needed to generate the noise function $`z(t)`$.

### III.2 Perturbation Approach for the Quadrature Unravelling

The perturbation method in the quadrature case is essentiality the same as the coherent case, but the memory function expressed in Eq. (47b) is too general. This is because the memory function for the quadrature unraveling must be consistent with the assumptions stated below Eq. (32). The most general memory function that satisfies these requirements is

$$\beta (ts)=\underset{j}{\overset{^{}}{}}\beta ^{(j,\mathrm{cos})}(ts),$$
(76a)
where
$$\beta ^{(j,\mathrm{cos})}(ts)=2|G_j|^2e^{\kappa _j|ts|/2}\mathrm{cos}[\mathrm{\Omega }_j(ts)].$$
(76b)
This presents a problem as $`_t\beta ^{(j,\mathrm{cos})}(ts)`$ is not proportional to $`\beta ^{(j,\mathrm{cos})}(ts)`$. To get around this we define a new function $`\beta ^{(j,\mathrm{sin})}(ts)`$ as

$$\beta ^{(j,\mathrm{sin})}(ts)=2|G_j|^2e^{\kappa _j|ts|/2}\mathrm{sin}(\mathrm{\Omega }_j(ts)),$$
(77)
and two functionals

$`{}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{cos})}(t)`$ $`=`$ $`{\displaystyle _0^t}\beta ^{(j,\mathrm{cos})}(ts)\widehat{q}_z(t,s)𝑑s,`$ (78a)
$`{}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{sin})}(t)`$ $`=`$ $`{\displaystyle _0^t}\beta ^{(j,\mathrm{sin})}(ts)\widehat{q}_z(t,s)𝑑s.`$ (78b)
The functional $`{}_{}{}^{(0)}\widehat{Q}_{z}^{}(t)`$ is the found by

$${}_{}{}^{(0)}\widehat{Q}_{z}^{}(t)=\underset{j}{}{}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{cos})}(t).$$
(79)
Taking the time derivative of Eqs. (78a) and (78b) we get

$`d_t{}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{cos})}(t)`$ $`=`$ $`\beta ^{(j,\mathrm{cos})}(t,t){}_{}{}^{(0)}\widehat{q}_{z}^{}(t,t)`$
$`+{\displaystyle _0^t}[_t\beta ^{(j,\mathrm{cos})}(ts)]{}_{}{}^{(0)}\widehat{q}_{z}^{}(t,s)𝑑s`$
$`+{\displaystyle _0^t}\beta ^{(j,\mathrm{cos})}(ts)_t{}_{}{}^{(0)}\widehat{q}_{z}^{}(t,s)ds,`$
$`d_t{}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{sin})}(t)`$ $`=`$ $`{\displaystyle _0^t}[_t\beta ^{(j,\mathrm{sin})}(ts)]{}_{}{}^{(0)}\widehat{q}_{z}^{}(t,s)𝑑s`$
$`+{\displaystyle _0^t}\beta ^{(j,\mathrm{sin})}(ts)_t{}_{}{}^{(0)}\widehat{q}_{z}^{}(t,s)ds.`$
As in the coherent case it can be shown that $`{}_{}{}^{(0)}\widehat{q}_{z}^{}(t,t)=\widehat{L}`$. The two terms involving the derivative of $`\beta ^{(j,\mathrm{cos})}(ts)`$ and $`\beta ^{(j,\mathrm{sin})}(ts)`$ by definition give

$`{\displaystyle _0^t}_t\beta ^{(j,\mathrm{cos})}(ts){}_{}{}^{(0)}\widehat{q}_{z}^{}(t,s)ds`$ $`=`$ $`{\displaystyle \frac{\kappa _j}{2}}{}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{cos})}(t)`$ (81a)
$`\mathrm{\Omega }_j{}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{sin})}(t)`$
$`{\displaystyle _0^t}_t\beta ^{(j,\mathrm{sin})}(ts){}_{}{}^{(0)}\widehat{q}_{z}^{}(t,s)ds`$ $`=`$ $`{\displaystyle \frac{\kappa _j}{2}}{}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{sin})}(t)`$ (81b)
$`+\mathrm{\Omega }_j{}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{cos})}(t).`$
The last two terms require finding $`_t{}_{}{}^{(0)}\widehat{q}_{z}^{}(t,s)`$. As in the coherent case this is found by the consistency condition

$$_t\frac{\delta }{\delta z(s)}|\stackrel{~}{\psi }_z(t)=\frac{\delta }{\delta z(s)}_t|\stackrel{~}{\psi }_z(t),$$
(82)
yielding

$`_t{}_{}{}^{(0)}\widehat{q}_{z}^{}(t,s)`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}[\widehat{H}_{\mathrm{int}}(t),{}_{}{}^{(0)}\widehat{q}_{z}^{}(t,s)]+z(t)[\widehat{L},`$ (83)
$`{}_{}{}^{(0)}\widehat{q}_{z}^{}(t,s)][\widehat{L}_x{}_{}{}^{(0)}\widehat{Q}_{z}^{}(t),{}_{}{}^{(0)}\widehat{q}_{z}^{}(t,s)]`$
$`\widehat{L}_x{\displaystyle \frac{\delta }{\delta z(s)}}{}_{}{}^{(0)}\widehat{Q}_{z}^{}(t).`$
Substituting these terms into Eq. (80) we get

$`d_t{}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{cos})}(t)`$ $`=`$ $`2|G_j|^2\widehat{L}{\displaystyle \frac{\kappa _j}{2}}{}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{cos})}(t)\mathrm{\Omega }_j{}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{sin})}(t){\displaystyle \frac{i}{\mathrm{}}}[\widehat{H}_{\mathrm{int}}(t),{}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{cos})}(t)]`$ (84a)
$`+z(t)[\widehat{L},{}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{cos})}(t)][\widehat{L}_x{}_{}{}^{(0)}\widehat{Q}_{z}^{}(t),{}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{cos})}(t)]`$
$`\widehat{L}_x{\displaystyle \underset{k}{}}{}_{}{}^{(1)}\widehat{Q}_{z}^{(j,k,\mathrm{cos},\mathrm{cos})}(t),`$
$`d_t{}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{sin})}(t)`$ $`=`$ $`{\displaystyle \frac{\kappa _j}{2}}{}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{sin})}(t)+\mathrm{\Omega }_j{}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{cos})}(t){\displaystyle \frac{i}{\mathrm{}}}[\widehat{H}_{\mathrm{int}}(t),{}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{sin})}(t)]+z(t)[\widehat{L},{}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{sin})}(t)]`$ (84b)
$`[\widehat{L}_x{}_{}{}^{(0)}\widehat{Q}_{z}^{}(t),^{(0)}{}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{sin})}(t)]\widehat{L}_x{\displaystyle \underset{k}{}}{}_{}{}^{(1)}\widehat{Q}_{z}^{(j,k,\mathrm{sin},\mathrm{cos})}(t),`$
where

$${}_{}{}^{(1)}\widehat{Q}_{z}^{(j,k,\mathrm{cos},\mathrm{cos})}(t)=_0^t\beta ^{(j,\mathrm{cos})}(t,s)\frac{\delta {}_{}{}^{(0)}\widehat{Q}_{z}^{(k,\mathrm{cos})}(t)}{\delta z(s)}𝑑s,$$
(85a)
$${}_{}{}^{(1)}\widehat{Q}_{z}^{(j,k,\mathrm{sin},\mathrm{cos})}(t)=_0^t\beta ^{(j,\mathrm{sin})}(t,s)\frac{\delta {}_{}{}^{(0)}\widehat{Q}_{z}^{(k,\mathrm{cos})}(t)}{\delta z(s)}𝑑s.$$
(85b)
The higher order functional differential equations are found in the same manner as in the coherent case, except the form of $`\beta (ts)`$ results in $`2^n`$ as many equations for order $`n`$.

The perturbation expansion is similar for this unravelling, the only difference being that we have $`2^n`$ operators to approximate. The $`0^{\mathrm{th}}`$ order approximation is to set the $`0^{\mathrm{th}}`$ order functionals to

$`{}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{cos})}(t)`$ $`=`$ $`{\displaystyle _0^t}\beta ^{(j,\mathrm{cos})}(t,s)𝑑s\widehat{L}`$ (86a)
$`{}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{sin})}(t)`$ $`=`$ $`{\displaystyle _0^t}\beta ^{(j,\mathrm{sin})}(t,s)𝑑s\widehat{L}.`$ (86b)
The first order approximation is to set the four first order functionals to

$`{}_{}{}^{(1)}\widehat{Q}_{z}^{(j,k,\mathrm{cos},\mathrm{cos})}(t)`$ $`=`$ $`{\displaystyle _0^t}\beta ^{(j,\mathrm{cos})}(t,s)𝑑s[\widehat{L},{}_{}{}^{(0)}\widehat{Q}_{z}^{(k,\mathrm{cos})}(t)],`$
$`{}_{}{}^{(1)}\widehat{Q}_{z}^{(j,k,\mathrm{sin},\mathrm{cos})}(t)`$ $`=`$ $`{\displaystyle _0^t}\beta ^{(j,\mathrm{sin})}(t,s)𝑑s[\widehat{L},{}_{}{}^{(0)}\widehat{Q}_{z}^{(k,\mathrm{cos})}(t)],`$
$`{}_{}{}^{(1)}\widehat{Q}_{z}^{(j,k,\mathrm{cos},\mathrm{sin})}(t)`$ $`=`$ $`{\displaystyle _0^t}\beta ^{(j,\mathrm{cos})}(t,s)𝑑s[\widehat{L},{}_{}{}^{(0)}\widehat{Q}_{z}^{(k,\mathrm{sin})}(t)],`$
$`{}_{}{}^{(1)}\widehat{Q}_{z}^{(j,k,\mathrm{sin},\mathrm{sin})}(t)`$ $`=`$ $`{\displaystyle _0^t}\beta ^{(j,\mathrm{sin})}(t,s)𝑑s[\widehat{L},{}_{}{}^{(0)}\widehat{Q}_{z}^{(k,\mathrm{sin})}(t)].`$
and we calculate the $`0^{\mathrm{th}}`$ order functionals via Eq. (84).

## IV Enlarged System Approach

To test the accuracy of our perturbation method we compare our results for the reduced state with the reduced state found via the enlarged system method of Imamo$`\overline{\mathrm{g}}`$lu Ima94 ; SteIma96 . An example of how this method is applied to a non-Markovian system can be found in Ref. CaoLonWeiCao01 .

For those who are not familiar with the enlarged system method, we provide a short proof that the reduced system dynamics are exactly reproduced by the enlarged system method provided that $`\alpha (ts)`$, called $`\mathrm{\Gamma }(\tau )`$ in Refs Ima94 ; SteIma96 , is of the form

$$\alpha (ts)=\underset{j}{}|G_j|^2e^{\kappa _j|ts|/2i\mathrm{\Omega }_j(ts)},$$
(88)
which is the same as Eq. (47).

The total Hamiltonian for the enlarged system is

$`\widehat{H}_{\mathrm{tot}}`$ $`=`$ $`\widehat{H}_{\mathrm{sys}}+\mathrm{}{\displaystyle \underset{j}{}}\omega _j\widehat{c}_j^{}\widehat{c}_j+\mathrm{}{\displaystyle \underset{j}{}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\omega `$
$`\times \omega \widehat{\nu }_j(\omega )^{}\widehat{\nu }_j(\omega )+i\mathrm{}{\displaystyle \underset{j}{}}[G_j^{}\widehat{L}\widehat{c}_j^{}G_j\widehat{L}^{}\widehat{c}_j]`$
$`+i\mathrm{}{\displaystyle \underset{j}{}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\omega \sqrt{{\displaystyle \frac{\kappa _j}{2\pi }}}[\widehat{\nu }_j^{}(\omega )\widehat{c}_j\widehat{\nu }_j(\omega )\widehat{c}_j^{}],`$
where $`\widehat{H}_{\mathrm{sys}}=\widehat{H}_\mathrm{\Omega }+\widehat{H}`$, $`\widehat{c}_j`$ is the annihilation operator for the $`j^{\mathrm{th}}`$ added oscillator and $`\widehat{\nu }_j(\omega )`$ is the Markovian bath operator with the correlation

$$[\widehat{\nu }_j(\omega ),\widehat{\nu }_k^{}(\omega )]=\delta _{j,k}\delta (\omega \omega ^{}).$$
(90)
If this is to be the same as Eq. (7), then the first two lines of Eq. (IV) must give $`\widehat{H}_{\mathrm{sys}}+\widehat{H}_{\mathrm{bath}}`$ and the final line $`\widehat{V}`$. Going to the same interaction picture as we did in Sec. II.1, that is with respect to the Hamiltonians $`\widehat{H}_\mathrm{\Omega }`$ and $`\widehat{H}_{\mathrm{bath}}`$, we get

$$\widehat{V}_{\mathrm{int}}(t)=i\mathrm{}\underset{j}{}[G_j^{}\widehat{L}e^{i\mathrm{\Omega }t}\widehat{c}_j(t)^{}G_j\widehat{L}^{}e^{i\mathrm{\Omega }t}\widehat{c}_j(t)].$$
(91)
Comparing with Eq. (15), for the enlarged system method to be correct we need $`\widehat{b}_{\mathrm{int}}(t)=_jG_j\widehat{c}_j(t)`$. To calculate $`\widehat{c}_j(t)`$ we use the fact that

$$d_t\widehat{c}_j(t)=i\omega _j\widehat{c}_j(t)\frac{\kappa _j}{2}\widehat{c}_j(t)\sqrt{\kappa _j}\widehat{\nu }_{\mathrm{in},j}(t)$$
(92)
where $`\widehat{\nu }_{\mathrm{in},j}(t)`$ is the input field which has a time commutator $`[\widehat{\nu }_{\mathrm{in},j}(t),\widehat{\nu }_{\mathrm{in},k}^{}(s)]=\delta _{j,k}\delta (ts)`$. For a derivation of equation Eq. (92) see Ref. GarCol85 . This can be integrated to give

$`\widehat{c}_j(t)`$ $`=`$ $`\sqrt{\kappa _j}{\displaystyle _0^t}e^{\kappa _j(ts)/2i\omega _j(ts)}\widehat{\nu }_{\mathrm{in},j}(s)𝑑s`$ (93)
$`+\widehat{c}_j(0)e^{\kappa _jt/2i\omega _jt}.`$
It not obvious that $`_jG_j\widehat{c}_j(t)`$ is the same as Eq. (16). However the time commutator for the bath operators is

$$[\widehat{b}_{\mathrm{int}}(t),\widehat{b}_{\mathrm{int}}^{}(s)]e^{i\mathrm{\Omega }(ts)}=\alpha (ts).$$
(94)
In terms of the enlarged system this means

$`{\displaystyle \underset{j,k}{}}G_jG_k^{}[\widehat{c}_j(t),\widehat{c}_k^{}(s)]e^{i\mathrm{\Omega }(ts)}`$
$`={\displaystyle \underset{j}{}}|G_j|^2e^{\kappa _j(t+s)/2i(\omega _j\mathrm{\Omega })(ts)}[1`$
$`+\kappa _j{\displaystyle _0^t}{\displaystyle _0^s}e^{+\kappa _j(t^{}+s^{})/2+i\omega _j(t^{}s^{})}\delta (t^{}s^{})dt^{}ds^{}]`$
$`={\displaystyle \underset{j}{}}|G_j|^2e^{\kappa _j|ts|/2i(\omega _j\mathrm{\Omega })(ts)}`$
$`=\alpha (ts),`$ (95)
provide $`\alpha (ts)`$ has the form depicted in Eq. (88). It should noted that this result is exact. It is not necessary to discard initial transients as in the derivation in Ref SteIma96 .

Since we have shown that the total Hamiltonian for the enlarged system is equivalent to the standard non-Markovian, then the total states $`|\mathrm{\Psi }_{Sch}(t)`$ must be the same. We can define a reduced state (in the Schrödinger picture) for the enlarged system as $`W_{\mathrm{Sch}}(t)`$ which has the Markovian master equation

$`d_tW_{\mathrm{Sch}}(t)`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}[\widehat{H}_\mathrm{\Omega }+\widehat{H}+\mathrm{}{\displaystyle \underset{j}{}}\omega _j\widehat{c}_j^{}\widehat{c}_j+i\mathrm{}{\displaystyle \underset{j}{}}`$ (96)
$`\times (G_j^{}\widehat{L}\widehat{c}_j^{}G_j\widehat{L}^{}\widehat{c}_j),W_{\mathrm{Sch}}(t)]`$
$`+{\displaystyle \underset{j}{}}\kappa _j𝒟[\widehat{c}_j]W_{\mathrm{Sch}}(t).`$
The reduced state for the system in the $`\mathrm{\Omega }`$-interaction picture is

$$\rho _{\mathrm{red}}(t)=e^{\frac{i}{\mathrm{}}\widehat{H}_\mathrm{\Omega }t}Tr_{\mathrm{enl}}[W_{\mathrm{Sch}}(t)]e^{\frac{i}{\mathrm{}}\widehat{H}_\mathrm{\Omega }t}=Tr_{\mathrm{enl}}[W_{\mathrm{red}}(t)],$$
(97)
where the trace is performed over the added oscillators and

$$W_{\mathrm{red}}(t)=e^{i_j\omega _j\widehat{c}_j^{}\widehat{c}_jt+\frac{i}{\mathrm{}}\widehat{H}_\mathrm{\Omega }t}W_{\mathrm{Sch}}(t)e^{i_j\omega _j\widehat{c}_j^{}\widehat{c}_jt\frac{i}{\mathrm{}}\widehat{H}_\mathrm{\Omega }t}.$$
(98)
This allows us to define a new master equation for the reduced state $`W_{\mathrm{red}}(t)`$ as

$`d_tW_{\mathrm{red}}(t)`$ $`=`$ $`[{\displaystyle \frac{i}{\mathrm{}}}\widehat{H}_{\mathrm{int}}+{\displaystyle \underset{j}{}}[G_j^{}\widehat{L}\widehat{c}_j^{}e^{i(\omega _j\mathrm{\Omega })t}`$
$`G_j\widehat{L}^{}\widehat{c}_je^{i(\omega _j\mathrm{\Omega })t}],W_{\mathrm{red}}(t)]`$
$`+{\displaystyle \underset{j}{}}\kappa _j𝒟[\widehat{c}_j]W_{\mathrm{red}}(t).`$
which can be solved by standard Markovian techniques, for example quantum trajectories DalCasMol92 ; GarParZol92 ; Car93 .

## V Numerical Example: The Driven Two Level Atom

In this section we apply our theory to a driven TLA with a simple non-Markovian memory function.

$$\alpha (ts)=\frac{\gamma \kappa }{4}e^{i(\omega _{\mathrm{env}}\mathrm{\Omega })(ts)}e^{\kappa |ts|/2},$$
(100)
where $`\omega _{\mathrm{env}}`$ is the central frequency of the environment, $`\kappa `$ represent the exponential decay of bath memory and $`\gamma `$ is the Markovian limit decay rate. That is, in the $`\kappa \mathrm{}`$ limit, $`\alpha (ts)=\gamma \delta (ts)`$, which is the Markovian limit of the memory function GamWis02 . We choose an interaction picture such the $`\mathrm{\Omega }=\omega _{\mathrm{env}}`$ so that this memory function is simplifies to

$$\alpha (ts)=\frac{\gamma \kappa }{4}e^{\kappa |ts|/2},$$
(101)
which is consistent with the quadrature unravelings assumptions. This results in $`\alpha (ts)=\beta (ts)`$. However before we apply our theory to the TLA let us revise the standard TLA model.

### V.1 The TLA

The TLA is one of the most simple quantum systems to envisage. It consists of two levels, an excited state $`|e`$ of energy $`\mathrm{}\omega _e`$ and a ground state $`|g`$ of energy $`\mathrm{}\omega _g`$. We define the difference in these energies as $`\mathrm{}\omega _0`$ and the zero point energy is taken to be the mid point energy $`\mathrm{}(\omega _e+\omega _g)/2=0`$. This allows us to define a system Hamiltonian as

$$\widehat{H}_{\mathrm{sys}}=\mathrm{}\frac{\omega _0}{2}\widehat{\sigma }_z$$
(102)
where $`\widehat{\sigma }_z=|ee||gg|`$ is one of the spin matrices for the TLA.

Since we are dealing with open quantum systems we consider the dynamics of the TLA immersed in the electromagnetic field (the bath). In the Schrödinger picture with the dipole and rotating wave approximation (RWA) approximation the interaction Hamiltonian is

$`\widehat{V}=i\mathrm{}{\displaystyle \underset{k}{}}(g_k^{}\widehat{\sigma }\widehat{a}_k^{}g_k\widehat{\sigma }^{}\widehat{a}_k),`$ (103)
where $`\widehat{\sigma }`$ is the lowering operator for the TLA. This is the same form as Eq. (9) with $`\widehat{L}=\widehat{\sigma }`$, so the above non-Markovian SSE theory is applicable to this system.

If we have a TLA driven by a classical electromagnetic field the system Hamiltonian for the TLA under the RWA approximation is

$$\widehat{H}_{\mathrm{sys}}=\mathrm{}\frac{\omega _0}{2}\widehat{\sigma }_z+\mathrm{}\frac{\chi }{2}[\widehat{\sigma }e^{i\omega _{\mathrm{dr}}t}+\widehat{\sigma }^{}e^{i\omega _{\mathrm{dr}}t}],$$
(104)
where $`\chi `$ is the Rabi frequency and $`\omega _{\mathrm{dr}}`$ is the driving frequency of the classical field. However as shown in Eq. (7) we can also write $`\widehat{H}_{\mathrm{sys}}`$ as $`\widehat{H}_\mathrm{\Omega }+\widehat{H}(t)`$. If $`\widehat{H}_\mathrm{\Omega }=\mathrm{\Omega }\widehat{\sigma }_z/2`$, then in the $`\mathrm{\Omega }`$ interaction picture gives

$$\widehat{H}_{\mathrm{int}}(t)=\mathrm{}\frac{\omega _0\mathrm{\Omega }}{2}\widehat{\sigma }_z+\mathrm{}\frac{\chi }{2}[\widehat{\sigma }e^{i(\omega _{\mathrm{dr}}t\mathrm{\Omega }t)}+\widehat{\sigma }^{}e^{i(\omega _{\mathrm{dr}}t\mathrm{\Omega }t)}],$$
(105)
For our purposes we assume $`\mathrm{\Omega }=\omega _{\mathrm{dr}}`$. So

$$\widehat{H}_{\mathrm{int}}(t)=\mathrm{}\frac{\mathrm{\Delta }}{2}\widehat{\sigma }_z+\mathrm{}\frac{\chi }{2}\widehat{\sigma }_x,$$
(106)
where $`\mathrm{\Delta }=\omega _0\mathrm{\Omega }`$ is the detuning.

For the TLA the reduced state can be written in terms of the real Bloch vector components $`x(t)`$, $`y(t)`$, $`z(t)`$ as

$$\rho _{\mathrm{red}}(t)=\frac{1}{2}[\widehat{I}+x(t)\widehat{\sigma }_x+y(t)\widehat{\sigma }_y+z(t)\widehat{\sigma }_z].$$
(107)
### V.2 Enlarged System Method

For the driven TLA with a memory function given by Eq. (100) the master equation for the enlarged systems is

$`d_tW_{\mathrm{red}}(t)`$ $`=`$ $`[{\displaystyle \frac{i\mathrm{\Delta }}{2}}\widehat{\sigma }_z{\displaystyle \frac{i\chi }{2}}\widehat{\sigma }_x+{\displaystyle \frac{\gamma \kappa }{4}}(\widehat{\sigma }\widehat{c}\widehat{\sigma }^{}\widehat{c}),`$ (108)
$`W_{\mathrm{red}}(t)]+\kappa 𝒟[\widehat{c}]W_{\mathrm{red}}(t).`$
Using $`\gamma =1`$, $`\kappa =1`$, $`\chi =5`$ and $`\mathrm{\Delta }=3`$ the reduced state is shown in Fig. 1. For this simple case it was noted that the truncation error involved in the enlarged system state method was negligible. Because of this we use this reduced state for comparison with the ensemble average of the non-Markovian SSEs.

### V.3 Coherent Unravelling-TLA

Applying the coherent non-Markovian SSE theory to the driven TLA, we find that we can rewrite the actual non-Markovian SSE as

$`d_t|\psi _z(t)`$ $`=`$ $`[i{\displaystyle \frac{\mathrm{\Delta }}{2}}\widehat{\sigma }_zi{\displaystyle \frac{\chi }{2}}\widehat{\sigma }_x(\widehat{\sigma }^{}\widehat{\sigma }^{}_t)`$ (109)
$`\times {}_{}{}^{(0)}\widehat{F}_{z}^{}(t)+(\widehat{\sigma }^{}\widehat{\sigma }^{}_t){}_{}{}^{(0)}\widehat{F}_{z}^{}(t)_t`$
$`+z^{}(t)(\widehat{\sigma }\widehat{\sigma }_t)]|\psi _z(t),`$
and the noise function for the TLA becomes

$$z(t)=z_\mathrm{\Lambda }(t)+_0^t\alpha (ts)\widehat{\sigma }_s𝑑s.$$
(110)
To calculate the complex amplitudes for the actual non-Markovian SSE we apply the system state $`|\psi _z(t)=C_e(t)|e+C_g(t)|g`$ to Eq. (109) and expand $`{}_{}{}^{(0)}\widehat{F}_{z}^{}(t)`$ as

$${}_{}{}^{(0)}\widehat{F}_{z}^{}(t)=\underset{m}{}\widehat{m}{}_{}{}^{(0)}F_{m,z}^{}(t)$$
(111)
where $`m=\{\sigma `$, $`\sigma ^{}`$, $`\sigma _z`$, $`I\}`$. This results in

$`d_tC_g`$ $`=`$ $`i{\displaystyle \frac{\mathrm{\Delta }}{2}}C_gi{\displaystyle \frac{\chi }{2}}C_e+z^{}C_e|C_e|^2{}_{}{}^{(0)}F_{\sigma ^{},z}^{}C_g^3C_{e}^{}{}_{}{}^{2}+{}_{}{}^{(0)}F_{\sigma ,z}^{}C_g|C_e|^2(1+|C_e|^2)`$ (112a)
$`{}_{}{}^{(0)}F_{\sigma _z,z}^{}C_g^2C_e^{}(1+2|C_e|^2)+{}_{}{}^{(0)}F_{I,z}^{}C_g^2C_e^{},`$
$`d_tC_e`$ $`=`$ $`i{\displaystyle \frac{\mathrm{\Delta }}{2}}C_ei{\displaystyle \frac{\chi }{2}}C_gz^{}C_e^2C_g^{}{}_{}{}^{(0)}F_{\sigma ,z}^{}C_e|C_g|^2(1+|C_e|^2)+{}_{}{}^{(0)}F_{\sigma ^{},z}^{}C_g^2C_e^{}|C_g|^2`$ (112b)
$`+{}_{}{}^{(0)}F_{\sigma _z,z}^{}C_g|C_g|^2(1+2|C_e|^2){}_{}{}^{(0)}F_{I,z}^{}C_g|C_g|^2.`$
In this equation the noise function is given by

$$z^{}(t)=z_\mathrm{\Lambda }^{}(t)+\frac{\gamma \kappa }{4}e^{\kappa t/2}_0^te^{\kappa s/2}C_g(s)C_e^{}(s)𝑑s,$$
(113)
where $`z_\mathrm{\Lambda }^{}(t)`$ is defined by the correlation

$$\stackrel{~}{E}[z_\mathrm{\Lambda }(t)z_\mathrm{\Lambda }^{}(s)]=\frac{\gamma \kappa }{4}e^{\kappa |ts|/2}.$$
(114)
This is generated by having $`z_\mathrm{\Lambda }^{}(t)`$ obey the following stochastic differential equation,

$$d_tz_\mathrm{\Lambda }(t)=\frac{\kappa }{2}z_\mathrm{\Lambda }(t)+\frac{\kappa }{2}\sqrt{\gamma }\zeta (t),$$
(115)
with $`z_\mathrm{\Lambda }^{}(0)`$ being a Gaussian random variable (GRV) satisfying

$$E[z_\mathrm{\Lambda }(0)z_\mathrm{\Lambda }^{}(0)]=\frac{\kappa \gamma }{4}$$
(116)
Here $`\zeta (t)`$ is standard complex white noise Gar83  and satisfies $`E[\zeta (t)\zeta ^{}(s)]=\delta (ts)`$.

#### V.3.1 $`0^{\mathrm{th}}`$ Order Approximation

For the simple memory function, $`J=1`$, which means $`{}_{}{}^{(0)}\widehat{F}_{z}^{}(t)={}_{}{}^{(0)}\widehat{F}_{z}^{(j)}(t)`$. The $`0^{\mathrm{th}}`$ order approximation occurs when we assume the form for $`{}_{}{}^{(0)}\widehat{F}_{z}^{}(t)`$ in Eq. (72). From Eq. (101) this implies

$${}_{}{}^{(0)}\widehat{F}_{z}^{}(t)=\frac{\gamma }{2}(1e^{\kappa t/2})\widehat{\sigma },$$
(117)
thus

$`{}_{}{}^{(0)}F_{\sigma ,z}^{}(t)`$ $`=`$ $`{\displaystyle \frac{\gamma }{2}}(1e^{\kappa t/2}),`$ (118a)
$`{}_{}{}^{(0)}F_{\sigma ^{},z}^{}(t)`$ $`=`$ $`{}_{}{}^{(0)}F_{\sigma _z,z}^{}(t)={}_{}{}^{(0)}F_{I,z}^{}(t)=0.`$ (118b)
#### V.3.2 $`1^{\mathrm{st}}`$ Order Approximation

The $`1^{\mathrm{st}}`$ first order approximation occurs when we assume a form for $`{}_{}{}^{(1)}\widehat{F}_{z}^{(j,k)}(t)`$, by Eqs. (73) and (101) this means

$${}_{}{}^{(1)}\widehat{F}_{z}^{}(t)=\frac{\gamma }{2}(1e^{\kappa t/2})[\widehat{\sigma },{}_{}{}^{(0)}\widehat{F}_{z}^{}(t)],$$
(119)
thus

$`{}_{}{}^{(1)}F_{\sigma ,z}^{}(t)`$ $`=`$ $`\gamma (1e^{\kappa t/2}){}_{}{}^{(0)}F_{\sigma _z,z}^{}(t),`$ (120a)
$`{}_{}{}^{(1)}F_{\sigma _z,z}^{}(t)`$ $`=`$ $`{\displaystyle \frac{\gamma }{2}}(1e^{\kappa t/2}){}_{}{}^{(0)}F_{\sigma ^{},z}^{}(t),`$ (120b)
$`{}_{}{}^{(1)}F_{\sigma ^{},z}^{}(t)`$ $`=`$ $`{}_{}{}^{(1)}F_{I,z}^{}(t)=0.`$ (120c)
The zero order functionals are found by applying the TLA operators to Eq. (III.1), giving

$`d_t{}_{}{}^{(0)}\widehat{F}_{z}^{}(t)`$ $`=`$ $`{\displaystyle \frac{\gamma \kappa }{4}}\widehat{\sigma }{\displaystyle \frac{\kappa }{2}}{}_{}{}^{(0)}\widehat{F}_{z}^{}(t)+z^{}(t)[\widehat{\sigma },{}_{}{}^{(0)}\widehat{F}_{z}^{}(t)]`$ (121)
$`i[{\displaystyle \frac{\mathrm{\Delta }}{2}}\widehat{\sigma }_z+{\displaystyle \frac{\chi }{2}}\widehat{\sigma }_x,{}_{}{}^{(0)}\widehat{F}_{z}^{}(t)][\widehat{\sigma }^{}{}_{}{}^{(0)}\widehat{F}_{z}^{}(t),`$
$`{}_{}{}^{(0)}\widehat{F}_{z}^{}(t)]\widehat{\sigma }^{}^{(1)}\widehat{F}_z(t).`$
Using Eq. (111) this gives the following four coupled nonlinear equations

$`d_t{}_{}{}^{(0)}F_{\sigma ,z}^{}(t)`$ $`=`$ $`\frac{1}{4}\gamma \kappa {\displaystyle \frac{\kappa }{2}}{}_{}{}^{(0)}F_{\sigma ,z}^{}(t)+i\mathrm{\Delta }{}_{}{}^{(0)}F_{\sigma ,z}^{}(t)i\chi {}_{}{}^{(0)}F_{\sigma _z,z}^{}(t)+2z^{}(t){}_{}{}^{(0)}F_{\sigma _z,z}^{}(t)`$ (122a)
$`+{}_{}{}^{(0)}F_{\sigma ,z}^{2}(t),`$
$`d_t{}_{}{}^{(0)}F_{\sigma ^{},z}^{}(t)`$ $`=`$ $`{\displaystyle \frac{\kappa }{2}}{}_{}{}^{(0)}F_{\sigma ^{},z}^{}(t)+i\chi {}_{}{}^{(0)}F_{\sigma _z,z}^{}(t)i\mathrm{\Delta }{}_{}{}^{(0)}F_{\sigma ^{},z}^{}(t)+2{}_{}{}^{(0)}F_{\sigma _z,z}^{}(t)[{}_{}{}^{(0)}F_{I,z}^{}(t){}_{}{}^{(0)}F_{\sigma _z,z}^{}(t)]`$ (122b)
$`{}_{}{}^{(0)}F_{\sigma ^{},z}^{}(t){}_{}{}^{(0)}F_{\sigma ,z}^{}(t)[{}_{}{}^{(1)}F_{I,z}^{}(t){}_{}{}^{(1)}F_{\sigma _z,z}^{}(t)],`$
$`d_t{}_{}{}^{(0)}F_{\sigma _z,z}^{}(t)`$ $`=`$ $`{\displaystyle \frac{\kappa }{2}}{}_{}{}^{(0)}F_{\sigma _z,z}^{}(t)+i{\displaystyle \frac{\chi }{2}}{}_{}{}^{(0)}F_{\sigma ^{},z}^{}(t)i{\displaystyle \frac{\chi }{2}}{}_{}{}^{(0)}F_{\sigma ,z}^{}(t){}_{}{}^{(0)}F_{\sigma ,z}^{}(t)[{}_{}{}^{(0)}F_{I,z}^{}(t){}_{}{}^{(0)}F_{\sigma _z,z}^{}(t)]`$ (122c)
$`z^{}(t){}_{}{}^{(0)}F_{\sigma ^{},z}^{}(t)\frac{1}{2}{}_{}{}^{(1)}F_{\sigma ,z}^{}(t),`$
$`d_t{}_{}{}^{(0)}F_{I,z}^{}(t)`$ $`=`$ $`{\displaystyle \frac{\kappa }{2}}{}_{}{}^{(0)}F_{I,z}^{}(t)\frac{1}{2}{}_{}{}^{(1)}F_{\sigma ,z}^{}(t).`$ (122d)
which can be solved in parallel with Eq. (112).

#### V.3.3 $`2^{\mathrm{nd}}`$ Order Approximation

The $`2^{\mathrm{nd}}`$ order approximation occurs when we assume a form for $`{}_{}{}^{(2)}\widehat{F}_{z}^{(j,k,l)}(t)`$, by Eqs. (III.1) and (101) this means

$${}_{}{}^{(2)}\widehat{F}_{z}^{}(t)=\frac{\gamma }{2}(1e^{\kappa t/2})[\widehat{\sigma },{}_{}{}^{(1)}\widehat{F}_{z}^{}(t)],$$
(123)
thus

$`{}_{}{}^{(2)}F_{\sigma ,z}^{}(t)`$ $`=`$ $`\gamma (1e^{\kappa t/2}){}_{}{}^{(1)}F_{\sigma _z,z}^{}(t),`$ (124a)
$`{}_{}{}^{(2)}F_{\sigma _z,z}^{}(t)`$ $`=`$ $`{\displaystyle \frac{\gamma }{2}}(1e^{\kappa t/2}){}_{}{}^{(1)}F_{\sigma ^{},z}^{}(t),`$ (124b)
$`{}_{}{}^{(2)}F_{\sigma ^{},z}^{}(t)`$ $`=`$ $`{}_{}{}^{(2)}F_{I,z}^{}(t)=0.`$ (124c)
The zero order functionals are given by Eqs. (122a) – (122d), however we now need equations for $`{}_{}{}^{(1)}\widehat{F}_{z}^{}(t)`$. The fist order functionals are found applying TLA operators to Eq. (66). With a memory function specified by Eq. (101) we get

$`d_t{}_{}{}^{(1)}\widehat{F}_{z}^{}(t)`$ $`=`$ $`{\displaystyle \frac{\gamma \kappa }{4}}[\widehat{\sigma },{}_{}{}^{(0)}\widehat{F}_{z}^{}(t)]\kappa {}_{}{}^{(1)}\widehat{F}_{z}^{}(t)+z^{}(t)`$ (125)
$`\times [\widehat{\sigma },{}_{}{}^{(1)}\widehat{F}_{z}^{}(t)]i[{\displaystyle \frac{\mathrm{\Delta }}{2}}\widehat{\sigma }_z+{\displaystyle \frac{\chi }{2}}\widehat{\sigma }_x,{}_{}{}^{(1)}\widehat{F}_{z}^{}(t)]`$
$`[\widehat{\sigma }^{}{}_{}{}^{(1)}\widehat{F}_{z}^{}(t),{}_{}{}^{(0)}\widehat{F}_{z}^{}(t)][\widehat{\sigma }^{}`$
$`\times {}_{}{}^{(0)}\widehat{F}_{z}^{}(t),{}_{}{}^{(1)}\widehat{F}_{z}^{}(t)]\widehat{\sigma }^{}^{(2)}\widehat{F}_z(t).`$
Using Eq. (124) this turns into the four equations

$`d_t{}_{}{}^{(1)}F_{\sigma ,z}^{}(t)`$ $`=`$ $`\frac{1}{2}\gamma \kappa {}_{}{}^{(0)}F_{\sigma _z,z}^{}(t)\kappa {}_{}{}^{(1)}F_{\sigma ,z}^{}(t)+i\mathrm{\Delta }{}_{}{}^{(1)}F_{\sigma ,z}^{}(t)i\chi {}_{}{}^{(1)}F_{\sigma _z,z}^{}(t)+2z^{}(t){}_{}{}^{(1)}F_{\sigma _z,z}^{}(t)`$ (126a)
$`+2{}_{}{}^{(0)}F_{\sigma ,z}^{}(t){}_{}{}^{(1)}F_{\sigma ,z}^{}(t),`$
$`d_t{}_{}{}^{(1)}F_{\sigma ^{},z}^{}(t)`$ $`=`$ $`\kappa {}_{}{}^{(1)}F_{\sigma ^{},z}^{}(t)+i\chi {}_{}{}^{(1)}F_{\sigma _z,z}^{}(t)i\mathrm{\Delta }{}_{}{}^{(1)}F_{\sigma ^{},z}^{}(t)+2{}_{}{}^{(1)}F_{\sigma _z,z}^{}(t)[{}_{}{}^{(0)}F_{I,z}^{}(t){}_{}{}^{(0)}F_{\sigma _z,z}^{}(t)]`$ (126b)
$`+2{}_{}{}^{(0)}F_{\sigma _z,z}^{}(t)[{}_{}{}^{(1)}F_{I,z}^{}(t){}_{}{}^{(1)}F_{\sigma _z,z}^{}(t)][{}_{}{}^{(1)}F_{\sigma ^{},z}^{}(t){}_{}{}^{(0)}F_{\sigma ,z}^{}(t)+{}_{}{}^{(0)}F_{\sigma ^{},z}^{}(t){}_{}{}^{(1)}F_{\sigma ,z}^{}(t)]`$
$`{}_{}{}^{(2)}F_{I,z}^{}(t)+{}_{}{}^{(2)}F_{\sigma _z,z}^{}(t),`$
$`d_t{}_{}{}^{(1)}F_{\sigma _z,z}^{}(t)`$ $`=`$ $`{\displaystyle \frac{\gamma \kappa }{4}}{}_{}{}^{(0)}F_{\sigma ^{},z}^{}(t)\kappa {}_{}{}^{(1)}F_{\sigma _z,z}^{}(t)+i{\displaystyle \frac{\chi }{2}}{}_{}{}^{(1)}F_{\sigma ^{},z}^{}(t)i{\displaystyle \frac{\chi }{2}}{}_{}{}^{(1)}F_{\sigma ,z}^{}(t){}_{}{}^{(1)}F_{\sigma ,z}^{}(t)[{}_{}{}^{(0)}F_{I,z}^{}(t)`$ (126c)
$`{}_{}{}^{(0)}F_{\sigma _z,z}^{}(t)]^{(0)}F_{\sigma ,z}(t)[{}_{}{}^{(1)}F_{I,z}^{}(t){}_{}{}^{(1)}F_{\sigma _z,z}^{}(t)]z^{}(t)^{(1)}F_{\sigma ^{},z}(t)`$
$`\frac{1}{2}{}_{}{}^{(2)}F_{\sigma ,z}^{}(t),`$
$`d_t{}_{}{}^{(1)}F_{I,z}^{}(t)`$ $`=`$ $`\kappa {}_{}{}^{(1)}F_{I,z}^{}(t)\frac{1}{2}{}_{}{}^{(2)}F_{\sigma ,z}^{}(t).`$ (126d)
To illustrate how accurate our perturbation method is, the difference between the reduced state calculated via the enlarged system method and the ensemble average from the coherent non-Markovian SSE is plotted in Fig. 2. The dotted line corresponds to the $`0^{\mathrm{th}}`$ order perturbation, the dashed is the $`1^{\mathrm{st}}`$ and the solid is the $`2^{\mathrm{nd}}`$. It is observed that the $`1^{\mathrm{st}}`$ and $`2^{\mathrm{nd}}`$ order perturbation are a lot more accurate then the $`0^{\mathrm{th}}`$ order perturbation. However, it can be seen that the $`2^{\mathrm{nd}}`$ order perturbation is not necessarily more accurate than the $`1^{\mathrm{st}}`$ order perturbation. This suggest that our perturbation method is an asymptotic expansion rather than a convergent series.

### V.4 Quadrature Unravelling-TLA

For the quadrature unravelling the actual non-Markovian SSE is

$`d_t|\psi _z(t)`$ $`=`$ $`[i{\displaystyle \frac{\mathrm{\Delta }}{2}}\widehat{\sigma }_zi{\displaystyle \frac{\chi }{2}}\widehat{\sigma }_x(\widehat{\sigma }_x\widehat{\sigma }_x_t){}_{}{}^{(0)}\widehat{Q}_{z}^{}(t)`$ (127)
$`+(\widehat{\sigma }_x\widehat{\sigma }_x_t){}_{}{}^{(0)}\widehat{Q}_{z}^{}(t)_t`$
$`+z(t)(\widehat{\sigma }\widehat{\sigma }_t)]|\psi _z(t),`$
and the noise function for the TLA is

$$z(t)=z_\mathrm{\Lambda }(t)+_0^t\beta (ts)\widehat{\sigma }_x_s𝑑s.$$
(128)
Again the coherent case we can calculate the complex amplitude equation via applying the state $`|\psi _z(t)=C_e(t)|e+C_g(t)|g`$ to Eq. (127) and expanding $`{}_{}{}^{(0)}\widehat{Q}_{z}^{}(t)`$ as

$${}_{}{}^{(0)}\widehat{Q}_{z}^{}(t)=\underset{m}{}\widehat{m}{}_{}{}^{(0)}Q_{m,z}^{}(t)$$
(129)
where $`m=\{\sigma `$, $`\sigma ^{}`$, $`\sigma _z`$, $`I\}`$. This results in a coupled set of differential equations for $`C_e(t)`$ and $`C_g(t)`$ that depend on $`{}_{}{}^{(0)}Q_{m,z}^{}(t)`$ and $`z(t)`$. In these equations the real-valued noise is given by

$`z(t)`$ $`=`$ $`z_\mathrm{\Lambda }(t)+{\displaystyle \frac{\gamma \kappa }{4}}e^{\kappa t/2}{\displaystyle _0^t}e^{\kappa s/2}[C_g(s)C_e^{}(s)`$ (130)
$`+C_g^{}(s)C_e(s)]ds,`$
where $`z_\mathrm{\Lambda }(t)`$ is found by

$$\stackrel{~}{E}[z_\mathrm{\Lambda }(t)z_\mathrm{\Lambda }(s)]=\frac{\gamma \kappa }{4}e^{\kappa |ts|/2}.$$
(131)
This is generated by

$$d_tz_\mathrm{\Lambda }(t)=\frac{\kappa }{2}z_\mathrm{\Lambda }(t)+\frac{\kappa }{2}\sqrt{\gamma }\xi (t)$$
(132)
with $`z_\mathrm{\Lambda }(0)`$ being a GRV satisfying $`E[z_\mathrm{\Lambda }(0)z_\mathrm{\Lambda }^{}(0)]=\kappa \gamma /4`$. Here $`\xi (t)`$ is standard white noise and satisfies $`E[\xi (t)\xi ^{}(s)]=\delta (ts)`$ Gar83 .

#### V.4.1 $`0^{\mathrm{th}}`$ Order Approximation

The situation is greatly simplified with the memory function in Eq. (100), as $`\beta (t,s)=\beta ^{(j,cos)}(t,s)=\beta ^{(j,cos)}(t,s)`$, which in turn implies $`{}_{}{}^{(0)}\widehat{Q}_{z}^{}(t)={}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{cos})}(t)={}_{}{}^{(0)}\widehat{Q}_{z}^{(j,\mathrm{sin})}(t)`$.

The $`0^{\mathrm{th}}`$ order approximation is to set

$${}_{}{}^{(0)}\widehat{Q}_{z}^{}(t)=\frac{\gamma }{2}(1e^{\kappa t/2})\widehat{\sigma },$$
(133)
thus

$`{}_{}{}^{(0)}Q_{\sigma ,z}^{}(t)`$ $`=`$ $`{\displaystyle \frac{\gamma }{2}}(1e^{\kappa t/2}),`$ (134a)
$`{}_{}{}^{(0)}Q_{\sigma ^{},z}^{}(t)`$ $`=`$ $`{}_{}{}^{(0)}Q_{\sigma _z,z}^{}(t)={}_{}{}^{(0)}Q_{I,z}^{}(t)=0.`$ (134b)
#### V.4.2 $`1^{\mathrm{th}}`$ Order Approximation

The first order approximation is to set

$${}_{}{}^{(1)}\widehat{Q}_{z}^{}(t)=\frac{\gamma }{2}(1e^{\kappa t/2})[\widehat{\sigma },{}_{}{}^{(0)}\widehat{Q}_{z}^{}(t)]$$
(135)
thus

$`{}_{}{}^{(1)}Q_{\sigma ,z}^{}(t)`$ $`=`$ $`\gamma (1e^{\kappa t/2}){}_{}{}^{(0)}Q_{\sigma _z,z}^{}(t),`$ (136a)
$`{}_{}{}^{(1)}Q_{\sigma _z,z}^{}(t)`$ $`=`$ $`{\displaystyle \frac{\gamma }{2}}(1e^{\kappa t/2}){}_{}{}^{(0)}Q_{\sigma ^{},z}^{}(t),`$ (136b)
$`{}_{}{}^{(1)}Q_{\sigma ^{},z}^{}(t)`$ $`=`$ $`{}_{}{}^{(1)}Q_{I,z}^{}(t)=0.`$ (136c)
The $`0^{\mathrm{th}}`$ order functionals are found by applying TLA operators to Eq. (84). With the simple memory function this gives

$`d_t{}_{}{}^{(0)}\widehat{Q}_{z}^{}(t)`$ $`=`$ $`{\displaystyle \frac{\gamma \kappa }{4}}\widehat{\sigma }{\displaystyle \frac{\kappa }{2}}{}_{}{}^{(0)}\widehat{Q}_{z}^{}(t)+z(t)[\widehat{\sigma },{}_{}{}^{(0)}\widehat{Q}_{z}^{}(t)]`$ (137)
$`i[{\displaystyle \frac{\mathrm{\Delta }}{2}}\widehat{\sigma }_z+{\displaystyle \frac{\chi }{2}}\widehat{\sigma }_x,{}_{}{}^{(0)}\widehat{Q}_{z}^{}(t)]`$
$`[\widehat{\sigma }_x{}_{}{}^{(0)}\widehat{Q}_{z}^{}(t),{}_{}{}^{(0)}\widehat{Q}_{z}^{}(t)]`$
$`\widehat{\sigma }_x{}_{}{}^{(1)}\widehat{Q}_{z}^{}(t).`$
Using Eq. (129) this gives,

$`d_t{}_{}{}^{(0)}Q_{\sigma ,z}^{}(t)`$ $`=`$ $`\frac{1}{4}\gamma \kappa {\displaystyle \frac{\kappa }{2}}{}_{}{}^{(0)}Q_{\sigma ,z}^{}(t)+i\mathrm{\Delta }{}_{}{}^{(0)}Q_{\sigma ,z}^{}(t)i\chi {}_{}{}^{(0)}Q_{\sigma _z,z}^{}(t)+2z(t){}_{}{}^{(0)}Q_{\sigma _z,z}^{}(t)+{}_{}{}^{(0)}Q_{\sigma ,z}^{2}(t)`$ (138a)
$`2{}_{}{}^{(0)}Q_{\sigma _z,z}^{}(t)[{}_{}{}^{(0)}Q_{I,z}^{}(t)+{}_{}{}^{(0)}Q_{\sigma _z,z}^{}(t)]{}_{}{}^{(0)}Q_{\sigma ^{},z}^{}(t){}_{}{}^{(0)}Q_{\sigma ,z}^{}(t)`$
$`[{}_{}{}^{(1)}Q_{I,z}^{}(t)+{}_{}{}^{(1)}Q_{\sigma _z,z}^{}(t)],`$
$`d_t{}_{}{}^{(0)}Q_{\sigma ^{},z}^{}(t)`$ $`=`$ $`{\displaystyle \frac{\kappa }{2}}{}_{}{}^{(0)}Q_{\sigma ^{},z}^{}(t)+i\chi {}_{}{}^{(0)}Q_{\sigma _z,z}^{}(t)i\mathrm{\Delta }{}_{}{}^{(0)}Q_{\sigma ^{},z}^{}(t)+2{}_{}{}^{(0)}Q_{\sigma _z,z}^{}(t)[{}_{}{}^{(0)}Q_{I,z}^{}(t){}_{}{}^{(0)}Q_{\sigma _z,z}^{}(t)]`$ (138b)
$`{}_{}{}^{(0)}Q_{\sigma ^{},z}^{}(t){}_{}{}^{(0)}Q_{\sigma ,z}^{}(t)+{}_{}{}^{(0)}Q_{\sigma ^{},z}^{2}(t){}_{}{}^{(1)}Q_{I,z}^{}(t)+{}_{}{}^{(1)}Q_{\sigma _z,z}^{}(t),`$
$`d_t{}_{}{}^{(0)}Q_{\sigma _z,z}^{}(t)`$ $`=`$ $`{\displaystyle \frac{\kappa }{2}}{}_{}{}^{(0)}Q_{\sigma _z,z}^{}(t)+i{\displaystyle \frac{\chi }{2}}{}_{}{}^{(0)}Q_{\sigma ^{},z}^{}(t)i{\displaystyle \frac{\chi }{2}}{}_{}{}^{(0)}Q_{\sigma ,z}^{}(t){}_{}{}^{(0)}Q_{\sigma ,z}^{}(t)[{}_{}{}^{(0)}Q_{I,z}^{}(t){}_{}{}^{(0)}Q_{\sigma _z,z}^{}(t)]`$ (138c)
$`+{}_{}{}^{(0)}Q_{\sigma ^{},z}^{}(t)[{}_{}{}^{(0)}Q_{I,z}^{}(t)+{}_{}{}^{(0)}Q_{\sigma _z,z}^{}(t)]z(t){}_{}{}^{(0)}Q_{\sigma ^{},z}^{}(t)`$
$`\frac{1}{2}[{}_{}{}^{(1)}Q_{\sigma ,z}^{}(t){}_{}{}^{(1)}Q_{\sigma ^{},z}^{}(t)],`$
$`d_t{}_{}{}^{(0)}Q_{I,z}^{}(t)`$ $`=`$ $`{\displaystyle \frac{\kappa }{2}}{}_{}{}^{(0)}Q_{I,z}^{}(t)\frac{1}{2}[{}_{}{}^{(1)}Q_{\sigma ,z}^{}(t)+{}_{}{}^{(1)}Q_{\sigma ^{},z}^{}(t)].`$ (138d)
which can be solved in parallel with $`C_e(t)`$ and $`C_g(t)`$.

To illustrate how accurate our perturbation method is for the quadrature unravelling. Fig. 3 shows the difference between the reduced state calculated via the enlarged system method and the ensemble average from the quadrature non-Markovian SSEs for the $`0^{\mathrm{th}}`$ (dotted) and $`1^{\mathrm{st}}`$ (dashed) order perturbation. As in the coherent case we find the $`1^{\mathrm{st}}`$ order perturbation is more accurate then the $`0^{\mathrm{th}}`$.

## VI Post-Markovian perturbation

In this section we extend the YDGS post-Markovian perturbation YuDioGisStr99  to include the quadrature unraveling and compare the post-Markovian method with our perturbation method.

The basis idea behind their perturbation method is to expand the operators $`{}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s)`$ in powers of $`(ts)`$ around the point $`t=s`$ (this is why it is called the post Markovian perturbation). That is

$`{}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s)`$ $`=`$ $`{}_{}{}^{(0)}\widehat{f}_{z}^{}(s,s)+[_t{}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s)|_{t=s}](ts)`$ (139)
$`+\frac{1}{2}[_t^2{}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s)|_{t=s}](ts)^2+\mathrm{},`$
where $`{}_{}{}^{(0)}\widehat{f}_{z}^{}(s,s)=\widehat{L}`$. To find the first order term we simply evaluate Eq. (57) at $`t=s`$

$`_t{}_{}{}^{(0)}\widehat{f}_{z}^{}(t,s)|_{t=s}`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}[\widehat{H}_{\mathrm{int}}(s),\widehat{L}][\widehat{L}^{}{}_{}{}^{(0)}\widehat{F}_{z}^{}(s),\widehat{L}]`$ (140)
$`\widehat{L}^{}[\widehat{L},{}_{}{}^{(0)}\widehat{F}_{z}^{}(s)].`$
Thus the functional $`{}_{}{}^{(0)}\widehat{F}_{z}^{}(t)`$ for this perturbation is given by

$`{}_{}{}^{(0)}\widehat{F}_{z}^{}(t)`$ $`=`$ $`g_0(t)\widehat{L}g_1(t){\displaystyle \frac{i}{\mathrm{}}}[\widehat{H}_{\mathrm{int}}(t),\widehat{L}]`$ (141)
$`{\displaystyle _0^t}\alpha (ts)(ts)[\widehat{L}^{}{}_{}{}^{(0)}\widehat{F}_{z}^{}(s),\widehat{L}]𝑑s`$
$`{\displaystyle _0^t}\alpha (ts)(ts)\widehat{L}^{}[\widehat{L},{}_{}{}^{(0)}\widehat{F}_{z}^{}(s)]𝑑s,`$
where

$`g_0(t)`$ $`=`$ $`{\displaystyle _0^t}\alpha (ts)𝑑s,`$ (142)
$`g_1(t)`$ $`=`$ $`{\displaystyle _0^t}\alpha (ts)(ts)𝑑s.`$ (143)
This equation can not be solved without the initial condition $`d_t{}_{}{}^{(0)}\widehat{F}_{z}^{}(0)`$. However if we make the approximate $`{}_{}{}^{(0)}\widehat{F}_{z}^{}(s)=_0^s\alpha (su)\widehat{L}𝑑u`$, Eq. (141) becomes

$${}_{}{}^{(0)}\widehat{F}_{z}^{}(t)=g_0(t)\widehat{L}g_1(t)\frac{i}{\mathrm{}}[\widehat{H}_{\mathrm{int}}(t),\widehat{L}]g_2(t)[\widehat{L}^{}\widehat{L},\widehat{L}],$$
(144)
where

$$g_2(t)=_0^t_0^s\alpha (ts)\alpha (su)(ts)𝑑u𝑑s,$$
(145)
which can be solved. The same could be done for the second order terms, but as well as making an approximation for $`{}_{}{}^{(0)}\widehat{F}_{z}^{}(s)`$ we would need to approximate $`d_s{}_{}{}^{(0)}\widehat{F}_{z}^{}(s)`$. For the purpose of this paper we will only go to first order.

To extend the idea to the quadrature case we Taylor expand the operator $`{}_{}{}^{(0)}\widehat{q}_{z}^{}(t,s)`$ in powers of $`(ts)`$ around the point $`t=s`$. To find the first order term we simply evaluate Eq. (83) at $`t=s`$. With the approximation $`{}_{}{}^{(0)}\widehat{Q}_{z}^{}(s)=_0^s\beta (su)\widehat{L}𝑑u`$ we get

$${}_{}{}^{(0)}\widehat{Q}_{z}^{}(t)=h_0(t)\widehat{L}h_1(t)\frac{i}{\mathrm{}}[\widehat{H}_{\mathrm{int}}(t),\widehat{L}]h_2(t)[\widehat{L}_x\widehat{L},\widehat{L}].$$
(146)
where

$`h_0(t)`$ $`=`$ $`{\displaystyle _0^t}\beta (ts)𝑑s,`$ (147)
$`h_1(t)`$ $`=`$ $`{\displaystyle _0^t}\beta (ts)(ts)𝑑s,`$ (148)
$`h_2(t)`$ $`=`$ $`{\displaystyle _0^t}{\displaystyle _0^s}\beta (ts)\beta (su)(ts)𝑑u𝑑s.`$ (149)
For the simple TLA system it is easy to generate these approximate expressions for $`{}_{}{}^{(0)}\widehat{F}_{z}^{}(t)`$ and $`{}_{}{}^{(0)}\widehat{Q}_{z}^{}(t)`$ for all time, hence we can obtain solution to the non-Markovian SSE. To compare YDGS post-Markovian non-Markovian SSE method with our perturbation method, we again plot the difference between YDGS method (when 1000 trajectories where used) and the enlarged systems method. The results of this are shown in Fig. 4, where it is observed that YDGS first order perturbation has a greater error than our perturbation method (Figs. 2 and 3). This is perhaps not surprising, as the system we modelled has $`\kappa =1`$, which implies it is very non-Markovian. Since one of the requirements of YDGS perturbation method is for the environment to be close the Markovian regime one would expect their method to fail in this regime.

In Ref. YuDioGisStr99  YDGS suggest an alternative perturbation method. The functional operator $`{}_{}{}^{(0)}\widehat{F}_{z}^{}(t)`$, which equals $`\overline{O}_z(t)`$ in their notation, is expanded by the functional expansion

$`\overline{O}_z(t)`$ $`=`$ $`\overline{O}^{(0)}(t)+{\displaystyle _0^t}\overline{O}^{(1)}(t,v)z(v)𝑑v+{\displaystyle _0^t}{\displaystyle _0^t}\overline{O}^{(2)}(t,v_1,v_2)z(v_1)z(v_2)𝑑v_1𝑑v_2+\mathrm{}`$ (150)
$`+{\displaystyle _0^t}\mathrm{}{\displaystyle _0^t}\overline{O}^{(n)}(t,v_1,\mathrm{},v_n)z(v_1)\mathrm{}z(v_n)𝑑v_1\mathrm{}𝑑v_n+\mathrm{},`$
It can be shown that one can establish a set of coupled differential equations for these operators provided $`\alpha (ts)`$ is given by Eq. (47). To truncated this perturbation at $`\overline{O}^{(n)}`$ one has to assume a value $`\overline{O}^{(n+1)}`$. It turns out that for all operators $`\overline{O}^{(n)}`$ other then $`\overline{O}^{(0)}`$ the only reason the operators change from their initial value $`0`$ at $`t=0`$ is if the assumed $`\overline{O}^{(n+1)}`$ is nonzero. This suggest that this method is highly dependent on the assumed value for $`\overline{O}^{(n+1)}`$.

## VII Conclusions

In this paper we presented a perturbation method for solving the coherent and quadrature non-Markovian SSEs. This perturbation method is easily extended to any order and is not limited to the post Markovian regime. However, the environment is restricted such that it has a correlation function satisfying Eq. (47). As shown in Ref. SteIma96  most non-Markovian environments can be simulated via this correlation function with a relative small $`J`$. This suggest that this perturbation method might be useful for simulating non-Markovian evolution for $`\rho _{\mathrm{red}}(t)`$.

One appealing feature of this method is that it provides a perturbative solution for $`\rho _{\mathrm{red}}(t)`$ which is positive by definition. However there is another method, namely Imamo$`\overline{\mathrm{g}}`$lu’s enlarged system method Ima94 ; SteIma96 , which provides a better solution for $`\rho _{\mathrm{red}}(t)`$. Imamo$`\overline{\mathrm{g}}`$lu’s enlarged system method requires fewer coupled differential equations to solve and the only approximation comes in by a truncation of the Hilbert space of the fictitious modes. As one increases the basis size for these modes this method will converge to the correct solution. By contrast, convergence has not been shown for our method.

This does not mean that our method is useless, as the primary interest in our method is not to simulate $`\rho _{\mathrm{red}}(t)`$, but to simulate the non-Markovian SSEs. This is interesting as a continuous in time interpretation of non-Markovian SSEs is not clear. In Ref. GamWis02  we showed that these non-Markovian SSE under standard quantum measurement theory do not have a continuous measurement interpretation. However Loubenets in Ref. Lou01 ; BarLou02  claimed that she has developed a new framework for continuous quantum measurements in which non-Markovian SSEs represent the evolution of a system state which is continuously monitored.

Future work on this topic is to look into this question. Another question that needs answering is whether it is possible to derive non-Markovian SSE based on a discrete basis such as photon number. We believe this question and the previous question will be related. Finally, there is the possible application of our method to strongly non-Markovian systems such as an atom laser HopMoyColSav00  or photon emission in a photonic bad-gap material Joh84 ; BayLamMol97 .

## Appendix A Derivation of $`{}_{}{}^{(0)}\widehat{f}_{z}^{}(t,t)=\widehat{L}`$

To show that $`{}_{}{}^{(0)}\widehat{f}_{z}^{}(t,t)=\widehat{L}`$ we start by discretizing the functional derivative. We divide the range $`[0,t)`$ into $`N`$ intervals of width $`\mathrm{\Delta }t`$, so the change in $`|\stackrel{~}{\psi }_z(t)`$ is

$`\delta |\stackrel{~}{\psi }_z(t)`$ $`=`$ $`{\displaystyle _0^t}{\displaystyle \frac{\delta |\stackrel{~}{\psi }_z(t)}{\delta z^{}(s)}}\delta z^{}(s)𝑑s`$ (151)
$`=`$ $`{\displaystyle \underset{i=0}{\overset{N1}{}}}\mathrm{\Delta }t\left[{\displaystyle \frac{|\stackrel{~}{\psi }_z(t_N)}{z^{}(t_i)\mathrm{\Delta }t}}\right]dz^{}(t_i),`$
thus

$`{\displaystyle \frac{\delta }{\delta z^{}(s)}}|\stackrel{~}{\psi }_z(t)={\displaystyle \frac{|\stackrel{~}{\psi }_z(t_N)}{z^{}(t_i)\mathrm{\Delta }t}},`$ (152)
if $`s`$ ($`t_i`$) is less than $`t`$ $`(t_N)`$, which is the only situation we are interested in, then taking the limit that $`st`$ ($`t_i=t_{N1}`$) this becomes

$$\underset{st}{lim}\frac{\delta |\stackrel{~}{\psi }_z(t)}{\delta z^{}(s)}=\frac{[|\stackrel{~}{\psi }_z(t_{N1})+\mathrm{\Delta }t_t|\stackrel{~}{\psi }_z(t_{N1})]}{z^{}(t_{N1})\mathrm{\Delta }t}.$$
(153)
Discretizing Eq. (24) we get

$`_t|\stackrel{~}{\psi }_z(t_{N1})`$ $`=`$ $`[{\displaystyle \frac{i}{\mathrm{}}}\widehat{H}_{\mathrm{int}}(t_{N1})+z^{}(t_{N1})\widehat{L}`$ (154)
$`\widehat{L}^{}{\displaystyle \underset{j=0}{\overset{N2}{}}}\alpha (t_{N1}t_j){\displaystyle \frac{}{z^{}(t_j)}}]`$
$`\times |\stackrel{~}{\psi }_z(t_{N1}).`$
Substituting this into Eq. (153) and using the fact that the state at time $`t_{N1}`$ only depends on the noise at time less then $`t_{N1}`$, we get the limit

$$\frac{\delta |\stackrel{~}{\psi }_z(t)}{\delta z^{}(t)}\widehat{L}|\stackrel{~}{\psi }_z(t).$$
(155)
Thus by Eq. (25) $`{}_{}{}^{(0)}\widehat{f}_{z}^{}(t,t)=\widehat{L}`$.