# Quantum trajectories for time-local non-Lindblad master equations

Tobias Becker<EMAIL_ADDRESS>Institut für Theoretische Physik,
Technische Universität Berlin, Hardenbergstrasse 36, 10623 Berlin, Germany
Ché Netzer  Institut für Theoretische Physik, Technische Universität Berlin,
Hardenbergstrasse 36, 10623 Berlin, Germany  André Eckardt  eckardt@tu-
berlin.de Institut für Theoretische Physik, Technische Universität Berlin,
Hardenbergstrasse 36, 10623 Berlin, Germany

###### Abstract

For the efficient simulation of open quantum systems we often use quantum jump
trajectories given by pure states that evolve stochastically to unravel the
dynamics of the underlying master equation. In the Markovian regime, when the
dynamics is described by a Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) master
equation, this procedure is known as Monte-Carlo wavefunction (MCWF) approach
. However, beyond ultraweak system-bath coupling, the dynamics of the system
is not described by an equation of GKSL type, but rather by the Redfield
equation, which can be brought into pseudo-Lindblad form. Here negative
dissipation strengths prohibit the conventional approach. To overcome this
problem, we propose a pseudo-Lindblad quantum trajectory (PLQT) unraveling. It
does not require an effective extension of the state space, like other
approaches, except for the addition of a single classical bit. We test the
PLQT for the eternal non-Markovian master equation for a single qubit and an
interacting Fermi Hubbard chain coupled to a thermal bath and discuss its
computational effort compared to solving the full master equation.

_Introduction.–_ Away from thermodynamic equilibrium the properties of an open
quantum system do not simply follow from the fundamental principles of
statistical mechanics, but depend on the very details of the surrounding
environment. This includes both transient dynamics, as the algorithm of a
quantum computer or the relaxation following a quantum quench, and non-
equilibrium steady states. Therefore, it is crucial to find an effective
equation of motion for the open system that accurately captures the impact of
the environment. At the same time, and equally importantly, the theoretical
description should allow for efficient numerical simulations. A powerful
approach for the latter are quantum trajectory simulations, where a stochastic
process for the evolution of pure states is considered, the ensemble average
of which describes the open system. Compared to the evolution of the full
density operator (scaling quadratically with the state-space dimension $D$),
these simulations require less memory, since pure states scale only linearly
with $D$. Moreover, such unravelings can also directly describe stochastic
processes of measured systems [1, 2, 3].

Quantum trajectory simulations are rather straightforward in the ultraweak-
coupling limit, where the system-bath coupling is weak compared to the
(quasi)energy level splitting in the system. In this case, the system is
described by a master equation of GKSL (Gorini-Kossakowski-Sudarshan-Lindblad)
form [4, 5] ($\hbar=1$),

$\displaystyle\dot{\varrho}=-\mathrm{i}[H,\varrho]+\sum\limits_{i}\gamma_{i}\left(L_{i}\varrho
L_{i}^{\dagger}-\frac{1}{2}\\{L_{i}^{\dagger}L_{i},\varrho\\}\right),$ (1)

with the coherent evolution captured by some Hamiltonian $H$ and dissipation
described by jump operators $L_{i}$ with associated non-negative strengths
$\gamma_{i}$. Here, $H$, $\gamma_{i}$, and $L_{i}$ can be time dependent.

From this equation, we can immediately obtain a stochastic process for the
evolution of pure states, known as Monte-Carlo wavefunction (MCWF) approach
[6, 7, 8, 9, 10, 11, 12, 13]. In each time step $\delta t$, the state either
evolves coherently according to $\ket{\psi(t+\delta
t)}\propto(1-\mathrm{i}\delta tH_{\mathrm{eff}}(t))\ket{\psi(t)}$ with
probability $1-\sum_{i}r_{i}(t)\delta t$ and effective Hamiltonian

$\displaystyle
H_{\mathrm{eff}}(t)=H-\frac{\mathrm{i}}{2}\sum_{i}\gamma_{i}L_{i}^{\dagger}L_{i},$
(2)

or a quantum jump occurs, $\ket{\psi(t+\delta t)}\propto L_{i}\ket{\psi(t)}$,
with probability $r_{i}(t)\delta t$, with jump rates
$r_{i}(t)=\gamma_{i}\matrixelement{\psi(t)}{L_{i}^{\dagger}L_{i}}{\psi(t)}$.
The state of the system is then given (or approximated) by the ensemble
average $\rho(t)=\overline{\outerproduct{\psi(t)}{\psi(t)}}$ over an
infinitely (or sufficiently) large number $N$ of trajectories
$\ket{\psi_{n}(t)}$, where $\overline{X}\equiv\frac{1}{N}\sum_{n=1}^{N}X_{n}$.

However, the assumption of ultraweak coupling is questionable in various
situations, for instance in large systems, with small finite-size gaps and
tiny avoided crossings between many-body states, as well as in Floquet systems
with driving frequency $\omega$, where the average quasi energy level spacing
is $\omega/D$ [14].

Beyond ultraweak coupling, master equations in pseudo-Lindblad form can be
found, which look like a GKSL master equation Eq. 1, except for the fact that
the coefficients $\gamma_{i}$ also take negative values. For instance, the
Redfield equation obtained in (Floquet)-Born-Markov approximation can be
brought to this form [15]. Generally, negative relaxation strengths are
relevant for non-Markovian dynamics [16], stochastic Hamiltonians with non-
Markovian noise [17], gauge transformed Lindbladians [18] and exact master
equations [19, 20, 21]. These negative values are incompatible with the
conventional MCWF, since the probability $r_{i}(t)\delta t$ for a quantum jump
would become negative. To overcome this problem different quantum jump
unravelings have been formulated, which, however, require significantly more
computational resources [22, 23, 24, 25, 26, 27, 28, 29, 30]. In many
approaches the system’s state space needs to be extended, so that its
dimensionality at least doubles [22, 23, 24, 31, 25, 26, 32, 33]. For
oscillating strengths between positive and negative values, moreover, an
alternative non-Markovian quantum jump method (NMQJ) has been proposed in
which jump processes are inverted [34, 27, 28, 29]. This method does not work,
if $\gamma_{i}<0$ for all times and does not admit independent (i.e. parallel)
evaluation of trajectories. For non-oscillatory strengths the rate operator
quantum jump approach can be applied [30], however, it requires a rather
costly diagonalization of a state-dependent operator in every time step of the
evolution. A generalization of the non-Markovian jump method to many-body
systems has been proposed in Ref. [35] and can be used to study measurement
induced phase transitions.

In this work we propose pseudo-Lindblad quantum trajectories (PLQT), which
work for arbitrary $\gamma_{i}$, where the trajectories evolve independently
and which does not require the doubling of state space. In the following, this
is realized by extending the system’s state space in a minimal (and for the
memory requirement of simulations practically irrelevant) fashion by a single
classical bit $s\in\\{-1,+1\\}$, $\ket{\psi(t)}\to\\{\ket{\psi(t)},s(t)\\}$.

