Time-local unraveling of non-Markovian stochastic Schr\"odinger equations

TL;DR

This paper introduces a non-perturbative numerical method for solving non-Markovian stochastic Schr"odinger equations, making them more accessible and applicable in quantum mechanics research.

Contribution

It proposes a novel approach to rewrite NMSSE as averages over auxiliary equations, enabling accurate simulations for both linear and non-linear cases.

Findings

Provides a non-perturbative simulation algorithm for linear NMSSE.

Extends method to non-linear NMSSE for isotropic complex noises.

Enables direct sampling of norm-preserving NMSSE solutions.

Abstract

Non-Markovian stochastic Schr\"odinger equations (NMSSE) are important tools in quantum mechanics, from the theory of open systems to foundations. Yet, in general, they are but formal objects: their solution can be computed numerically only in some specific cases or perturbatively. This article is focused on the NMSSE themselves rather than on the open-system evolution they unravel and aims at making them less abstract. Namely, we propose to write the stochastic realizations of linear NMSSE as averages over the solutions of an auxiliary equation with an additional random field. Our method yields a non-perturbative numerical simulation algorithm for generic linear NMSSE that can be made arbitrarily accurate for reasonably short times. For isotropic complex noises, the method extends from linear to non-linear NMSSE and allows to sample the solutions of norm-preserving NMSSE directly.

Equations

$$\rho (t)=$$

$$=$$

$$\displaystyle-\theta_{\tau s}A^{i}_{L}(\tau)A^{j}_{L}(s)-\theta_{s\tau}A^{j}_{R}(s)A^{i}_{R}(\tau)\Big{]}\bigg{\}}\,\rho(0)$$

$$\begin{array}[]{ll}A_{L}\cdot\rho=\hat{A}\rho\;;&A_{R}\cdot\rho=\rho\hat{A}\end{array}.$$

$$H_{\text{int}}={\hat{A}}^k\otimes \int \mathrm{d}\omega {\kappa }_k^{\ell }{a}_{\ell }(\omega )+{\kappa }_k^{\ell \ast }{a}_{\ell }^{†}(\omega ),$$

$$D_{ij}(\tau ,s)=\int \mathrm{d}\omega {\kappa }_i^k(\omega ){\kappa }_j^{\ell \ast }(\omega )e^{-i\omega (\tau -s)}.$$

$$\rho(t)=\mathds{E}_{\xi}\big{[}\,|\psi_{\xi}(t)\rangle\langle\psi_{\xi}(t)|\,\big{]},$$

$$\begin{split}\mathds{E}_{\xi}\left[\xi_{i}(\tau)\xi^{*}_{j}(s)\right]&=C_{ij}(\tau,s)\\ \mathds{E}_{\xi}\big{[}\xi_{i}(\tau)\xi_{j}\,(s)\big{]}&=S_{ij}(\tau,s)\end{split}$$

$$\text{mathds}{E}_{\xi }\left[\exp \left\{-i{\int }_0^t\mathrm{d}s(a_s^k{\xi }_k(s)-b_s^k{\xi }_k^{\ast }(s))\right\}\right]=\exp \left\{{\int }_0^t{\int }_0^t\mathrm{d}\tau \mathrm{d}s{C}_{ij}(\tau ,s)a_{\tau }^i{b}_s^j-\frac{S_{ij}(\tau ,s)a_{\tau }^i{a}_s^j+S_{ij}^{\ast }(\tau ,s)b_{\tau }^i{b}_s^j}{2}\right\},$$

$$\text{mathds}{E}_{\xi }\left[\mathcal{T}\exp \left\{-i{\int }_0^t\mathrm{d}s({\hat{a}}_s^k{\xi }_k(s)-{\hat{b}}_s^k{\xi }_k^{\ast }(s))\right\}\right]=\mathcal{T}\exp \left\{{\int }_0^t{\int }_0^t\mathrm{d}\tau \mathrm{d}s{C}_{ij}(\tau ,s){\hat{a}}_{\tau }^i{\hat{b}}_s^j-\frac{S_{ij}(\tau ,s){\hat{a}}_{\tau }^i{\hat{a}}_s^j+S_{ij}^{\ast }(\tau ,s){\hat{b}}_{\tau }^i{\hat{b}}_s^j}{2}\right\}.$$