_Algorithm.–_ To unravel the dynamics of a pseudo-Lindblad quantum master
equation by quantum trajectories $\\{\ket{\psi(t)},s(t)\\}$, first choose a
time step $\delta t$, which is sufficiently short for the first-order time
integration, and jump rates $r_{i}(t)>0$ for each jump operator $L_{i}$ (to be
specified below). Within one time step a quantum jump occurs described by

$\displaystyle\begin{aligned} \ket*{\psi^{(i)}(t+\delta
t)}&=\frac{\sqrt{|\gamma_{i}|}L_{i}\ket{\psi(t)}}{\sqrt{r_{i}(t)}},\\\
s^{(i)}(t+\delta t)&=\frac{\gamma_{i}}{|\gamma_{i}|}\,s(t),\end{aligned}$ (3)

with probability $r_{i}(t)\delta t$ or, alternatively, with the remaining
probability $1-\sum_{i}r_{i}(t)\delta t$ the state evolves coherently, with
$H_{\mathrm{eff}}$ [Eq. 2] [36]

$\displaystyle\ket*{\psi^{(0)}(t+\delta t)}$
$\displaystyle=\frac{(1-\mathrm{i}\delta
tH_{\mathrm{eff}}(t))\ket{\psi(t)}}{\sqrt{1-\delta t\sum_{i}r_{i}(t)}},$ (4)
$\displaystyle s^{(0)}(t+\delta t)$ $\displaystyle=s(t).$ (5)

We now show that the ensemble of pure states obtained by this PLQT approach,
converges to the correct density operator solving Eq. 1,

$\displaystyle\varrho(t)=\overline{s(t)\outerproduct{\psi(t)}{\psi(t)}}.$ (6)

For a pure initial state, $\sigma(t)=s(t)\outerproduct{\psi(t)}{\psi(t)}$, on
average the update scheme is the weighted sum of the processes described
above:

$\displaystyle\begin{aligned} \overline{\sigma(t+\delta
t)}=&\sum_{i}r_{i}(t)\delta t\ \sigma^{(i)}(t+\delta t)\\\
&+\bigg{(}1-\sum_{i}r_{i}(t)\delta t\bigg{)}\sigma^{(0)}(t+\delta
t),\end{aligned}$ (7)

with $\sigma^{(i)}=s^{(i)}\outerproduct{\psi^{(i)}}{\psi^{(i)}}$ and
$\sigma^{(0)}=s^{(0)}\outerproduct{\psi^{(0)}}{\psi^{(0)}}$. By inserting Eqs.
3 and 4, the jump rates $r_{i}(t)$ cancel out and one arrives at

$\displaystyle\begin{split}\overline{\sigma(t+\delta t)}=&\sigma(t)+\delta
t\bigg{(}\sum_{i}\gamma_{i}L_{i}\sigma(t)L_{i}^{\dagger}\\\
&-\mathrm{i}H_{\mathrm{eff}}(t)\sigma(t)+\mathrm{i}\sigma(t)H_{\mathrm{eff}}(t)^{\dagger}\bigg{)},\end{split}$
(8)

almost corresponding to the action of the master Eq. 1. The final step to
arrive at Eq. 1 is to average Eq. 8 also over an ensemble of pure states at
time $t$, so that $\overline{\sigma(t+\delta)}\to\varrho(t+\delta t)$ and
$\sigma(t)\to\varrho(t)$. As will be discussed below, one consequence of the
presence of negative weights $\gamma_{i}<0$ is that individual wave functions
$\psi_{n}$ are not normalized. As a result, the ensemble averaged trace is
preserved only in the limit $N\to\infty$ of an infinite ensemble 111The
fluctuations of the trace and their dependence on $N$ are shown in Fig. 1 of
the supplemental material [38] for an explicit example.. Therefore, in a
finite ensemble, one obtains better convergence by explicit normalization,
$\varrho_{N}=\frac{1}{\mathcal{N}}\sum_{n}^{N}s_{n}\outerproduct{\psi_{n}}{\psi_{n}}$,
with $\mathcal{N}=\sum_{n}^{N}s_{n}\innerproduct{\psi_{n}}{\psi_{n}}$ at every
time $t$. A rigorous proof of our algorithm using the Ito-formalism is
outlined in the supplemental material [38]. In case that all $\gamma_{i}$ are
positive, the sign bits do not change and the algorithm corresponds to the
conventional MCWF approach.

Note that, recently, also another unraveling of non-Lindblad master equations
was proposed in Ref. [39]. It is different from our approach, but also
involves an effective classical degree of freedom, given by a real number of
constant average, rather than our single bit, whose average is time-dependent,
as will be seen below.

For the PLQT approach, as for other unraveling schemes [23], the jump rates
$r_{i}(t)>0$ can, in principle, be chosen arbitrarily. In practice there is,
however, a trade off. Whereas for too small rates $r_{i}$, large ensembles of
trajectories are required to sample each jump process $i$ sufficiently, we
also have to require that the probability $1-\sum_{i}r_{i}\delta t$ remains
positive and large enough for the given time step $\delta t$. A typical choice
is [13]

$\displaystyle r_{i}(t)=|\gamma_{i}|\frac{\lVert
L_{i}\ket{\psi(t)}\rVert^{2}}{\lVert\ket{\psi(t)}\rVert^{2}},$ (9)

for which the quantum jump does not alter the norm
$\lVert\ket{\psi}\rVert\equiv\innerproduct{\psi}{\psi}^{1/2}$ of the state,
i.e. $\lVert\ket{\psi^{(i)}(t+\delta t)}\rVert=\lVert\ket{\psi(t+\delta
t)}\rVert$. Note, however, that for $\gamma_{i}<0$ this choice implies that
the norm increases during the coherent time evolution with $H_{\mathrm{eff}}$,
$\lVert\ket{\psi^{(0)}(t+\delta t)}\rVert=(1+\delta
t\sum_{\gamma_{i}<0}r_{i}(t))\lVert\ket{\psi(t)}\rVert$ [40]. This is not the
case for the conventional MCWF approach, where $\gamma_{i}\geq 0$.

Non-Markovian dephasing for single qubit.– As a proof of principle we
implement the PLQT algorithm for a qubit subjected to purely dissipative
dynamics,

Figure 1: Non-Markovian dynamics [Eq. 10] for density matrix elements
$\varrho_{00}$ (a), $\mathrm{Re}\varrho_{01}$ (solid),
$\mathrm{Im}\varrho_{01}$ (dashed) (b) (and the Bloch-vector in the
$x$-$y$-plane (d). Analytical solution (black) and unraveling with $N=10^{5}$
PLQTs with time step $\delta t=0.01$ in blue ($N=10,10^{3}$ in thin and
intermediate grey lines) for an initial Bloch state with $\phi=\Theta=\pi/4$.
(c) Shows the averaged sign bit.

$\displaystyle\begin{aligned}
\dot{\varrho}(t)&=\frac{1}{2}\Big{[}\mathcal{L}_{x}+\mathcal{L}_{y}-\tanh(t)\mathcal{L}_{z}\Big{]}\varrho(t),\end{aligned}$
(10)

with GKSL channels
$\mathcal{L}_{i}\varrho=\sigma_{i}\varrho\sigma_{i}-\varrho$, where
$\sigma_{i}$ are Pauli operators, with
$\sigma_{i}^{\dagger}\sigma_{i}=\sigma_{i}^{2}=1$. This equation is known as
the eternal non-Markovian master equation [41, 16]. The existence of a
negative relaxation rate makes it inaccessible to the standard MCWF, while
also the NMQJ approach fails, since $-\tanh(t)<0$ for all $t$.

However, for this model the PLQT approach is easily implemented and, because
the jump operators are unitary the jump rates are state independent, i.e.
$\lVert\sigma_{i}\ket{\psi(t)}\rVert^{2}=\lVert\ket{\psi(t)}\rVert^{2}$ leads
to a simplification in Eq. 9, and one has $r_{x}=r_{y}=1/2$,
$r_{z}(t)=\tanh(t)/2$. Also the effective Hamiltonian
$H_{\mathrm{eff}}=-\frac{\mathrm{i}}{2}(1-\tanh(t)/2)$ entering Eq. 4 is not
state-dependent (since it is proportional to the identity).

On average the sign follows the rate equation
$\dot{\overline{s(t)}}=-2r_{z}(t)\overline{s(t)}$, which is solved by
$\overline{s(t)}=1/\cosh(t)$, as shown in Fig. 1 (c) [38]. Quantum jumps are
realized in the Bloch-vector representation by reflections at the
$y$-$z$-plane and $x$-$z$-plane for $\sigma_{x}$ and $\sigma_{y}$,
respectively [Fig. 1 (d)]. The $\sigma_{z}$ quantum jump is a reflection at
the $x$-$y$-plane and, due to the negative relaxation strength, the sign flip
is accounted by an additional point reflection at the origin.

By simulating $N=10^{5}$ trajectories in Fig. 1 (a), (b), we obtain accurate
results for the transient dynamics until, at $t_{\mathrm{R}}\sim 2$, the
system reaches the steady state regime. Besides this physical relaxation time,
we also find an _algorithmic relaxation time_ $t_{A}\sim 4$, at which the
number of negative and positive trajectories become equal and the averaged
sign decays to zero [Fig. 1 (c)]. Beyond this algorithmic relaxation time,
fluctuations are typically increased [Fig. 1 (a) and (b)]. This effect can be
understood by noting that a stochastic process of a real variable $x_{n}$ with
positive mean $\overline{x}$ will have bounded fluctuations $\Delta
x=\overline{(x-\overline{x})^{2}}^{-1/2}\leq\overline{x}$ as long as
$x_{n}>0$, whereas $\Delta x$ is not bounded, when $x_{n}$ can also take
negative values. Thus, ideally, $t_{A}$ should be large compared to the time
span of the simulation (which is $t_{R}$, if we are interested in computing
the steady state). The algorithmic relaxation time is determined by the
inverse sign-flip rate $r_{\text{sf}}=\sum_{i,\gamma_{i}<0}r_{i}$, e.g.
$t_{A}=r_{\text{neg}}^{-1}$ for time-independent $r_{\text{sf}}$. Thus, we can
increase $t_{A}$ simply by lowering the rates for negative processes with
“rates” $\gamma_{i}<0$ relative to positive ones with $\gamma_{i}>0$. However,
this will also increase the number of trajectories needed for properly
sampling those negative-“rate” processes. Thus, before doing this, one should
first attempt to rewrite the master equation, so that the relative weight of
negative processes is reduced. This can be done for pseudo-Lindblad equations
derived from the Redfield equation [15], as we will recapitulate now.

Redfield dynamics.– For a microscopic model a master equation is often derived
within the Born-Markov-Redfield formalism [42, 43]. We consider a system-bath
Hamiltonian of the form $H_{\mathrm{tot}}=H+\sum_{i}S_{i}\otimes B_{i}+H_{i}$
with system Hamiltonian $H$ that couples to individual baths $H_{i}$ where
$S_{i}$ and $B_{i}$ denote the system and bath coupling operators,
respectively. The Redfield equation can then be written in pseudo-Lindblad
form [15]

$\displaystyle\begin{aligned}
\dot{\varrho}&=-\mathrm{i}[H+H^{\mathrm{LS}},\varrho]\\\
&+\sum\limits_{i,\sigma=\pm}\sigma\left(L_{i\sigma}\varrho
L_{i\sigma}^{\dagger}-\frac{1}{2}\\{L_{i\sigma}^{\dagger}L_{i\sigma},\varrho\\}\right),\end{aligned}$
(11)

with Lamb-shift Hamiltonian
$H^{\mathrm{LS}}=(1/2\mathrm{i})\sum_{i}S_{i}\mathbb{S}_{i}+\text{H.c.}$,
convolution operators
$\mathbb{S}_{i}=\int_{0}^{\infty}d\tau\expectationvalue{B(\tau)B}e^{\mathrm{i}H\tau}S_{i}e^{-\mathrm{i}H\tau}$
and Lindblad-like jump operators

$\displaystyle
L_{i\sigma}=\frac{1}{\sqrt{2}}\Bigg{[}\lambda_{i}(t)S_{i}+\sigma\frac{1}{\lambda_{i}(t)}\mathbb{S}_{i}\Bigg{]}$
(12)

with arbitrary, time-dependent real parameters $\lambda_{i}(t)$. We see that
due to the negative relaxation rates with $\sigma=-1$, the Redfield equation
is generally not of GKSL form, unless further approximations are employed in
the limit of ultra weak coupling [44, 43, 42] or for high bath temperatures
[45, 15]. For a purely Ohmic bath, the choice [15]

$\displaystyle\lambda_{i,\mathrm{glob}}(t)^{2}=\sqrt{\frac{\tr\
\mathbb{S}_{i}^{\dagger}\mathbb{S}_{i}}{\tr\ S_{i}S_{i}}}$ (13)

minimizes the norm of the negative channels in the pseudo-Lindblad equation
globally, i.e. on average for all states. A further reduction of negative
processes can be achieved by state-dependent optimization. Namely, according
to Eq. 9, where (assuming, without loss of generality a normalized state) the
rates for negative quantum jumps with $L_{i-}$ are given by
$r_{i-}(t;\lambda_{i}(t))=\frac{1}{2}\Big{(}\lambda_{i}(t)^{2}\lVert
S_{i}\ket{\psi(t)}\rVert^{2}+\frac{1}{\lambda_{i}(t)^{2}}\lVert\mathbb{S}_{i}\ket{\psi(t)}\rVert^{2}-2\mathrm{Re}\expectationvalue{S_{i}\mathbb{S}_{i}}{\psi(t)}\Big{)}$.
Thus, the choice,

$\displaystyle\begin{aligned}
\lambda_{i,\mathrm{loc}}(t)^{2}=\frac{\lVert\mathbb{S}_{i}\ket{\psi(t)}\rVert}{\lVert
S_{i}\ket{\psi(t)}\rVert},\end{aligned}$ (14)

minimizes the rates for negative quantum jumps in the unraveling of the
Redfield equation. Since the states in the numerator and the denominator of
Eq. 14 have to be evaluated for evolving the state anyway, this local
optimization (which is not described in Ref. [15]) can be implemented
efficiently.

Figure 2: Dynamics of scaled interaction energy of extended Hubbard chain of 2
[8] particles on 4 [13] sites (a) [(b)] with $V/J=7$. We compare the dynamics
of the isolated (grey line) and open system with $\gamma/J=0.02$ and $T/J=1$
(blue line for PLQT, thin orange for Redfield equation). The decrease of
interactions is related to bath-induced doublon-breaking processes.

We test our method, using state-dependent minimization Eq. 14, for the
extended Hubbard chain of spinless fermions,

$\displaystyle
H=-J\sum\limits_{\ell=0}^{M-1}\Big{(}a_{\ell}^{\dagger}a_{\ell+1}+a_{\ell+1}^{\dagger}a_{\ell}\Big{)}+V\sum\limits_{\ell=0}^{M-1}a_{\ell}^{\dagger}a_{\ell}a_{\ell+1}^{\dagger}a_{\ell+1},$
(15)

with fermionic operators $a_{\ell}$, tunneling strength $J$ and nearest-
neighbour interaction strength $V$. For the dissipator, we have

$\displaystyle\matrixelement{n}{\mathbb{S}_{\ell}}{m}=\frac{J(\Delta_{nm})}{e^{\Delta_{nm}/T}-1}\matrixelement{n}{S_{\ell}}{m},$
(16)

with system operator $S_{\ell}=a_{\ell}^{\dagger}a_{\ell}$, level splitting
$\Delta_{nm}=E_{n}-E_{m}$ and bath temperature $T$. We consider a purely Ohmic
bath, with spectral density $J(E)=\gamma E$ and coupling strength $\gamma$.

In Fig. 2 we depict the decay of the interaction energy for an initial state
in which pairs of adjacent sites are occupied
$\ket{\psi(0)}=\ket{011011\dots}$. Quench dynamics for such a charge density
wave in a spin polarized Fermi Hubbard model have been recently observed
experimentally in the Bakr group [46]. We assume strong interactions, $V/J=7$,
for which the dublon pairs can only be broken when the system exchanges energy
with the bath. This leads to a decay of the energy of the open system as
depicted in Fig. 2 (a), where the transient oscillations are well reproduced.

Numerical implementation.– Let us now discuss the numerical implementation of
the PLQT approach. Since the trajectories are independent, we run them in
parallel. Depending on the physical quantity of interest, let’s say observable
$A$, it is often reasonable not to store the actual time dependent state as
large vectors with complex entries, but rather expectation values
$\expectationvalue{A}{\psi_{n}(t)}$ together with the norm
$\lVert\psi_{n}(t)\rVert$ and the sign $s_{n}(t)$. While the storage of the
trajectory data boils down to a few real numbers, the time evolution requires
the full state vector. The memory needed for the time integration of a quantum
trajectory would grow linear with the state-space dimension $D$, if not only
the Hamiltonian, but also the jump operators were sparse. The latter is the
case, however, mainly in phenomenological master equations with local jump
operators and not for the Redfield equation, so that the memory needed usually
scales like $D^{2}$. The memory needed for integrating the Redfield equation
scales equally like $D^{2}$ (since it is sufficient to store and apply the
jump operators rather than the full superoperator). Nevertheless, we find that
the memory requirement for quantum trajectories to be much lower than that for
integrating the master equation. In Fig. 3 the required memory, (a) and (b),
and single-time-step run time, (c) and (d), is compared for solving the full
Redfield master equation (blue) and a single trajectory (red).

Figure 3: Required memory and single-time-step run time for both a single
trajectory (blue) and the Redfield master equation (red). The data was
obtained for an Intel Core i9-10900 processor with up to $5.2\,\mathrm{GHz}$.

We find that the required memory is noticeably reduced for the quantum
trajectory simulation, even though, as discussed above, it still scales like
$D^{2}$. (The latter is not specific to our approach, but generically the case
also for other forms of quantum trajectory simulations). For the run time the
relative reduction is even stronger and shows different scaling with $D$.
Essentially the difference are two matrix-matrix products needed for the
Redfield integration and one matrix-vector product for the PLQTs. Note that
the unravelling can also be combined with matrix product states (e.g. [47,
48]). It is interesting to see in how far such an approach would compare to a
representation of the density operator by matrix product operators (e.g. Refs.
[49, 50]).

Finally, we would like to mention that while finalizing this manuscript, a new
approach for microscopically deriving accurate local quantum master equations
was proposed [51]. It is a very promising perspective to apply the the PLQT to
this equation, which can be broad into pseudo-Lindblad form [51]. Namely, it
will give rise to equations of motion involving spars matrices for Hamiltonian
and jump operators only. As a result, the overall memory requirement will
scale linearly with the Hilbert-space dimension $D$ only. Moreover, this could
allow for efficient matrix-product state approaches for individual
trajectories.

Conclusion.– We have developed an efficient unraveling of master equations of
pseudo-Lindblad form, which includes the Redfield equation as an important
case [15]. Different from previous approaches, it requires a minimal extension
of the state space by one classical sign bit only, is applicable also for
dissipation strengths that are always negative, does not require matrix
diagonalization during the time integration, and allows for a parallel
implementation, since all trajectories are independent from each other. We
believe that it will be a useful tool for the simulation of open many-body
systems beyond ultra weak system-bath coupling. In future work it will be
interesting to systematically investigate the impact of negative dissipation
strengths $\gamma_{i}$ on the required ensemble size and the optimal choice of
the corresponding rates for efficient simulation. Moreover, our algorithm
should be compared to the influence-martingale approach [39]. Finally, the
combination of PLQT simulations with the recently derived local quantum master
equation [51] opens a promising route for the efficient simulation of open
many-body quantum systems.

###### Acknowledgements.

We thank Francesco Intravaia, Charlotte Maurer and Alexander Schnell for
fruitful discussions. This research was funded by by the Deutsche
Forschungsgemeinschaft (DFG) via the Research Unit FOR 2414 under project No.
277974659.

## References

  * Basché  _et al._ [1995] T. Basché, S. Kummer, and C. Bräuchle, Direct spectroscopic observation of quantum jumps of a single molecule, Nature 373, 132 (1995).
  * Gleyzes  _et al._ [2007] S. Gleyzes, S. Kuhr, C. Guerlin, J. Bernu, S. Deléglise, U. Busk Hoff, M. Brune, J.-M. Raimond, and S. Haroche, Quantum jumps of light recording the birth and death of a photon in a cavity, Nature 446, 297 (2007).
  * Minev  _et al._ [2019] Z. K. Minev, S. O. Mundhada, S. Shankar, P. Reinhold, R. Gutiérrez-Jáuregui, R. J. Schoelkopf, M. Mirrahimi, H. J. Carmichael, and M. H. Devoret, To catch and reverse a quantum jump mid-flight, Nature 570, 200 (2019).
  * Gorini  _et al._ [1976] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Completely positive dynamical semigroups of n‐level systems, J. Math. Phys. 17, 821 (1976).
  * Lindblad [1976] G. Lindblad, On the generators of quantum dynamical semigroups, Comm. Math. Phys. 48, 119 (1976).
  * Dalibard  _et al._ [1992] J. Dalibard, Y. Castin, and K. Mølmer, Wave-function approach to dissipative processes in quantum optics, Phys. Rev. Lett. 68, 580 (1992).
  * Mølmer  _et al._ [1993] K. Mølmer, Y. Castin, and J. Dalibard, Monte carlo wave-function method in quantum optics, J. Opt. Soc. Am. B 10, 524 (1993).
  * Gardiner  _et al._ [1992] C. W. Gardiner, A. S. Parkins, and P. Zoller, Wave-function quantum stochastic differential equations and quantum-jump simulation methods, Phys. Rev. A 46, 4363 (1992).
  * Dum  _et al._ [1992] R. Dum, A. S. Parkins, P. Zoller, and C. W. Gardiner, Monte carlo simulation of master equations in quantum optics for vacuum, thermal, and squeezed reservoirs, Phys. Rev. A 46, 4382 (1992).
  * Gisin and Percival [1992a] N. Gisin and I. C. Percival, Wave-function approach to dissipative processes: are there quantum jumps?, Phys. Lett. A 167, 315 (1992a).
  * Gisin and Percival [1992b] N. Gisin and I. C. Percival, The quantum-state diffusion model applied to open systems, J. Phys. A 25, 5677 (1992b).
  * Carmichael [1999] H. J. Carmichael,  _Statistical Methods in Quantum Optics 1_ (Springer-Verlag, Berlin, 1999).
  * Daley [2014] A. J. Daley, Quantum trajectories and open many-body quantum systems, Adv. Phys. 63, 77 (2014).
  * Eckardt [2017] A. Eckardt, Colloquium: Atomic quantum gases in periodically driven optical lattices, Rev. Mod. Phys. 89, 011004 (2017).
  * Becker  _et al._ [2021] T. Becker, L.-N. Wu, and A. Eckardt, Lindbladian approximation beyond ultraweak coupling, Phys. Rev. E 104, 014110 (2021).
  * Megier  _et al._ [2017] N. Megier, D. Chruściński, J. Piilo, and W. T. Strunz, Eternal non-markovianity: from random unitary to markov chain realisations, Sci. Rep. 7, 6379 (2017).
  * Groszkowski  _et al._ [2023] P. Groszkowski, A. Seif, J. Koch, and A. A. Clerk, Simple master equations for describing driven systems subject to classical non-Markovian noise, Quantum 7, 972 (2023).
  * McDonald and Clerk [2023] A. McDonald and A. A. Clerk, Third quantization of open quantum systems: new dissipative symmetries and connections to phase-space and keldysh field theory formulations,  (2023), arXiv:2302.14047 .
  * Karrlein and Grabert [1997] R. Karrlein and H. Grabert, Exact time evolution and master equations for the damped harmonic oscillator, Phys. Rev. E 55, 153 (1997).
  * Shaji and Sudarshan [2005] A. Shaji and E. Sudarshan, Who’s afraid of not completely positive maps?, Phys. Lett. A 341, 48 (2005).
  * Alipour  _et al._ [2020] S. Alipour, A. T. Rezakhani, A. P. Babu, K. Mølmer, M. Möttönen, and T. Ala-Nissila, Correlation-picture approach to open-quantum-system dynamics, Phys. Rev. X 10, 041024 (2020).
  * Breuer  _et al._ [1999] H.-P. Breuer, B. Kappler, and F. Petruccione, Stochastic wave-function method for non-markovian quantum master equations, Phys. Rev. A 59, 1633 (1999).
  * Kondov  _et al._ [2003] I. Kondov, U. Kleinekathöfer, and M. Schreiber, Stochastic unraveling of redfield master equations and its application to electron transfer problems, J. Chem. Phys. 119, 6635 (2003).
  * Breuer [2004] H.-P. Breuer, Genuine quantum trajectories for non-markovian processes, Phys. Rev. A 70, 012106 (2004).
  * Suess  _et al._ [2014] D. Suess, A. Eisfeld, and W. T. Strunz, Hierarchy of stochastic pure states for open quantum system dynamics, Phys. Rev. Lett. 113, 150403 (2014).
  * Hush  _et al._ [2015] M. R. Hush, I. Lesanovsky, and J. P. Garrahan, Generic map from non-lindblad to lindblad master equations, Phys. Rev. A 91, 032113 (2015).
  * Piilo  _et al._ [2009] J. Piilo, K. Härkönen, S. Maniscalco, and K.-A. Suominen, Open system dynamics with non-markovian quantum jumps, Phys. Rev. A 79, 062112 (2009).
  * Härkönen [2010] K. Härkönen, Jump probabilities in the non-markovian quantum jump method, J. Phys. A 43, 065302 (2010).
  * Luoma  _et al._ [2020] K. Luoma, W. T. Strunz, and J. Piilo, Diffusive limit of non-markovian quantum jumps, Phys. Rev. Lett. 125, 150403 (2020).
  * Smirne  _et al._ [2020] A. Smirne, M. Caiaffa, and J. Piilo, Rate operator unraveling for open quantum system dynamics, Phys. Rev. Lett. 124, 190402 (2020).
  * Chin  _et al._ [2010] A. W. Chin, Á. Rivas, S. F. Huelga, and M. B. Plenio, Exact mapping between system-reservoir quantum models and semi-infinite discrete chains using orthogonal polynomials, J. Math. Phys. 51, 092109 (2010).
  * Strasberg  _et al._ [2016] P. Strasberg, G. Schaller, N. Lambert, and T. Brandes, Nonequilibrium thermodynamics in the strong coupling and non-Markovian regime based on a reaction coordinate mapping, New J. Phys. 18, 073007 (2016).
  * Landi  _et al._ [2022] G. T. Landi, D. Poletti, and G. Schaller, Nonequilibrium boundary-driven quantum systems: Models, methods, and properties, Rev. Mod. Phys. 94, 045006 (2022).
  * Piilo  _et al._ [2008] J. Piilo, S. Maniscalco, K. Härkönen, and K.-A. Suominen, Non-markovian quantum jumps, Phys. Rev. Lett. 100, 180402 (2008).
  * Chiriacò  _et al._ [2023] G. Chiriacò, M. Tsitsishvili, D. Poletti, R. Fazio, and M. Dalmonte, Diagrammatic method for many-body non-markovian dynamics: Memory effects and entanglement transitions, Phys. Rev. B 108, 075151 (2023).
  * [36] By Taylor expanding the denominator in Eq. [4] to first order in the small time step $\delta t$, $1/\sqrt{1-\delta t\sum_{i}r_{i}(t)}\simeq\delta t/2\sum_{i}r_{i}(t)$ the deterministic evolution can equivalently be formulated as first order differential equation $\partial_{t}\ket*{\psi^{(0)}(t)}=\left(-\mathrm{i}H_{\mathrm{eff}}(t)+1/2\sum_{i}r_{i}(t)\right)\ket*{\psi^{(0)}(t)}$.
  * Note [1] The fluctuations of the trace and their dependence on $N$ are shown in Fig. 1 of the supplemental material [38] for an explicit example.
  * [38] _In the supplemental material, we provide a proof of the PLQT algorithm using the Ito formalism, discuss the impact of the choice of rates on the norm of the states, and derive an equation of motion for the averaged sign bit._
  * Donvil and Muratore-Ginanneschi [2022] B. Donvil and P. Muratore-Ginanneschi, Quantum trajectory framework for general time-local master equations, Nat. Commun. 13, 4140 (2022).
  * [40] The norm after the coherent time step follows from, $\innerproduct*{\psi^{(0)}(t)}{\psi^{(0)}(t)}=\bra{\psi(t)}\frac{1-\mathrm{i}\delta t(H_{\mathrm{eff}}(t)-H_{\mathrm{eff}}^{\dagger})}{1-\delta t\sum_{i}r_{i}(t)}\ket{\psi(t)}=\frac{1-\delta t\sum_{i}r_{i}(t)\mathrm{sgn}(\gamma_{i})}{1-\delta t\sum_{i}r_{i}(t)}\innerproduct{\psi(t)}{\psi(t)}$ and Taylor expansion to first order in the small time step $\delta t$.
  * Hall  _et al._ [2014] M. J. W. Hall, J. D. Cresser, L. Li, and E. Andersson, Canonical form of master equations and characterization of non-markovianity, Phys. Rev. A 89, 042120 (2014).
  * Hartmann and Strunz [2020] R. Hartmann and W. T. Strunz, Accuracy assessment of perturbative master equations: Embracing nonpositivity, Phys. Rev. A 101, 012103 (2020).
  * Tupkary  _et al._ [2022] D. Tupkary, A. Dhar, M. Kulkarni, and A. Purkayastha, Fundamental limitations in lindblad descriptions of systems weakly coupled to baths, Phys. Rev. A 105, 032208 (2022).
  * Breuer and Petruccione [2007] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, New York, 2007).
  * Wichterich  _et al._ [2007] H. Wichterich, M. J. Henrich, H.-P. Breuer, J. Gemmer, and M. Michel, Modeling heat transport through completely positive maps, Phys. Rev. E 76, 031115 (2007).
  * Guardado-Sanchez  _et al._ [2021] E. Guardado-Sanchez, B. M. Spar, P. Schauss, R. Belyansky, J. T. Young, P. Bienias, A. V. Gorshkov, T. Iadecola, and W. S. Bakr, Quench dynamics of a fermi gas with strong nonlocal interactions, Phys. Rev. X 11, 021036 (2021).
  * Verstraete  _et al._ [2008] F. Verstraete, V. Murg, and J. Cirac, Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems, Adv. Phys. 57, 143 (2008).
  * Schollwöck [2011] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Ann. Phys. 326, 96 (2011).
  * Xu  _et al._ [2019] X. Xu, J. Thingna, C. Guo, and D. Poletti, Many-body open quantum systems beyond lindblad master equations, Phys. Rev. A 99, 012106 (2019).
  * Weimer  _et al._ [2021] H. Weimer, A. Kshetrimayum, and R. Orús, Simulation methods for open quantum many-body systems, Rev. Mod. Phys. 93, 015008 (2021).
  * Schnell [2023] A. Schnell, Global becomes local: Efficient many-body dynamics for global master equations (2023), arxiv:2309.07105 .