$$|{\psi }_{\xi }(t)⟩=$$

$$:=$$

$$\displaystyle\!\!\!\!-\!\!\int_{0}^{t}\!\!\int_{0}^{t}\!\!\mathrm{d}\tau\,\mathrm{d}s\,\theta_{\tau s}\,[D-S]_{ij}(\tau,s)\hat{A}^{i}(\tau)\hat{A}^{j}(s)\bigg{\}}|\psi_{0}\rangle,$$

$$\begin{split}\frac{\mathrm{d}}{\mathrm{d}t}|\psi_{\xi}(t)\rangle=&-i\hat{A}^{i}(t)\bigg{[}\xi_{i}(t)\\ &+\int_{0}^{t}\!\!\mathrm{d}s\,[D-S]_{ij}(t,s)\frac{\delta}{\delta\xi_{j}(s)}\bigg{]}|\psi_{\xi}(t)\rangle.\end{split}$$

$$|\psi_{\xi}(t)\rangle=\mathds{E}_{\eta}\big{[}\,|\psi_{\xi,\eta}(t)\rangle\,\big{]}$$

$${\hat{G}}_{\xi }(t)=\text{mathds}{E}_{\eta }\left[\mathcal{T}\exp \left\{-i{\int }_0^t\mathrm{d}s{\hat{A}}^k(s)[{\xi }_k(s)+{\eta }_k(s)]\right\}\right].$$

$$\begin{array}{rl}\text{mathds}{E}_{\eta }[{\eta }_i(\tau ){\eta }_j(s)]& =K_{ij}(\tau ,s)\\ \text{mathds}{E}_{\eta }[{\eta }_i(\tau ){\eta }_j^{\ast }(s)]& =J_{ij}(\tau ,s)\end{array}$$

$$K_{ij}(\tau ,s)={\theta }_{\tau s}[D-S{]}_{ij}(\tau ,s)+{\theta }_{s\tau }[D-S{]}_{ji}(s,\tau )$$

$$\Gamma=\left(\begin{array}[]{cc}J&K\\ K^{*}&J^{*}\end{array}\right)$$

$$\frac{\mathrm{d}}{\mathrm{d}t}|{\psi }_{\xi ,\eta }(t)⟩=-i{\hat{A}}^k(t)[{\xi }_k(t)+{\eta }_k(t)]|{\psi }_{\xi ,\eta }(t)⟩.$$

$$\rho(t)=\mathds{E}_{\xi}\Big{[}\,\mathds{E}_{\eta}\big{[}\,|\psi_{\xi,\eta}(t)\rangle\,\big{]}\;\,\mathds{E}_{\eta}\big{[}\,\langle\psi_{\xi,\eta}(t)|\,\big{]}\,\Big{]}.$$

$$\rho(t)=\mathds{E}_{\xi,\eta^{(1)},\eta^{(2)}}\big{[}\,|\psi_{\xi,\eta^{(1)}}(t)\rangle\langle\psi_{\xi,\eta^{(2)}}(t)|\,\big{]},$$

$$|{\tilde{\psi }}_{\xi }(t)⟩=\frac{|{\psi }_{\xi }(t)⟩}{\sqrt{⟨{\psi }_{\xi }(t)|{\psi }_{\xi }(t)⟩}}$$

$$\mathrm{d}{\mu }_t(\xi )=⟨{\psi }_{\xi }(t)|{\psi }_{\xi }(t)⟩\mathrm{d}{\mu }_{\circ }(\xi )$$

$$\mathds{E}^{t}_{\xi}[f(\xi)]=\mathds{E}_{\xi}\Big{[}f\big{(}\xi^{[t]}(\xi)\big{)}\Big{]}.$$

$${\xi }_k^{[t]}(u)={\xi }_k(u)+i{\int }_0^t\mathrm{d}s{D}_{k\ell }(u,s)⟨A^{\ell }(s){⟩}_s,$$