Supplemental Material for
“Quantum trajectories for time-local non-Lindblad master equations”

##  I Pseudo-Lindblad quantum trajectory unraveling in Ito-formalism

The starting point is the pseudo-Lindblad master equation

$\displaystyle\dot{\varrho}=-i[H,\varrho]+\sum\limits_{i}\gamma_{i}\left(L_{i}\varrho
L_{i}^{\dagger}-\frac{1}{2}\\{L_{i}^{\dagger}L_{i},\varrho\\}\right),$ (A1)

with jump operators $L_{i}$ and relaxation strength $\gamma_{i}$, which we
explicitly allow to be negative. For convenience, in the following we absorb
the magnitude of relaxation strength into the jump operator
$L_{i}\propto\sqrt{|\gamma_{i}|}$ and denote the sign by
$\sigma_{i}=\gamma_{i}/|\gamma_{i}|=\pm 1$,

$\displaystyle\dot{\varrho}=-i[H,\varrho]+\sum\limits_{i}\sigma_{i}\left(L_{i}\varrho
L_{i}^{\dagger}-\frac{1}{2}\\{L_{i}^{\dagger}L_{i},\varrho\\}\right).$ (A2)

Each channel $i$ gives rise to a stochastic process, which we can unravel with
pseudo-Lindblad quantum trajectories (PLQT). In the Ito-formalism the change
of the trajectory is described as follows:

$\displaystyle\ket{\mathrm{d}\psi}=\sum_{i}\mathrm{d}N_{i}\Big{(}\frac{\sigma_{i}L_{i}}{\sqrt{r_{i}(t)}}\ket{\psi}-\ket{\psi}\Big{)}+\mathrm{d}t\Big{(}\sum_{i}\frac{r_{i}(t)}{2}-\mathrm{i}H_{\mathrm{eff}}\Big{)}\ket{\psi},$
(A3)

with $\mathrm{d}N_{i}$ describing independent jump processes. For quantum
jumps $\mathrm{d}N_{i}$ take the values $1$ or $0$, depending on whether a
jump takes place or not. The average value

$\displaystyle\overline{\mathrm{d}N_{i}}=r_{i}(t)\mathrm{d}t$ (A4)

defines the jump rate $r_{i}(t)>0$, which is a free parameter because it later
cancels out in the ensemble average. The second term in Eq. A3 describes
deterministic evolution,

$\displaystyle\partial_{t}\ket{\psi^{(1)}}=\Big{(}\sum_{i}\frac{r_{i}(t)}{2}-\mathrm{i}H_{\mathrm{eff}}\Big{)}\ket{\psi^{(1)}},$
(A5)

where
$H_{\mathrm{eff}}=H-\frac{\mathrm{i}}{2}\sum_{i}\sigma_{i}L_{i}^{\dagger}L_{i}$
with anti Hermitian contribution is often regarded as effective Hamiltonian.