$${\xi }_k^{[v+\mathrm{d}v]}(u)={\xi }_k^{[v]}(u)+i{D}_{k,\ell }(u,v)⟨A^{\ell }(v){⟩}_v\mathrm{d}v+O(\mathrm{d}{v}^2),$$

$$\text{mathds}E$$

$$=\text{mathds}{E}_{\xi }^t\left[\phi \left(\frac{|{\psi }_{\xi }(t)⟩}{\sqrt{⟨{\psi }_{\xi }(t)|{\psi }_{\xi }(t)⟩}}\right)\right]$$

$$=\text{mathds}{E}_{\xi }\left[\phi \left(\frac{|{\psi }_{\xi }(t)⟩}{\sqrt{⟨{\psi }_{\xi }(t)|{\psi }_{\xi }(t)⟩}}\right)⟨{\psi }_{\xi }(t)|{\psi }_{\xi }(t)⟩\right].$$

$$\frac{\mathrm{d}}{\mathrm{d}t}|{\psi }_{\xi }(t)⟩=-i\hat{A}\left[\xi (t)+{\int }_0^t\mathrm{d}sD(t,s)\frac{\delta }{\delta \xi (s)}\right]|{\psi }_{\xi }(t)⟩.$$

$$\frac{\delta }{\delta \xi (s)}|{\psi }_{\xi }(t)⟩=-i\hat{A}|{\psi }_{\xi }(t)⟩.$$

Full text

1 Introduction

Stochastic pure state representations have a wide range of applications in quantum mechanics: they serve as computational tools to unravel open-system evolutions [4, 19, 32], as modelling tools to describe continuous measurement situations [22, 36] and as foundational tools to solve the measurement problem in models where the superposition principle breaks down [2, 3]. In the Markovian limit, the solutions of stochastic Schrödinger equations can be efficiently computed numerically, justifying their wide use in the three aforementioned fields. In the non-Markovian case however, there exists no general purpose method to compute a realization of the random pure state process: plotting a single trajectory is in general impossible. In this article, we propose to write the solutions of general stochastic Schrödinger equations as averages over the solutions of time-local stochastic Schrödinger equations with an auxiliary noise. We shall use the same trick that usually allows to replace openness by stochasticity, this time to get rid of non-Markovianity.

Although stochastic Schrödinger equations have a long history dating back to the eighties [18, 8], the field took off when their potential in numerics was understood by Dalibard, Castin and Mølmer [4, 26] (see also Dum, Zoller and Ritsch [12]) in the jump case and by Gisin and Percival [19] in the diffusive case. The connection with actual quantum trajectories coming from realistic measurement setups was understood roughly at the same time by Milburn and Wiseman [35] (see e.g. [22] for an introduction). Interestingly, the introduction of continuous stochastic pure state equations in foundations is slightly anterior with the Continuous Spontaneous Localization (CSL) model of Pearle, Ghirardi and Rimini [27, 17] and the gravity related collapse model of Diósi [6], both aimed at solving the measurement problem (see [2, 3] for a review). The generalization of the formalism to the non-Markovian realm was carried out by Diósi and Strunz [30, 11]. In that case, the measurement interpretation exists but is admittedly more subtle: a given trajectory only has a local point by point measurement interpretation [15] but no real time meaning [9, 10, 36], at least within orthodox interpretations of the formalism [16]. Non-Markovian stochastic approaches have also infected foundations with various extensions of the initial CSL proposal [1, 13]. In all these applications, the practical problem is always the same: the equations are very formal, involving functional derivatives, and plotting their solution is typically as hard as solving a general non-Markovian open-system evolution. There are, of course, quite a few cases where the equations can be solved [31, 23, 37, 24, 25] as well as perturbative methods [5] to deal with a broader class of problems, but general non-perturbative approaches have so far been lacking.