For the deterministic time evolution the sign is directly encoded in the
effective Hamiltonian, whereas for the quantum jump one cannot differentiate
between positive and negative relaxation strengths. This can be seen from the
fact that $-L_{i}\ket{\psi}\simeq L_{i}\ket{\psi}$ differ by a global phase
factor only and, thus, are equivalent. In other words for a pure state
$\outerproduct{\psi}{\psi}$ the relative sign in the jumps
$(-L_{i})\outerproduct{\psi}{\psi}(-L_{i})^{\dagger}$ would simply cancel. In
order to recover the relative sign the idea is to create an asymmetric
ensemble $\outerproduct*{\psi}{\tilde{\psi}}=s\outerproduct{\psi}{\psi}$ where
$\ket{\tilde{\psi}}$ is a copy of the trajectory $\ket{\psi}$, except that the
sign after a quantum jump is not changed,

$\displaystyle\ket{\mathrm{d}\tilde{\psi}}=\sum_{i}\mathrm{d}N_{i}\Big{(}\frac{L_{i}}{\sqrt{r_{i}(t)}}\ket{\tilde{\psi}}-\ket{\tilde{\psi}}\Big{)}+\mathrm{d}t\Big{(}\sum_{i}\frac{r_{i}(t)}{2}-\mathrm{i}H_{\mathrm{eff}}\Big{)}\ket{\tilde{\psi}}.$
(A6)

This is the same asymmetric unraveling proposed by Breuer and coworkers in
reference [22] for the special case that the trajectories only differ by a
relative sign. That is why it is not necessary to simulate
$\ket{\tilde{\psi}}$ but that it is sufficient to store the information about
the additional sign bit.

The rest of the proof follows the standard arguments of any quantum jump
unraveling. A nice overview is given by I. Kondov and coworkers in reference
[23]. To show that the unraveling yields on average the pseudo-Lindblad
equation (A2), we calculate an infinitesimal change of one signed trajectory,

$\displaystyle\mathrm{d}(s\outerproduct{\psi}{\psi})$
$\displaystyle=\mathrm{d}\outerproduct*{\psi}{\tilde{\psi}}=\outerproduct*{\mathrm{d}\psi}{\tilde{\psi}}+\outerproduct*{\psi}{\mathrm{d}\tilde{\psi}}+\outerproduct*{\mathrm{d}\psi}{\mathrm{d}\tilde{\psi}}.$
(A7)

Note that, for stochastic processes the quadratic term also contributes to
lowest order. The first two terms give

$\displaystyle\begin{split}\outerproduct*{\mathrm{d}\psi}{\tilde{\psi}}+\outerproduct*{\psi}{\mathrm{d}\tilde{\psi}}&=\sum_{i}\mathrm{d}N_{i}\Big{(}\frac{\sigma_{i}L_{i}\outerproduct*{\psi}{\tilde{\psi}}}{\sqrt{r_{i}(t)}}+\frac{\outerproduct*{\psi}{\tilde{\psi}}L_{i}^{\dagger}}{\sqrt{r_{i}(t)}}-\outerproduct*{\psi}{\tilde{\psi}}\Big{)}\\\
&-\sum_{i}\mathrm{d}N_{i}\outerproduct*{\psi}{\tilde{\psi}}+\mathrm{d}t\Bigg{(}\sum_{i}r_{i}(t)\outerproduct*{\psi}{\tilde{\psi}}-\mathrm{i}H_{\mathrm{eff}}\outerproduct*{\psi}{\tilde{\psi}}+\outerproduct*{\psi}{\tilde{\psi}}\mathrm{i}H_{\mathrm{eff}}^{\dagger}\Bigg{)},\end{split}$
(A8)

and essentially contribute to the effective Hamiltonian. Note that on average
the first terms in the second line cancels, since
$\overline{\mathrm{d}N_{i}}=r_{i}\mathrm{d}t$.

For the quadratic term
$\outerproduct*{\mathrm{d}\psi}{\mathrm{d}\tilde{\psi}}$ only stochastic
processes $\mathrm{d}N_{i}$ contribute. Since the processes are independent
and for single realizations one has $\mathrm{d}N_{i}=0,1$, it follows that
$\mathrm{d}N_{i}\mathrm{d}N_{j}=\delta_{ij}\mathrm{d}N_{i}$. We then get

$\displaystyle\begin{split}\outerproduct*{\mathrm{d}\psi}{\mathrm{d}\tilde{\psi}}&=\sum_{i}\mathrm{d}N_{i}^{2}\Big{(}\frac{\sigma_{i}L_{i}\ket{\psi}}{\sqrt{r_{i}(t)}}-\ket{\psi}\Big{)}\Big{(}\frac{\bra{\tilde{\psi}}L_{i}^{\dagger}}{\sqrt{r_{i}(t)}}-\bra{\tilde{\psi}}\Big{)}\\\
&=\sum_{i}\mathrm{d}N_{i}\
\sigma_{i}\frac{L_{i}\outerproduct*{\psi}{\tilde{\psi}}L_{i}^{\dagger}}{r_{i}(t)}-\sum_{i}\mathrm{d}N_{i}\Big{(}\frac{\sigma_{i}L_{i}\outerproduct*{\psi}{\tilde{\psi}}}{\sqrt{r_{i}(t)}}+\frac{\outerproduct*{\psi}{\tilde{\psi}}L_{i}^{\dagger}}{\sqrt{r_{i}(t)}}-\outerproduct*{\psi}{\tilde{\psi}}\Big{)},\end{split}$
(A9)

where the latter terms cancel the first line in Eq. A8, and we arrive at

$\displaystyle\mathrm{d}(\outerproduct*{\psi}{\tilde{\psi}})=\mathrm{d}t\,\Big{[}-\mathrm{i}H_{\mathrm{eff}}\outerproduct*{\psi}{\tilde{\psi}}+\outerproduct*{\psi}{\tilde{\psi}}\mathrm{i}H_{\mathrm{eff}}^{\dagger}\Big{]}+\sum_{i}\frac{\mathrm{d}N_{i}}{r_{i}(t)}\sigma_{i}L_{i}\outerproduct*{\psi}{\tilde{\psi}}L_{i}^{\dagger}.$
(A10)

By taking the ensemble average we replace
$\overline{\mathrm{d}N_{i}}=r_{i}(t)\mathrm{d}t$ and arrive at the pseudo-
Lindblad equation (A2) for the density matrix
$\varrho\equiv\overline{\outerproduct*{\psi}{\tilde{\psi}}}$.

##  II Norm conservation and choice of jump rates

As shown above the jump rates $r_{i}(t)$ are free parameters of the unraveling
because they cancel out in the ensemble average [23]. This is similar to the
relaxation strengths in the master equation Eq. A2, which can always be
redefined by rescaling the jump operators, as we employed it above to bring
the master equation to the form Eq. A2. However, the jump rates influence the
dynamics on the level of single trajectories and typically one chooses state-
dependent jump rates (see main text)

$\displaystyle r_{i}(t)=\frac{\lVert
L_{i}\psi(t)\rVert^{2}}{\lVert\psi(t)\rVert^{2}},$ (A11)