In the following, we will introduce stochastic Schrödinger equations as unravelings of non-Markovian open-system dynamics making an extensive use of the general framework introduced recently by Diósi and Ferialdi [7]. Yet, our objective will not primarily be to find a numerically efficient way to compute the open-system evolution: we will be interested in the stochastic pure state trajectories themselves, having in mind their broader range of applications.

2 General framework

2.1 Gaussian master equations

We consider an open quantum system evolving according to a general Gaussian Master Equation (GME), one of the broadest non-Markovian generalizations of the Lindblad equation that is analytically tractable. We suppose the system has a proper Hamiltonian $H_{0}$ and will use the interaction picture for all operators $\hat{\mathcal{O}}(t)=e^{iH_{0}t}\hat{\mathcal{O}}e^{-iH_{0}t}$. The GME for the density matrix in interaction picture can then be written as a time ordered exponential [7]:

[TABLE]

where $\mathcal{T}$ is the time ordering operator, $\{\hat{A}^{k}\}_{1\leq k\leq n}$ are arbitrary system Hermitian operators, $D_{ij}(\tau,s)$ is a complex positive semi-definite kernel, $\theta_{\tau s}=\theta(\tau-s)$ is the Heaviside function, and we have used summation on repeated indices as well as the left-right notation for superoperators:

[TABLE]

Equation (2.1) is a direct non-Markovian generalization of the Lindblad master equation and reduces to it in the limit $D_{ij}(\tau,s)\rightarrow d_{ij}\,\delta(\tau-s)$. It is also an equivalent operator rewriting of the well known quadratic Feynman-Vernon influence functional for systems linearly coupled to harmonic oscillators [14], a representation used e.g. by Hu, Paz and Zhang [20, 21] for the Quantum Brownian Motion and which can be derived from a wide class of microscopic models [34].

At least if $D$ is time translation invariant, equation (2.1) can be obtained from the linear coupling of the system of interest with a general bosonic bath (without the Born-Markov approximation). To avoid being excessively abstract, let us recall the explicit derivation of [7]. We consider a bath constituted of a continuum of frequency of $n$ types of harmonic oscillators. We write the corresponding creation and annihilation operators $a_{k}(\omega)$, $a_{k}^{\dagger}(\omega)$ which verify $[a_{i}(\omega),a_{j}^{\dagger}(\omega^{\prime})]=\delta_{ij}\delta(\omega-\omega^{\prime})$. We then couple our system of interest with the bath through the interaction Hamiltonian:

[TABLE]

where $\kappa_{k}^{\ell}(\omega)$ is an arbitrary matrix of complex coefficients. If the initial state of the total system is a product with the bath ground state, the reduced system density matrix can be expressed exactly using Wick’s theorem [7] and is given by formula (2.1). In that case, $D$ is given by:

[TABLE]

Here, we are not interested in the open system dynamics per se and thus step back from its implementation details to take the very general equation (2.1) as our starting point.

2.2 Stochastic unraveling

A stochastic unraveling of equation (2.1) is a set of pure state trajectories $|\psi_{\xi}(t)\rangle$ indexed by a set of random fields $\xi=\{\xi_{k}\}_{1\leq k\leq n}$ such that:

[TABLE]

where $\mathds{E}_{\xi}[\cdot]=\int\cdot\,\mathrm{d}\mu_{\circ}(\xi)$ denotes averaging over the set of complex stochastic fields $\xi$. The interest of this rewriting is that it decouples the left and right parts of the time ordered exponential $\Phi_{t}$ of equation (2.1), a trick which is sometimes called a Hubbard-Stratonovich transformation. Let us now find an explicit candidate for $|\psi_{\xi}(t)\rangle$. Following e.g. [7], we first introduce a set of Gaussian complex fields of zero average and thus fully characterized by its following two point correlation functions:

[TABLE]

where $S$ is the relation function, a symmetric complex kernel which can be chosen freely provided the full correlation kernel is positive semi-definite.

For two sets of fields $a_{t}^{k}$ and $b_{t}^{k}$, the following generalized characteristic function is obtained by Gaussian integration:

[TABLE]

an equality that can be further generalized for $a$’s and $b$’s promoted to operators provided both sides are time ordered:

[TABLE]

The similarity between the r.h.s of (7) and (2.1) suggests to identify $\mathds{E}_{\xi}[|\psi_{\xi}\rangle\langle\psi_{\xi}|]$ with the l.h.s of (7) with $\hat{a}^{k}_{s}=\hat{A}_{L}^{k}(s)$, $\hat{b}^{k}_{s}=\hat{A}^{k}_{R}(s)$ and $C_{ij}(\tau,s)=D_{ij}(\tau,s)$. Actually, this does not fully do the trick, a counter term remains, and one sees following [7] that the correct prescription is:

[TABLE]

which fulfills the unraveling condition (4). Equation (2.2) can subsequently be written in differential form to yield a non-Markovian stochastic Schrödinger equation:

[TABLE]

To our knowledge, this equation is the most general linear non-Markovian stochastic Schrödinger equation ever proposed [7]. Unfortunately, equation (9) is quite formal unless it is possible to write the functional derivative111The reader puzzled by the rigorous definition of the functional derivative can as well use (2.2) in all situations. as a simple local operator acting on $|\psi_{\xi}(t)\rangle$.

3 Main result

The core objective of this article is to show how the solutions of (9) can nonetheless be computed in the general case by reusing a stochastic unraveling trick, i.e. by writing this time:

[TABLE]

where $\eta$ is an auxiliary set of classical Gaussian fields. More precisely, we can try to write the evolution operator of (2.2) in the following way:

[TABLE]

Using again the Gaussian integration formula (7), one sees that the new unraveling condition (10) is satisfied provided $\eta$ is a set of Gaussian fields with zero mean and two-point correlation functions given by:

[TABLE]

with:

[TABLE]

and where the other correlation function $J$ can be chosen freely provided the total kernel $\Gamma$:

[TABLE]

is positive semi-definite222For a kernel $K$ given by equation (13), there always exists a kernel $J$ such that $\Gamma$ is positive semi-definite. The existence of a random field $\eta$ verifying (12) is then a consequence of the Bochner-Minlos theorem..

Equation (11) can be written in an explicit linear differential form free of functional derivatives:

[TABLE]

This latter equation, together with its unraveling interpretation (10) is the main result of this article. The evolution given by equation (15) is time-local in the sense that once the random fields $\xi$ and $\eta$ are fixed, the state can be evolved without reference to the past. Although our interest was primarily in the stochastic pure state $|\psi_{\xi}\rangle$ itself, we can write the open-system density matrix as a double average over our successive unravelings (4) and (10):

[TABLE]

This can be written equivalently:

[TABLE]

where $\eta^{(1)}$ and $\eta^{(2)}$ are independent. This is the two state unraveling proposal of Stockburger and Grabert [29, 28]. We have thus found a connection between the standard stochastic pure state representations and the more exotic two state vector method used in numerical approaches to open-systems. This allows us to identify at least part of the freedom in the noise present in this latter method as coming from the kernel $S$ which encodes different stochastic Schrödinger evolutions corresponding to the same open-system evolution. We can further risk the following heuristic interpretation of the two noises. The noise $\xi$ is classical in the sense that it is averaged over at the density matrix level. The noise $\eta$, on the other hand, is quantum or coherent as all the possible contributions are summed over coherently at the pure state level.

4 Norm preserving unravelings

4.1 Context and objectives

The linear stochastic differential equation (9) does not preserve the norm of the state vector $|\psi_{\xi}\rangle$. Normalized pure states are usually preferred in foundations but also for numerics: the norm of $|\psi_{\xi}\rangle$ typically diffuses and for large times, most trajectories give a vanishing contribution to the average (4) while nearly all the weight is provided by rare events.

To obtain a norm preserving evolution from the linear one, one needs to define a normalized state $|\widetilde{\psi}_{\xi}\rangle$:

[TABLE]

and also to introduce a new “cooked” probability measure:

[TABLE]