which yield quantum jumps that do not change the norm of the state:

$\displaystyle\ket{\psi(t)}\rightarrow\frac{L_{i}\ket{\psi(t)}}{\sqrt{r_{i}(t)}}=\frac{L_{i}\ket{\psi(t)}}{\lVert
L_{i}\psi(t)\rVert}\lVert\psi(t)\rVert.$ (A12)

Also for state-dependent jump rates the deterministic time evolution Eq. A5 is
non-linear, which for the MCWF compensates for the loss of norm under the
effective Hamiltonian. In turn, for master equations with negative signs the
norm increases according to

$\displaystyle\begin{split}\partial_{t}\lVert\psi^{(1)}\rVert^{2}&=\innerproduct*{\psi^{(1)}}{\partial_{t}\psi^{(1)}}+\innerproduct*{\partial_{t}\psi^{(1)}}{\psi^{(1)}}\\\
&=\underbrace{\matrixelement*{\psi^{(1)}}{-\mathrm{i}H}{\psi^{(1)}}+\matrixelement*{\psi^{(1)}}{\mathrm{i}H}{\psi^{(1)}}}_{=0}+\sum_{i}\matrixelement*{\psi^{(1)}}{(r_{i}(t)-\sigma_{i}L_{i}^{\dagger}L_{i})}{\psi^{(1)}}\\\
&=2\sum_{i,\sigma_{i}=-1}\matrixelement*{\psi^{(1)}}{L_{i}^{\dagger}L_{i}}{\psi^{(1)}},\end{split}$
(A13)

where the sum is taken over all channels with negative relaxation rate.
Alternatively, one may choose different jump rates $r_{i}(t)>0$ such that the
deterministic time evolution is norm preserving. The condition follows from
Eq. A13 to be

$\displaystyle\sum_{i}r_{i}(t)=\sum_{i}\sigma_{i}\frac{\lVert
L_{i}\psi\rVert^{2}}{\lVert\psi\rVert^{2}}.$ (A14)

For unitary jump operators the rates can be chosen state-independent.

For the eternal non-Markovian master equation discussed in the main text the
condition reads $\sum_{i}r_{i}(t)=1-\tanh(t)/2$, so that one possible choice
is

$\displaystyle r_{x}=r_{y}=\frac{1}{4},\qquad r_{z}=\frac{1-\tanh(t)}{2},$
(A15)

as compared to the choice in the main text [$r_{x}=r_{y}=1/2$,
$r_{z}=\tanh(t)/2$]. Again, the rates for $\sigma_{x}$ and $\sigma_{y}$-jumps
are taken to be equal. The $\sigma_{z}$ jump, which corresponds to the
negative relaxation strength $\gamma_{i}<0$, happens more often in the
beginning but vanishes in the static limit. Note, if the deterministic time
evolution is norm conserving, the quantum jumps are not. For the qubit this
would give quantum jumps $\ket{\psi}\rightarrow 2\sigma_{x}\ket{\psi}$,
$\ket{\psi}\rightarrow 2\sigma_{y}\ket{\psi}$ and
$\ket{\psi}\rightarrow\sqrt{2}\sigma_{z}\ket{\psi}/\sqrt{1-\tanh(t)}$. One
concludes that for negative relaxation strengths, there is no choice for jump
rates $r_{i}(t)>0$ for which the trajectories remain normalized. The increase
of the norm, however, is cancelled by trajectories with opposing sign, such
that on average the sign is preserved. This is shown in Fig. 4 (a), where we
plot the dynamics of the ensemble average
$\frac{1}{N}\sum_{n}^{N}s_{n}\langle\psi_{n}|\psi_{n}\rangle$ for increasing
number $N$ of trajectories. For proving this statement, note that on average
the ensemble of trajectories evolves according to the master equation, which
is trace-preserving, i.e.
$\tr\overline{s\,\outerproduct{\psi}{\psi}}=\overline{s\,\langle\psi|\psi\rangle}=1$.
For finite $N$ the fluctuations are on the order of $1/\sqrt{N}$ [Fig. 4 (b)].
They can be removed by explicitly normalizing the ensemble as mentioned in the
main text.

Figure 4: We simulate PLQTs for the eternal non-Markovian master equation
discussed in the main text. (a) Dynamics of ensemble average
$\overline{s\,\langle\psi|\psi\rangle}$ for increasing number $N$ of
trajectories (thicker lines). (b) Maximal deviation of the average from the
value $1$ during the time evolution as a function of ensemble size $N$.

##  III Dynamics of the averaged sign bit

Following the arguments in the main text, a trajectory
$\\{\ket{\psi(t)},s(t)\\}$ probabilistically changes its sign with the rate

$\displaystyle R(t)\equiv\sum_{i,\gamma_{i}<0}r_{i}(t),$ (A16)

where the sum runs over all jump rates $r_{i}(t)$ that correspond to a
negative $\gamma_{i}<0$. The dynamics of the sign then follows the update
scheme

$\displaystyle s(t+\delta t)=\begin{cases}-s(t)&\text{with probability
}R(t)\delta t\\\ s(t)&\text{with probability }1-R(t)\delta t.\end{cases}$
(A17)

For the averaged sign $\overline{s}(t)$ one takes the weighted average

$\displaystyle\overline{s}(t+\delta t)=-R(t)\delta
t\,\overline{s}(t)+\left(1-R(t)\delta t\right)\,\overline{s}(t)$ (A18)

from which we can deduce the first order differential equation

$\displaystyle\dot{\overline{s}}(t)=-2R(t)\ \overline{s}(t),$ (A19)

where we have used the limit $\dot{\overline{s}}(t)=\lim\limits_{\delta t\to
0}(s(t+\delta t)-s(t))/\delta t$.

Provided $R(t)$ is a known, smooth and state-independent function, the
averaged sign follows an exponential decay

$\displaystyle\overline{s}(t)=\exp[-2\int_{0}^{t}R(\tau)\mathrm{d}\tau\Big{]}.$
(A20)

This is the case, e.g., for the eternal master equation discussed in the main
text. However, even if the decay rate depends on the particular state,
qualitatively the net sign decays exponentially on a time scale that is given
by the negative relaxation rate. Therefore, eventually for sufficiently long
times the number of positive and negative trajectories will be balanced (up to
statistical fluctuations). As a result $\bar{s}$ vanishes for times that are
long compared to an _algorithmic relaxation time_.

$p_{+}$$p_{-}$rate $R$rate $R$ Figure 5: Let $R$ denote the rate under which
trajectories change their sign. For long times, the algorithm will always
equilibrate to a state having (up to statistical fluctuations) the same number
of trajectories with positive and negative sign. This can be readily seen from
the simple rate equation, $\dot{p}_{+}=-R\,(p_{+}-p_{-})$ and
$\dot{p}_{-}=-\dot{p}_{+}$, for the relative amount of positive and negative
trajectories, respectively.