which insures that the unraveling condition $\rho(t)=\mathds{E}^{t}_{\xi}\left[\,|\widetilde{\psi}_{\xi}(t)\rangle\langle\widetilde{\psi}_{\xi}(t)|\right]$ –where $\mathds{E}^{t}_{\xi}$ is the expectation value taken with the new measure $\mu_{t}$– is still valid. It is important to note at that stage that the state $|\widetilde{\psi}_{\xi}(t)\rangle$ with the probability measure $\mathrm{d}\mu_{t}(\xi)$ unravels the open-system evolution only for a single instant $t$. If one wants a trajectory that unravels the open-system evolution at all times, it is necessary to dynamically redefine the field variable $\xi$. It can be done by introducing $\xi^{[t]}$, a complete random field trajectory depending on $\xi$ such that for any functional $f$:

[TABLE]

This way, starting from a field realization $\xi$ drawn from the Gaussian measure $\mathrm{d}\mu_{\circ}$, one can compute a trajectory $t\mapsto|\widetilde{\psi}_{\xi^{[t]}}(t)\rangle$ that unravels the open evolution for all times. By “computing the norm preserving unraveling” one may thus want to sample two different things:

the state $|\widetilde{\psi}_{\xi^{[t]}}(t)\rangle$ for a single time, 2. 2.

the full trajectory $t\mapsto|\widetilde{\psi}_{\xi^{[t]}}(t)\rangle$.

While knowing the state at a single time may be enough in most cases, it is not sufficient to compute correlation functions of the stochastic states at different times (something one might want to do for collapse models in foundations). In any case, the ideal solution would be to find a way to write an explicit stochastic differential equation sampling the states of interest directly. This is unfortunately not a trivial endeavour and the method we shall propose is arguably less appealing than in the linear case.

4.2 Computation of the full trajectory for S=0

When $S=0$, there exists an efficient way to construct the full norm preserving trajectory $t\mapsto|\widetilde{\psi}_{\xi^{[t]}}(t)\rangle$, exploiting the fact that $\xi$ and $\xi^{[t]}$, defined implicitly in (20), can be connected explicitly in that case. Indeed, after relatively painful computations (see [19, 35, 33]) one can show that for $S=0$:

[TABLE]

with $\langle A_{\ell}(s)\rangle_{s}=\langle\widetilde{\psi}_{\xi^{[s]}}(s)|A^{\ell}(s)|\widetilde{\psi}_{\xi^{[s]}}(s)\rangle$. Using this result, it is possible to compute $|\widetilde{\psi}_{\xi^{[t]}}(t)\rangle$ inductively.

We start from $\xi^{[0]}=\xi$ sampled from the Gaussian measure $\mathrm{d}\mu_{\circ}$. At a time $v\leq t$ we assume that:

• –

$\forall u\leq v$, $|\widetilde{\psi}_{\xi^{[u]}}(u)\rangle$ is known,

• –

$\forall u\leq t$, $\xi^{[v]}_{u}$ is known.

We first compute $\xi^{[v+\mathrm{d}v]}$ for $u\leq t$ using (21):

[TABLE]

We then compute $|\widetilde{\psi}_{\xi^{[v+\mathrm{d}v]}}(v+\mathrm{d}v)\rangle$ solving the linear equation (9) up to time $v+\mathrm{d}v$ with the new complete noise trajectory $\xi^{[v+\mathrm{d}v]}$ and normalizing the result. This is done using the unraveling (15). By induction, we thus obtain $|\widetilde{\psi}_{\xi^{[t]}}(t)\rangle$ for all $t$ starting from a single realization of the noise $\xi$ (see figure 1). For that matter, it is necessary to solve the linear stochastic Schrödinger equation from the beginning with a new field at every time step. Hence this adds a linear a overhead $\propto t/\mathrm{d}t$ to the time needed to compute linear trajectories. As the errors coming from the discretization are typically much smaller than the statistical errors coming from the averaging of (15), there is however no loss in precision.

4.3 General case, single time

In the general case, the previous method cannot be used as we do not know of a generalization of formula (20) for $S\neq 0$. Finding a way to sample from norm-preserving trajectories has eluded us. In most cases however, one is only interested in the properties of the stochastic state at a single time. In such a situation, the linear unraveling is sufficient to compute everything. Indeed, let us consider an arbitrary function $\varphi$ of the norm-preserving trajectory $|\widetilde{\psi}_{\xi^{[t]}}(t)\rangle$ at a fixed time $t$. We have:

[TABLE]

That is, expectation values of all time-local functionals of the norm preserving trajectory can be computed by averaging over the solutions of the linear NMSSE.

5 Example

5.1 Linear evolution

We now illustrate our method on an example where the stochastic Schrödinger equation is fully explicit. For that matter, we consider a single Hermitian operator $\hat{A}$ commuting with $H_{0}$ (s.t. $\hat{A}(t)=\hat{A}$) and purely isotropic noise, i.e. $S=0$, which gives the following stochastic Schrödinger equation:

[TABLE]

For a given realization of $\xi$, $|\psi_{\xi}\rangle$ obeys an explicit differential equation. Indeed, using the integral form (2.2), we can replace the functional derivative:

[TABLE]

This example thus offers a nice test bed for our stochastic approach. In this example, the relaxation timescale responsible for the non-Markovianity and the timescale of the system-bath coupling are both of order $1$ and we would thus be deep in the non-Markovian regime had the system Hamiltonian not been trivial. To compute the linear trajectories, we average over the solutions of (15) with $\eta$ a real noise such that $\mathds{E}[\eta_{t}\eta_{s}]=D(t,s)$. Numerical results for $\hat{A}=\mathrm{diag}(1,0,-1)$, $D(t,0)=e^{-t}$, and $|\psi(0)\rangle=(1,1,1)/\sqrt{3}$ are shown in figure 2.

Our method is efficient for evolution times of the order of the the system-bath coupling timescale but becomes expensive for larger times. Indeed, we are summing states $|\psi_{\xi,\eta}(t)\rangle$ with a phase with respect to $|\psi_{\xi}(t)\rangle$ that is more and more random as time grows. This is an important limit of this scheme for large time.

5.2 Norm preserving evolution

We also sample the norm-preserving non-linear trajectories using the method of section 4.2. The precision for the sampling of a single trajectory $|\psi_{\xi^{[t]}}(t)\rangle$ is very similar to that of the linear case (see figure 2).

To test the robustness of the method, we also compute the probability distribution of $\langle\psi_{\xi^{[t]}}(t)|\hat{A}|\psi_{\xi^{[t]}}(t)\rangle$ for two different times using norm-preserving trajectories computed with the numerical method of 4.2 and using exact linear trajectories and the re-weighting of 4.3. Up to statistical and discretization errors, the second method is exact. The results are displayed in figure 3 and show an excellent agreement between the two methods.

6 Conclusion

In this article, we have introduced a rewriting (or “unraveling”) of the solutions of formal NMSSE as averages over the solutions of explicit time-local stochastic differential equations. This new formulation may provide a basis for further analytical studies, notably as it gives an extremely simple new way to define NMSSE. More importantly, it immediately offers a method to compute numerically their solutions. Although we have proposed a heuristic interpretation of the quantum and classical character of the two noises involved, the auxiliary state $|\psi_{\xi,\eta}\rangle$ does not yet have a clear operational characterization. The latter would certainly help understand the freedom in the noise kernel $J$. We have also extended our method to non-linear norm-preserving NMSSE in the case where $S=0$, i.e. where the complex noise $\xi$ is isotropic. Computing expectation values of time non-local functionals of norm-preserving trajectories when $S\neq 0$ is still an open problem. Our work has been focused on the NMSSE as an end in themselves rather than on the open-system evolution they unravel. It is however possible that our methods, especially the one to sample non-linear NMSSE directly, might ultimately help improve our ability to deal with non-Markovian open-systems numerically.

I thank Lajos Diósi and Luca Ferialdi for helpful discussions. Part of this work was supported by the Alexander von Humboldt foundation and the Agence Nationale de la Recherche (ANR) contract ANR-14-CE25-0003-01.