# How to design quantum-jump trajectories via distinct master equation
representations

Dariusz Chruściński<EMAIL_ADDRESS>Institute of Physics, Faculty of
Physics, Astronomy and Informatics, Nicolaus Copernicus University,
Grudziadzka 5/7, 87-100 Toruń, Poland  Kimmo Luoma<EMAIL_ADDRESS>Institut für Theoretische Physik, Technische Universität Dresden, D-01062,
Dresden, Germany  Turku Center for Quantum Physics, Department of Physics and
Astronomy, University of Turku, FI-20014, Turun Yliopisto, Finland  Jyrki
Piilo<EMAIL_ADDRESS>Turku Center for Quantum Physics, Department of
Physics and Astronomy, University of Turku, FI-20014, Turun Yliopisto, Finland
Andrea Smirne<EMAIL_ADDRESS>Dipartimento di Fisica “Aldo
Pontremoli”, Università degli Studi di Milano, Via Celoria 16, I-20133 Milan,
Italy  Istituto Nazionale di Fisica Nucleare, Sezione di Milano, Via Celoria
16, I-20133 Milan, Italy

###### Abstract

Every open-system dynamics can be associated to infinitely many stochastic
pictures, called unravelings, which have proved to be extremely useful in
several contexts, both from the conceptual and the practical point of view.
Here, focusing on quantum-jump unravelings, we demonstrate that there exists
inherent freedom in how to assign the terms of the underlying master equation
to the deterministic and jump parts of the stochastic description, which leads
to a number of qualitatively different unravelings. As relevant examples, we
show that a fixed basis of post-jump states can be selected under some
definite conditions, or that the deterministic evolution can be fixed by a
chosen time-independent non-Hermitian Hamiltonian, even in the presence of
external driving. Our approach relies on the definition of rate operators,
whose positivity equips each unraveling with a measurement-scheme
interpretation and is related to a long known but so far not widely used
property to classify quantum dynamics, known as dissipativity. Starting from
formal mathematical concepts, our results allow us to get fundamental insights
into open quantum system dynamics and to ease their numerical simulations.

##  1 Introduction

Quantum jumps provide a powerful and insightful tool to describe the dynamics
of open quantum systems [1, 2]. In the quantum-jump unraveling, the open-
system state satisfying an assigned master equation is expressed as the
average of – in principle, infinitely many – trajectories of pure states,
which consist of deterministic evolutions interrupted by discontinuous,
randomly distributed jumps [3]. Quantum jumps have been observed in several
experimental platforms [4, 5, 6, 7, 8, 9] and are computationally convenient
to solve high-dimensional master equations, and so they are routinely used to
describe, e.g., quantum optical [10] or open many-body systems [11].
Interestingly, if the jumps are frequent and small enough, one can recover a
different stochastic description of the open-system evolution, characterized
by diffusive trajectories [12, 13, 14, 15, 16].

In the standard quantum-jump method, named Monte Carlo wave function (MCWF),
the state transformations induced by the jumps and their occurrence
probabilities are directly fixed by the operators and coefficients of the
master equation [3, 10]. This allows one to associate each trajectory with a
continuous selective measurement performed on the open system [17], so that
the master equation can be interpreted as due to the action of a non-selective
observer, replacing the environment. However, the very definition of MCWF
calls for a master equation with positive coefficients. Under some regularity
conditions, this requirement is equivalent to the completely-
positive(CP)-divisibility of the dynamics [18, 19, 20], meaning that the
dynamics can be decomposed into intermediate completely positive maps. CP-
divisibility has been introduced within the context of the definition of
quantum Markovianity [21, 22] and its validity implies the absence of memory
effects [23, 18, 24].

While generalizations of MCWF for master equations with negative coefficients
have been introduced [25, 26], the possibility to extend the continuous-
measurement interpretation beyond the realm of the CP-divisible evolutions has
been extensively debated [27, 28, 29]. Only very recently, a systematic
approach to read the quantum-jump unraveling of non-CP-divisible evolutions in
terms of continuous measurements has been put forward [30]. The approach
relies on the definition of a proper rate operator [31, 32, 33] and has hence
been named rate operator quantum jumps (ROQJ). It is associated with a
continuous-measurement interpretation applying to any positive(P)-divisible
dynamics [34, 35, 36, 24], i.e., dynamics that can be decomposed into
intermediate maps which are positive, but not necessarily completely positive.

In this paper, we define a class of jump unravelings interpolating between and
combining the advantages of MCWF and ROQJ. We exploit the freedom in dividing
any master equation into a deterministic and a jump part, possibly even mixing
the contributions from the Hamiltonian and the dissipative terms. Besides its
conceptual interest, this allows us to simplify the numerical and experimental
implementation of the trajectories. On the one hand, we can select a fixed set
of post-jump states for a well-defined class of master equations. On the
other, we can instead set a desired time-independent, non-Hermitian linear
operator ruling the deterministic evolution, also in the presence of an
external time-dependent driving on the open system. Importantly, our analysis
also shows that, under a constraint on P-divisible dynamics known as
dissipativity [37], there exist (infinitely) many unravelings that can be
associated with positive rate operators, bringing along different continuous-
measurement interpretations.

The rest of the paper is organized as follows. In Sec. 2, we present the
notation and the general notions of quantum-jump unravellings that will be
used throughout the paper. In Sec. 3, we introduce a whole class of rate-
operator unravelings, which interpolates between MCWF and ROQJ and allows for
a rather versatile control of the resulting trajectories. The positivity of
such unravelings is discussed in Sec. 4, where we prove that it is guaranteed
by dissipativity of the dynamics. Sec. 5 is devoted to an extended analysis of
a case study, a two-level system dynamics where CP-divisibility is broken at
any time, for which we compare the descriptions obtained via different ROQJ
unravelings. Finally, the conclusions of our work are discussed in Sec. 6.

##  2 Quantum-jump unravelings

We start off by recalling briefly the formalism of quantum jumps to describe
the dynamics of open quantum systems; in particular, we focus on the standard
quantum-jump unraveling, i.e., the MCWF method [3, 10], and the recently
introduced ROQJ [30].

###  2.1 Open quantum system dynamics

We consider the evolution of finite-dimensional open quantum systems, given by
the time-local master equation $d\rho(t)/(dt)=\mathcal{L}_{t}(\rho(t))$, where
$\mathcal{L}_{t}$ is the generator

$\mathcal{L}_{t}(\rho)=-i[H(t),\rho]+\mathcal{J}_{t}(\rho)-\frac{1}{2}\\{\Gamma(t),\rho\\},$
(1)

with
$\mathcal{J}_{t}(\rho)=\sum_{\alpha=1}^{N^{2}-1}c_{\alpha}(t)L_{\alpha}(t)\rho
L_{\alpha}^{\dagger}(t)$,
$\Gamma(t)=\mathcal{J}_{t}^{\dagger}(\mathbbm{1})=\sum_{\alpha=1}^{N^{2}-1}c_{\alpha}(t)L_{\alpha}^{\dagger}(t)L_{\alpha}(t)$
($\mathbbm{1}$ is the identity operator and $\mathcal{J}_{t}^{\dagger}$ is the
dual map of $\mathcal{J}_{t}$), $N$ the dimension of the open-system Hilbert
space $\mathcal{H}_{S}$, $H(t)=H^{\dagger}(t)$ and $L_{\alpha}(t)$ linear
operators on $\mathcal{H}_{S}$, and $c_{\alpha}(t)$ real functions of time.
The structure of the generator follows from the trace- and Hermiticity-
preservation properties of the dynamics [38]
$\Lambda_{t}=T\exp\left(\int_{0}^{t}d\tau\mathcal{L}_{\tau}\right)$ ($T$ is
the time-ordering operator); moreover, the functions $c_{\alpha}(t)$ can take
negative values, yet with the resulting evolution being positive and even
completely positive [1, 39]. On the other hand, under some regularity
conditions [18, 19, 20], the positivity of the coefficients,
$c_{\alpha}(t)\geq 0$ for every $\alpha$ and $t\geq 0$, is equivalent to the
CP-divisibility of the dynamics, i.e, for any $t\geq s\geq 0$ one has the
decomposition $\Lambda_{t}=V_{t,s}\Lambda_{s}$, where the so-called propagator
$V_{t,s}$ is completely positive and trace preserving; CP-divisibility has
been identified with quantum Markovianity in [21, 22].

A crucial remark for our following analysis is that the representation of the
generator $\mathcal{L}_{t}$ as in Eq.(1) is highly non-unique. In fact, one
can define a new map

$\mathcal{J}^{\prime}_{t}(\rho)=\mathcal{J}_{t}(\rho)+\frac{1}{2}\left(\mathbf{C}(t)\rho+\rho\mathbf{C}^{\dagger}(t)\right),$
(2)

with an arbitrary linear operator $\mathbf{C}(t)$ on $\mathcal{H}_{S}$, so
that decomposing

$\mathbf{C}(t)=A(t)+iB(t)$ (3)

(with Hermitian $A(t)$ and $B(t)$), together with

$H^{\prime}(t)=H(t)+\frac{1}{2}B(t),\quad\Gamma^{\prime}(t)=\Gamma(t)+A(t),$
(4)

the same generator $\mathcal{L}_{t}$ can be written as

$\mathcal{L}_{t}(\rho)=-i[H^{\prime}(t),\rho]+\mathcal{J}^{\prime}_{t}(\rho)-\frac{1}{2}\\{\Gamma^{\prime}(t),\rho\\};$
(5)

we also note that still
$\Gamma^{\prime}(t)=\mathcal{J}_{t}^{\prime\dagger}(\mathbbm{1})$.

###  2.2 Monte Carlo wave function vs rate operator quantum jumps

In the case of CP-divisible dynamics, the master equation fixed by Eq.(1) can
be obtained via a proper unraveling, the well-known MCWF, consisting of
discontinuous pure-state trajectories in $\mathcal{H}_{S}$.

The deterministic parts of the trajectories are fixed by the non-Hermitian
operator

$K(t)=H(t)-\frac{i}{2}\Gamma(t),$ (6)

according to

$|\psi(t)\rangle\to|\psi(t+dt)\rangle=\frac{(1-iK(t)dt)|\psi(t)\rangle}{\|(1-iK(t)dt)|\psi(t)\rangle\|},$
(7)

while the discontinuous parts, the jumps, are given by

$|\psi(t)\rangle\to|\psi(t+dt)\rangle=\frac{L_{\alpha}(t)|\psi(t)\rangle}{\|L_{\alpha}(t)|\psi(t)\rangle\|},$
(8)

and each jump occurs between $t$ and $t+dt$ with probability

$p_{\psi(t),\alpha}=c_{\alpha}(t)\|L_{\alpha}(t)|\psi(t)\rangle\|^{2}\,dt,$
(9)

where $dt$ is an infinitesimal time increment. It is clear that the previous
formulation requires all the rates $c_{\alpha}(t)$ to be positive.

Recently [30], a different quantum-jump unraveling, named ROQJ, has been
defined, weakening considerably such a requirement. ROQJ relies on the rate
operator

$\mathbf{W}_{\psi}=\sum_{\alpha=1}^{N^{2}-1}c_{\alpha}(t)(L_{\alpha}(t)-\ell_{\psi,\alpha})P_{\psi}(L_{\alpha}(t)-\ell_{\psi,\alpha}(t))^{\dagger},$
(10)

with $\ell_{\psi,\alpha}(t)=\langle\psi|L_{\alpha}(t)|\psi\rangle$ and the
projector $P_{\psi}=|\psi\rangle\langle\psi|$. We observe that
$\mathbf{W}_{\psi}$ is directly associated with the time-local generator
$\mathcal{L}_{t}$ and does not depend on its specific representation via
$H(t),\mathcal{J}_{t}$ and $\Gamma(t)$: in fact, it can be equivalently
written as

$\mathbf{W}_{\psi}=(\mathbbm{1}-P_{\psi})\mathcal{L}_{t}(P_{\psi})(\mathbbm{1}-P_{\psi}).$
(11)

The building block of ROQJ is the observation that $\mathbf{W}_{\psi}\geq 0$
for any state vector $|\psi\rangle$ if and only if the corresponding dynamics
is P-divisible [40], i.e., it can be decomposed as
$\Lambda_{t}=V_{t,s}\Lambda_{s}$, and the propagator $V_{t,s}$ is trace
preserving and positive, but not necessarily completely positive;
P-divisibility, which is a significantly weaker requirement than CP-
divisibility, has been identified with quantum Markovianity in [36, 24].
Focusing on P-divisible evolutions, one has then the spectral resolution

$\mathbf{W}_{\psi(t)}=\sum_{k=1}^{N}\lambda_{\psi(t),k}|\varphi_{\psi(t),k}\rangle\langle\varphi_{\psi(t),k}|,$
(12)

with $\lambda_{\psi(t),k}\geq 0$ for all $k$ and $t\geq 0$. Introducing the
non-Hermitian state-dependent operator

${K}_{\psi(t)}=K(t)+\Delta_{\psi(t)},$ (13)

with the non-linear correction

$\Delta_{\psi(t)}=\frac{i}{2}\sum_{\alpha=1}^{N^{2}-1}c_{\alpha}(t)(2L_{\alpha}(t)\ell^{*}_{\psi,\alpha}-|\ell^{*}_{\psi,\alpha}|^{2}),$

one realizes the jump unraveling as follows [30]: the free evolution

$|\psi(t)\rangle\to|\psi(t+dt)\rangle=\frac{(1-iK_{\psi(t)}dt)|\psi(t)\rangle}{\|(1-iK_{\psi(t)}dt)|\psi(t)\rangle\|}$
(14)

is interrupted by the sudden jumps

$|\psi(t)\rangle\to|\psi(t+dt)\rangle=\frac{V_{\psi(t),j}|\psi(t)\rangle}{\|V_{\psi(t),j}|\psi(t)\rangle\|},$
(15)

where

$V_{\psi(t),j}=\sqrt{\lambda_{\psi(t),j}}|\varphi_{\psi(t),k}\rangle\langle\psi(t)|;$
(16)

the probability that the jump $j$ occurs between $t$ and $t+dt$ is now

${p}^{\prime}_{\psi(t),j}=\|V_{\psi(t),j}|\psi(t)\rangle\|^{2}dt=\lambda_{\psi(t),j}dt.$
(17)

Compared to MCWF, the operators $L_{\alpha}(t)$ and rates $c_{\alpha}(t)$ are
replaced in ROQJ by the eigenvectors and eigenvalues of the rate operator
$\mathbf{W}_{\psi(t)}$. This is the key that allows one to formulate a well-
defined unraveling, with positive probabilities ${p}^{\prime}_{\psi(t),j}$,
for the set of P-divisible dynamics, rather than just for CP-divisible
dynamics as in MCWF. Both in MCWF [17] and in ROQJ [30] unraveling, the
trajectories can be interpreted as due to a continuous measurement on the open
system, i.e., the state transformations and corresponding probabilities can be
associated with a quantum instrument [41], mapping the set of outcomes into
the set of open-system completely positive trace non-increasing maps, but,
crucially, in ROQJ the continuous-measurement interpretation can be defined
for the whole set of P-divisible dynamics. Finally, both methods can be
extended to deal with, respectively, non-CP-divisible and non-P-divisible
evolutions, via the use of reversed jumps connecting different trajectories
[25, 26, 30]; however, ROQJ applies to general not necessarily P-divisible
evolutions [30], and it can be used to treat certain non-CP-divisible dynamics
[42, 43, 44] where the non-Markovian version of MCWF cannot be used.

Finally, we stress that the ROQJ involves the diagonalization of the rate
operator, which adds numerical complexity to the actual implementation of the
unraveling, but such a diagonalization needs not to be done at each time step.
Rather, one has to diagonalize the rate operator only if a jump occurs and
jumps are indeed rare events, whose probability is proportional to the
infinitesimal time step $dt$. Crucially, whether a jump occurs or not can be
fixed by looking only at the deterministic part of the evolution, as the
probability of having any jump is

$\displaystyle\sum_{j}{p}^{\prime}_{\psi(t),j}$
$\displaystyle=\sum_{j}\lambda_{\psi(t),j}dt=\mbox{Tr}\left\\{\mathbf{W}_{\psi(t)}\right\\}dt,$
(18)

which can be written as $1-\|(1-iK(t)dt)|\psi(t)\rangle\|^{2}$ (neglecting as
usual terms of order $dt^{2}$). Therefore, since we know the operator $K(t)$
giving the deterministic evolution, we need to construct the rate operator and
diagonalize it only at those (rare) times when a jump takes place. Note that
this applies not only to $\mathbf{W}_{\psi}$, but also to the rate operators
$\mathbf{R}_{\psi}$ described in the following.

##  3 A new class of rate operator quantum jumps

We are now ready to present a novel class of quantum-jump unravelings, which
is based on a family of rate operators defined starting from Eqs.(2)-(4). This
class combines features of both MCWF, allowing for a linear effective
Hamiltonian, and the ROQJ discussed above, being positive for dynamics with
possibly negative rates. In addition, we show how the freedom in choosing the
rate operator allows us to control and manipulate some basic properties of the
trajectories of the unraveling, which can be useful for the numerical
simulation of the dynamics and the actual experimental implementation of the
jumps.

###  3.1 Definition of the unravelings

The basic observation to define the class of ROQJs is that $\mathbf{W}_{\psi}$
in Eq.(11) can also be written as

$\mathbf{W}_{\psi}=(\mathbbm{1}-P_{\psi})\mathbf{R}_{\psi}(\mathbbm{1}-P_{\psi}),$
(19)

where

$\mathbf{R}_{\psi(t)}=\mathcal{J}_{t}(|\psi(t)\rangle\langle\psi(t)|);$ (20)

in other terms, decomposing the action of the generator $\mathcal{L}_{t}$ into
a block structure fixed by the projector $P_{\psi}$ into the current state
$\ket{\psi}$ and the orthogonal projector, $(\mathbbm{1}-P_{\psi})$, the rate
operator $\mathbf{W}_{\psi}$ equals the off-diagonal block, further restricted
to the sole contribution of the first term in the dissipative part of
$\mathcal{L}_{t}$. Now, if $\mathbf{R}_{\psi(t)}\geq 0$ for any
$|\psi(t)\rangle$ one can define a jump unraveling that merges (7)-(8) with
(14)-(15); that is, the free evolution is governed by (7), but the jumps are
realized via

$|\psi(t)\rangle\to|\psi(t+dt)\rangle=\frac{R_{\psi(t),k}|\psi(t)\rangle}{\|R_{\psi(t),k}|\psi(t)\rangle\|},$
(21)

with

$R_{\psi(t),k}=\sqrt{r_{\psi(t),k}}|\phi_{\psi(t),k}\rangle\langle\psi(t)|,$
(22)

and the probability that the jump $k$ occurs between $t$ and $t+dt$ reads

${p}^{\prime\prime}_{\psi(t),k}=\|R_{\psi(t),k}|\psi(t)\rangle\|^{2}dt=r_{\psi(t),k}\,dt;$
(23)

indeed, $r_{\psi(t),k}$ and $|\phi_{\psi(t),k}\rangle$ are the eigenvalues and
eigenvectors of $\mathbf{R}_{\psi(t)}$, i.e.,

$\mathbf{R}_{\psi(t)}=\sum_{k=1}^{N}r_{\psi(t),k}|\phi_{\psi(t),k}\rangle\langle\phi_{\psi(t),k}|.$
(24)

General conditions guarantying the positivity of $\mathbf{R}_{\psi(t)}$ will
be presented in Sec. 4.

Before proving that this construction provides us with a well-defined
unraveling of the master equation (1) whenever the rate operator is positive,
let us stress that it avoids the non-linear correction in the non-Hermitian
operator fixing the deterministic evolution, still potentially allowing for
positive probabilities in the case of non-CP-divisible evolutions;
${p}^{\prime\prime}_{\psi(t),k}\geq 0$ is in fact equivalent to the
requirement that the map $\mathcal{J}_{t}$ is positive for all $t\geq 0$,
which is significantly weaker than CP-divisibility as we will see in Sec. 4.
Moreover, the definition of the rate operator $\mathbf{R}_{\psi(t)}$ in
Eq.(20) does depend on the specific choice of the map $\mathcal{J}_{t}$ in the
representation of the generator $\mathcal{L}_{t}$ as in Eq.(1); as a
consequence, Eqs.(7), (21) and (23) define a whole family of unravelings,
corresponding to different choices of the operator $\mathbf{C}(t)$ in Eqs.(2)
and (4). We denote these unravelings with $\mathbf{R}$-ROQJ, specifying the
explicit form of the rate operator $\mathbf{R}$ when needed; furthermore, we
will refer to the ROQJ unraveling discussed in Sec. 2.2 as $\mathbf{W}$-ROQJ.

###### Proof.

Given the pure state $|\psi(t)\rangle\langle\psi(t)|$ at time $t$, the
deterministic evolution will occur with probability

$\displaystyle P^{det}_{\psi(t)}$
$\displaystyle=1-\sum_{k=1}^{N}{p}^{\prime\prime}_{\psi(t),k}=1-\sum_{k=1}^{N}r_{\psi(t),k}\,dt$
$\displaystyle=1-\mbox{Tr}\left\\{\mathbf{R}_{\psi(t)}\right\\}dt,$ (25)

where in the last equality Tr denotes the trace and we used Eq.(24). But using
Eq.(20) along with

$\Gamma(t)=\mathcal{J}_{t}^{\dagger}(\mathbbm{1}),$ (26)

we get

$\displaystyle P^{det}_{\psi(t)}$
$\displaystyle=1-\mbox{Tr}\left\\{\mathcal{J}_{t}(|\psi(t)\rangle\langle\psi(t)|)\right\\}dt$
$\displaystyle=1-\mbox{Tr}\left\\{|\psi(t)\rangle\langle\psi(t)|\mathcal{J}^{\dagger}_{t}(\mathbbm{1})\right\\}dt$
$\displaystyle=1-\langle\psi(t)|\Gamma(t)|\psi(t)\rangle dt.$ (27)

In addition, the deterministic evolution will map the pure state
$|\psi(t)\rangle\langle\psi(t)|$ into the pure state (see Eq.(7))

$\displaystyle\frac{(1-iK(t)dt)|\psi(t)\rangle\langle\psi(t)|(1+iK^{\dagger}(t)dt)}{\|(1-iK(t)dt)|\psi(t)\rangle\|^{2}}$
$\displaystyle=\frac{(1-iK(t)dt)|\psi(t)\rangle\langle\psi(t)|(1+iK^{\dagger}(t)dt)}{1-\langle\psi(t)|\Gamma(t)|\psi(t)\rangle
dt}$
$\displaystyle=\frac{(1-iK(t)dt)|\psi(t)\rangle\langle\psi(t)|(1+iK^{\dagger}(t)dt)}{P^{det}_{\psi(t)}},$
(28)

where in the first equality we used Eq.(6) in the denominator (and neglected
the terms of order $dt^{2}$), while in the second equality we used Eq.(27).

On the other hand, as said, given the state $|\psi(t)\rangle\langle\psi(t)|$
at time $t$, we will have the jump described by $R_{\psi(t),k}$ in Eq.(22),
i.e. (compare with Eq.(21))

$|\psi(t)\rangle\langle\psi(t)|\to\frac{R_{\psi,k}|\psi(t)\rangle\langle\psi(t)R^{\dagger}_{\psi,k}}{\|R_{\psi,k}|\psi(t)\rangle\|^{2}}=|\phi_{\psi(t),k}\rangle\langle\phi_{\psi(t),k}|,$
(29)

with probability as in Eq.(23).

All in all, if we average the state at time $t+dt$ over the trajectories where
the state at time $t$ is $|\psi(t)\rangle\langle\psi(t)|$, we get the mixture
of the states obtained via the deterministic evolution or one of the jumps,
each weighted with the corresponding probability: using Eqs.(28), (29) and
(23), such a (conditioned) average corresponds to

$\displaystyle
P^{det}_{\psi(t)}\frac{(1-iK(t)dt)|\psi(t)\rangle\langle\psi(t)|(1+iK^{\dagger}(t)dt)}{P^{det}_{\psi(t)}}$
$\displaystyle+\sum_{k=1}^{N}r_{\psi(t),k}|\phi_{\psi(t),k}\rangle\langle\phi_{\psi(t),k}|\,dt$
$\displaystyle=|\psi(t)\rangle\langle\psi(t)|-i\left[H,|\psi(t)\rangle\langle\psi(t)|\right]dt$
$\displaystyle-\frac{1}{2}\\{\Gamma(t),|\psi(t)\rangle\langle\psi(t)|\\}dt+\mathcal{J}_{t}(|\psi(t)\rangle\langle\psi(t)|)dt,$
(30)

where the equality is due to Eqs.(6), (20) and (24), and we neglected the
terms of order $dt^{2}$. Finally, we perform a second average, this time with
respect to the pure states $|\psi(t)\rangle\langle\psi(t)|$ we fixed at time
$t$, so that we get the average over all the trajectories; the previous
expression then yields the first order expansion of the time-local master
equation $d\rho(t)/(dt)=\mathcal{L}_{t}(\rho(t))$ with $\mathcal{L}_{t}$ as in
Eq.(1), which concludes the proof. ∎

Indeed, the proof above does not depend on the specific representation of the
generator $\mathcal{L}_{t}$ and then it holds for any choice of the operator
$\mathbf{C}$, as long as the corresponding rate operator
$\mathbf{R}_{\psi(t)}$ is positive. Quite remarkably, we can exploit the
different legitimate rate operators to control the deterministic and jump
parts of the unraveling.

###  3.2 $\mathbf{R}$-ROQJ with fixed post-jump states

First of all, a proper choice of $\mathbf{C}$ can be used to simplify to a
significant extent the jumps in the unraveling. An explicit condition for this
can be derived for two-dimensional open quantum systems, i.e.,
$\mathcal{H}_{S}=\mathbbm{C}^{2}$.

Given any time-local generator $\mathcal{L}_{t}$ on the set of linear
operators $\mathcal{B}(\mathbbm{C}^{2})$ in the GKLS form as in Eq.(1),
$\mathcal{J}_{t}$ can be represented via the 16 parameters $J^{kl}_{ij}(t)$
defined by

$J^{kl}_{ij}(t)=\bra{i}\mathcal{J}_{t}[\ket{k}\bra{l}]\ket{j}\quad i,j=1,2,$
(31)

with respect to any orthonormal basis $\left\\{\ket{1},\ket{2}\right\\}$; note
that the Hermiticity-preservation condition implies
$\left(J^{kl}_{ij}(t)\right)^{*}=J^{lk}_{ji}(t)$. Now, if there is a basis
such that

$J^{kl}_{ij}(t)=(J^{kl}_{ij}(t))^{*}\quad\forall i,j,k,l=1,2,$ (32)

i.e., the Choi matrix associated with $\mathcal{J}_{t}$ is real, and the
deterministic evolution in-between the jumps does not introduce a relative
phase when acting on the basis elements, it is possible to define a rate
operator $\mathbf{R}_{\psi(t)}$ via Eqs.(2) and (20) whose eigenvectors are
independent from the pre-jump state $\ket{\psi(t)}$. This implies that the
post-jump states are always the same, so that the trajectories are fixed by at
most 3 deterministically-evolving states, which results in a dramatic
simplification of both the numerical simulation and the experimental
implementation of the corresponding unraveling, as will also be shown by means
of an example in Sec.5. More explicitely, we have in fact the following.

######  Proposition 1.

Given the master equation (1) for $N=2$, if there is an orthonormal basis of
$\mathbbm{C}^{2}$, denoted as
$\left\\{\ket{\varphi_{1}},\ket{\varphi_{2}}\right\\}$, and a time-dependent
linear operator $\mathbf{C}(t)\in\mathcal{B}(\mathbbm{C}^{2})$ with matrix
representation in this basis

$\displaystyle\mathbf{C}(t)=$ (33)
$\displaystyle\left(\begin{array}[]{cc}J^{11}_{22}(t)-J^{11}_{11}(t)+ix(t)&y(t)\\\
y(t)+2(J^{12}_{11}(t)-J^{12}_{22}(t))&J^{22}_{11}(t)-J^{22}_{22}(t)+ix(t)\end{array}\right)$
(36)

with $x(t),y(t)$ two real functions, such that

  * •

the jump operator $\mathcal{J}^{\prime}_{t}$ fixed by Eq.(2) satisfies
Eq.(32);

  * •

if the state $\ket{\psi(t)}$ at time $t$ is in the form

$\ket{\psi(t)}=c(t)\ket{\varphi_{1}}\pm\sqrt{1-c(t)^{2}}\ket{\varphi_{2}}$
(37)

for some $-1\leq c(t)\leq 1$, the deterministic evolution fixed by Eqs.(4),
(6) (with $H^{\prime}$ instead of $H$ and $\Gamma^{\prime}$ instead of
$\Gamma$) and (7) implies that the state at time $t+dt$ is in the form

$\displaystyle\ket{\psi(t+dt)}=a(t,dt)\ket{\varphi_{1}}\pm\sqrt{1-a(t,dt)^{2}}\ket{\varphi_{2}}$
(38)

for some $-1\leq a(t,dt)\leq 1$;

  * •

the rate operator
$\mathbf{R}_{\psi(t)}=\mathcal{J}^{\prime}_{t}(|\psi(t)\rangle\langle\psi(t)|)$
is positive for any state $\ket{\psi(t)}$ as in Eq.(37);

and we restrict to

  * •

initial states $\ket{\psi_{0}}=\ket{\psi(0)}$ in the form as in Eq.(37),

then there exists a jump unraveling that involves only the 3 families of
states
$\left\\{\ket{\psi_{0}(t)},\ket{\varphi_{+}(t,s)},\ket{\varphi_{-}(t,s)}\right\\}$,
which are deterministically evolved from
$\left\\{\ket{\psi_{0}},\ket{\varphi_{\pm}}=\frac{1}{\sqrt{2}}(\ket{\varphi_{1}}\pm\ket{\varphi_{2}})\right\\}$
via

$\displaystyle D(t,s)=$ $\displaystyle
T\exp\left(-i\int_{s}^{t}d\tau\left(H^{\prime}(\tau)-\frac{i}{2}\Gamma^{\prime}(\tau)\right)\right).$
(39)

The conditions in Eqs.(37) and (38) mean that the unraveling will involve
exclusively pure states without a relative phase in the basis
$\left\\{\ket{\varphi_{1}},\ket{\varphi_{2}}\right\\}$, if this is the case at
the initial time; this is also why the positivity requirement on the rate
operator $\mathbf{R}_{\psi(t)}$ can be now restricted only to states of this
form. In practice, the deterministic evolution between two times $s$ and $t$
will be given by the operator in Eq.(39), which indeed corresponds to the
effective non-Hermitian Hamiltonian fixed by Eq.(4) for the operator
$\mathbf{C}(t)$ as in Eq.(33). The jumps at time $t$, given by the
eigenvectors of $\mathbf{R}_{\psi(t)}$, will always end up in one of the
eigenvectors $\left\\{\ket{\varphi_{+}},\ket{\varphi_{-}}\right\\}$, and then
all in all the trajectories will consist of piecewise deterministic evolutions
among
$\left\\{\ket{\psi_{0}(t)},\ket{\varphi_{+}(t,s)},\ket{\varphi_{-}(t,s)}\right\\}$.
Note that the probabilities of having the jumps between time $t$ and $t+dt$
are fixed by the eigenvalues of $\mathbf{R}_{\psi(t)}$, so that different
matrices $\mathbf{C}(t)$ in Eq.(33) will generally define distinct
unravelings, with the same trajectories but different associated
probabilities.

In addition, the existence of a class of ROQJ can be exploited to control to a
large extent the deterministic evolution in-between the jumps. From Eqs.(3),
(4) and (6) we see how any effective non-Hermitian Hamiltonian can be enforced
by means of a proper choice of $\mathbf{C}$, which will also result in a
modified set of jump states and probabilities, according to Eq.(2), (20) and
(23). As will be shown by means of a significant example in Sec.5, this choice
can be made while keeping the positivity of the rate operator and it is
particularly convenient for driven systems.

##  4 Positivity of the rate operators

We can now complete the definition of the $\mathbf{R}$-ROQJ by deriving
general conditions ensuring that it is associated with a positive rate
operator. Besides justifying the construction in Eqs.(21) and (23), this also
guarantees the existence of a fully consistent continuous-measurement
interpretation [30].

Clearly, if the evolution is CP-divisible, not only both MCWF and
$\mathbf{W}$-ROQJ are well defined, but one can always find a completely
positive map $\mathcal{J}_{t}$ giving rise to a positive $\mathbf{R}$-ROQJ.
Instead, any P-divisible evolution does guarantee positivity of
$\mathbf{W}_{\psi(t)}$, but not necessarily the positivity of
$\mathbf{R}_{\psi(t)}$. Still, we note that P-divisibility constrains the
number of possible negative eigenvalues of the rate operator
$\mathbf{R}_{\psi(t)}$. We have in fact the following:

######  Proposition 2.

For any P-divisible evolution, $\mathbf{R}_{\psi(t)}$ can have at most one
negative eigenvalue.

###### Proof.

The proof easily follows from the well known min-max principle for Hermitian
matrices: let $A$ be a Hermitian $n\times n$ matrix, and let

$\lambda_{n}\geq\lambda_{n-1}\geq\ldots\geq\lambda_{1}$

be the real eigenvalues of $A$. Then

$\lambda_{k}=\max_{\Sigma}\,\min_{x\in\Sigma}\langle x|A|x\rangle$ (40)

where $x$ is normalized, and $\Sigma$ is a $(n-k+1)$-dimensional subspace of
$\mathbb{C}^{n}$.

Now, since
$\mathbf{W}_{\psi}=(\mathbbm{1}-P)\mathbf{R}_{\psi}(\mathbbm{1}-P)\geq 0$, one
has for any $x\in\Sigma=(\mathbbm{1}-P)\mathbb{C}^{n}$

$\langle x|\mathbf{R}_{\psi}|x\rangle=\langle x|\mathbf{W}_{\psi}|x\rangle\geq
0.$

Hence, from (40) one finds $\lambda_{2}\geq 0$ (since in (40) one maximizes
over all $\Sigma$s), and hence only $\lambda_{1}$ may be negative. ∎

We stress that the proof does not depend on the specific $\mathcal{J}_{t}$
used to represent the generator $\mathcal{L}_{t}$ and hence applies to any
$\mathbf{R}$-ROQJ.

Now, since $\mathbf{R}_{\psi(t)}$ does depend upon the specific
$\mathcal{J}_{t}$, see Eq.(20), a natural strategy to ensure the positivity of
$\mathbf{R}_{\psi(t)}$ arises: is it possible to use the freedom introduced by
(2) so that the rate operator defined in terms of $\mathcal{J}^{\prime}_{t}$
is positive? Interestingly, in all the examples we analyzed this is the case;
however, we can give a positive answer in generality only for a class of
P-divisible evolutions that enjoys an additional property. In the Heisenberg
picture P-divisibility means that the propagators $V_{t,s}^{\dagger}$ are
positive and unital. It is well known [45, 46] that this implies the Kadison-
Schwarz inequality

$V_{t,s}^{\dagger}(X^{\dagger}X)\geq
V_{t,s}^{\dagger}(X^{\dagger})V_{t,s}^{\dagger}(X),$ (41)

for all Hermitian operators $X=X^{\dagger}\in{\mathcal{B}(\mathcal{H}_{S})}$.
A more restricted class of evolutions is thus defined by propagators that
satisfy (41) for all $X$, not necessarily Hermitian. This leads us to the
identification of the condition ensuring the positivity of the
$\mathbf{R}$-ROQJ:

######  Proposition 3.

If the propagators of the dynamics satisfy (41) for all
$X\in{\mathcal{B}(\mathcal{H}_{S})}$, there exists a representation of the
master equation (1) with a positive map $\mathcal{J}_{t}$, which we denote as
$\mathbf{J}_{t}$.

The proof, which is reported in Appendix B, is based on some techniques
introduced by Lindblad in his seminal paper [37], based on the fact that the
Kadison-Schwarz inequality may be rephrased by the following condition for the
time-local generator (in the Heisenberg picture):

$\mathcal{L}^{\dagger}_{t}(X^{\dagger}X)\geq\mathcal{L}^{\dagger}_{t}(X^{\dagger})X+X^{\dagger}\mathcal{L}^{\dagger}_{t}(X),$
(42)

again for all $X\in{\mathcal{B}(\mathcal{H}_{S})}$ and $t\geq 0$. The
generators satisfying Eq.(42) are called dissipative [37]. We stress that the
dissipativity condition is equivalent to P-divisibility if one restricts to
Hermitian operators $X$; however, assuming that the condition (42) holds for
all $X\in{\mathcal{B}(\mathcal{H}_{S})}$ the corresponding evolution is not
only P-divisible but in addition the propagator $V_{t,s}$ satisfies (41).
Indeed, if $X=P_{\psi}$ is a rank-1 projector, then $P_{\psi}^{2}=P_{\psi}$
and (42) implies $({\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}-P_{\psi})\mathcal{L}^{\dagger}_{t}(P_{\psi})({\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}-P_{\psi})\geq 0$, which coincides with the condition for
P-divisibility (in Heisenberg picture). From the physical point of view, the
dissipativity condition can be understood with the following remark: let
$\varrho(t)$ be an instantaneous invariant state, i.e.
$\mathcal{L}_{t}(\varrho(t))=0$. One has ${\rm
Tr}(\varrho(t)\mathcal{L}^{\dagger}_{t}(X^{\dagger}X))={\rm
Tr}(\mathcal{L}_{t}(\varrho(t))(X^{\dagger}X))=0$, and hence introducing the
following inner product in the space of Hermitian operators

$(X,Y)_{\varrho}={\rm Tr}(\varrho X^{\dagger}Y)$

the formula (42) implies [20]

${\rm Re}\,(X,\mathcal{L}^{\dagger}_{t}(X))_{\varrho(t)}\leq 0,$ (43)

for any operator $X\in{\mathcal{B}(\mathcal{H}_{S})}$. In particular, if
$|i(t)\rangle$ defines an eigen-basis of $\varrho(t)$, then taking
$X=|i(t)\rangle\langle j(t)|$ the condition (43) clearly shows that
P-divisibility provides constraints for populations $(i=j)$, whereas
dissipativity is much more restrictive providing constraints also for
coherences $(i\neq j$).

Summarizing, the condition in Eq.(41), or equivalently Eq.(42), for all
$X\in{\mathcal{B}(\mathcal{H}_{S})}$ allows for a representation of (1) with a
positive map $\mathbf{J}_{t}$ (see Appendix B for technical details). In other
terms, dissipativity guarantees the existence of a rate operator
$\mathbf{R}_{\psi(t)}$, defined from $\mathbf{J}_{t}$ by Eq.(20), such that
the unraveling given by Eqs.(7), (21) and (23) is well-defined, as
$\mathbf{R}_{\psi(t)}\geq 0$ for any $\psi(t)$. It should be stressed that the
representation with a positive map $\mathbf{J}_{t}$ is not unique,
corresponding to a non-unique continuous-measurement scheme, as we will see in
the examples in the next section.

##  5 Eternally non-Markovian qubit master equation

To explore in an explicit case study the different possible $\mathbf{R}$-ROQJ
unravelings and compare them among each other, as well as with the
$\mathbf{W}$-ROQJ, we consider the two-level system dynamics fixed by the
following generator

$\mathcal{L}_{t}(\rho)=i\frac{b(t)}{2}[\sigma_{z},\rho]+\frac{1}{2}\sum_{k=1}^{3}\gamma_{k}(t)(\sigma_{k}\rho\sigma_{k}-\rho),$
(44)

where $\sigma_{k}$ are the Pauli spin operators. Such a master equation has
been studied extensively in the literature [42, 43, 44], since, despite its
simplicity, it possesses several interesting features. In particular, for some
choices of the rates $\gamma_{k}(t)$ the resulting dynamics is P-divisible,
while CP-divisibility is broken for any $t>0$. This kind of dynamics is
usually referred to as eternally non-Markovian and, most importantly for us,
it cannot be treated by means of the MCWF.

###  5.1 Undriven master equation

We set at first $b(t)=0$ and study the case where we have a P-divisible
evolution even if one of the rates is temporally negative, i.e., CP-
divisibility is broken; we suppose in particular that $\gamma_{3}(t)<0$ for
any $t>0$. Then $\Lambda_{t}$ is P-divisible provided that
$\gamma_{1}(t),\gamma_{2}(t)\geq|\gamma_{3}(t)|$. Dissipativity requires the
stronger condition [47] that $\gamma_{1}(t),\gamma_{2}(t)\geq
2|\gamma_{3}(t)|$. Interestingly, it turns out that whenever $\gamma_{3}(t)<0$
the map
$\mathcal{J}_{t}(\rho)=\sum_{k=1}^{3}\gamma_{k}(t)\sigma_{k}\rho\sigma_{k}/2$
is not positive. However, if the generator (44) gives rise to a P-divisible
evolution, one can define different positive $\mathbf{R}$-ROQJ by means of
different operators $\mathbf{C}(t)$ in Eqs.(2)-(4).

#### Unraveling with $\mathbf{R1}$

The first choice we make is to use

$\mathbf{C}(t)=\frac{\gamma(t)}{2}{\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}$ (45)

with

$\gamma(t)=\sum_{k}\gamma_{k}(t),$

as the corresponding map

$\mathbf{J}_{t}(\rho)=\mathcal{J}_{t}(\rho)+\frac{\gamma(t)}{2}\rho$ (46)

is positive, as shown in Appendix C. In particular for the eternally non-
Markovian evolution [42, 44] defined by $\gamma_{1}=\gamma_{2}=1$, and
$\gamma_{3}(t)=-{\rm tanh}\,t$, one has a P-divisible evolution, but the
corresponding generator $\mathcal{L}_{t}^{\dagger}=\mathcal{L}_{t}$ is not
dissipative. Nevertheless, the map $\mathbf{J}_{t}$ is positive and one can
define the positive rate operator
$\mathbf{R1}_{\psi(t)}=\mathbf{J}_{t}(|\psi(t)\rangle\langle\psi(t)|)\geq 0$.
The deterministic evolution is fixed by
$K1(t)=\frac{i}{2}\gamma(t)\mathbbm{1}$. Figure 1 (a) shows 7 example pure-
state trajectories, with initial state
$|\psi(0)\rangle=\sqrt{0.1}|1\rangle+\sqrt{0.9}|2\rangle$, and the probability
$\rho_{11}$ of the state $|1\rangle$. Figure 1 (b) displays the final
distribution of the Bloch vector X and Z components over the trajectories
while the inset shows the agreement between the analytical results and
simulations for the coherence $\rho_{12}$.

Figure 1: Simulation results using rate operator $\mathbf{R1}$ [(a), (b)] and
$\mathbf{R2}$ [(c), (d)] for model (44) without driving. Example trajectories
are displayed in (a) and (c). Final steady state distributions for Bloch
vector X and Z components are shown in (b) and (d). The insets in (b) and (d)
demonstrate the agreement between the analytical results (solid line) and
simulation results (dots) with $10^{4}$ trajectories. The error bars are
similar to the size of the dots.

#### Unraveling with $\mathbf{R2}$

The second choice we make takes advantage of the general result for qubits
dynamics in Proposition 1. In fact, as we show explicitly in Appendix D, the
eternal non-Markovian evolution satisfies the assumptions of the mentioned
proposition, meaning that it is possible to define a $\mathbf{R}$-ROQJ with a
fixed set of (deterministically evolving) states after the jumps. Actually, in
Appendix D we define a continuous family of rate operators with fixed post-
jump states for the eternal non-Markovian dynamics, also ensuring their
positivity. Quite interestingly, this means that there exists a continuous
family of continuous measurement schemes for the same open-system dynamics
originating from the non-unique decomposition (5) of the generator. We stress
that this freedom is inherently different from the well-known [1] unitary
freedom in defining the operators that appear in the master equation (1).
Moreover, the instruments defining the different measurement schemes, as well
as the associated probabilities, will generally be different, while the post-
measurement states will be the same.

A special instance of the $\mathbf{R}$-ROQJs with fixed post-jump states is
given by

$\mathbf{R2}_{\psi(t)}=\mathbf{R1}_{\psi(t)}-\gamma_{3}(t)\ket{\psi(t)}\bra{\psi(t)}$
(47)

associated with the deterministic evolution
$K2(t)=\frac{i}{2}[\gamma_{1}(t)+\gamma_{2}(t)]\mathbbm{1}$. The corresponding
simulation results are shown in Fig. 1 (c) and (d). The example trajectories
in Fig. 1 (c) illustrate that the ensemble consists now a discrete set of pure
states. In the final distribution of states, see Fig. 1 (d), we have only two
states $|\pm\rangle=\frac{1}{\sqrt{2}}(|1\rangle\pm|2\rangle)$. As a matter of
fact, during the whole simulation, only three states appear: the initial state
and $|\pm\rangle$; in other terms, in this case the 3 states fixing the
unraveling according to Proposition 1 are even time independent. Remarkably,
$\mathbf{R2}_{\psi(t)}$ explicitly demonstrates that the measurement basis can
be fixed once and for all, without any state and time dependence.

#### Unraveling with $\mathbf{R3}$ and $\mathbf{W}$

Figure 2: Simulation results using rate operator $\mathbf{R3}$ [(a), (b)] and
$\mathbf{W}$ [(c), (d)]. Seven example pure state trajectories are displayed
in (a) and (c) by plotting the probability of state $|1\rangle.$ Final steady
state distributions for Bloch vector X and Z components, with 500
trajectories, are shown in (b) and (d). Here, each dot corresponds to a single
trajectory. In (d), all the 500 dots have the value $Z=0$ and $X=-1/2$ or
$X=+1/2$, i.e., we have only two distinctively visible dots. The insets in (b)
and (d) demonstrate the agreement between the analytical and simulation
results by using the coherence $\rho_{12}$ and having $10^{4}$ trajectories.
In the insets, the solid line is the analytical result and the dots are the
simulation results. The error bars are similar to the size of the dots. In all
the panels, the initial state is
$|\psi(0)\rangle=\sqrt{0.1}|1\rangle+\sqrt{0.9}|2\rangle$ and the used time
step size is $dt=0.002$.

Finally, we consider the rate operator

$\mathbf{R3}_{\psi}=\frac{1}{2}\sum_{k=1}^{3}\gamma_{k}(t)\sigma_{k}|\psi(t)\rangle\langle\psi(t)|\sigma_{k}-\frac{1}{2}\gamma_{3}(t)|\psi(t)\rangle\langle\psi(t)|,$
(48)

which is positive as shown in Appendix C. The corresponding deterministic
evolution is fixed by the linear operator
$K3(t)=\frac{i}{4}[\gamma_{1}(t)+\gamma_{2}(t)]\mathbbm{1}$. We show the
simulation results in Fig. 2 (a) and (b). The example realizations in Fig. 2
(a) demonstrate that, similarly to $\mathbf{R1}$ (cf. Fig. 1), the jumps
continue even though the steady state has been reached. However, contrary to
$\mathbf{R1}$, now it holds all the time that $\rho_{11}\leqslant 0.1$ or
$\rho_{11}\geqslant 0.9$. This can be seen more clearly in Fig. 2 (b). While
the distribution of trajectories in the steady state for $\mathbf{R1}$
contained arcs on the western and eastern sides of the XZ-plane of the Bloch
sphere, now with $\mathbf{R3}$ the distribution covers arcs on the northern
and southern sides of the circle. As before, there is excellent agreement
between the analytical and simulation results [see the inset of Fig. 2 (b).]

The last operator we consider is $\mathbf{W}$, see Eq. (10). We display the
simulation results in Fig. 2 (c) and (d). With the example realizations in
Fig. 2 (c), one can clearly see that now also the deterministic evolution
changes the states. Moreover, the jumps happen between a pair of states only
and terminate when the steady state is reached. Fig. 2 (d) shows that
similarly to $\mathbf{R2}$ [cf. Fig. 1 (d)], all the trajectories eventually
end up being on one of two states on the equator of the Bloch sphere. However,
how they reach these points is totally different with respect to
$\mathbf{R2}$. In terms of the measurement scheme, there is another crucial
difference between $\mathbf{W}$ and $\mathbf{R2}$. With $\mathbf{W}$ the post-
measurement states are time dependent, while with $\mathbf{R2}$ they are time
independent.

Table 1 collects and compares the properties of all the used four operators in
terms of (i) whether the quantum jumps continue also in the asymptotic regime
or terminate when the steady state is reached, (ii) whether the deterministic
evolution changes the state of the trajectories, (iii) whether the post-
measurement states are time dependent or time independent, and (iv) what is
the effective ensemble size in the simulation (how many different kinds of
state vectors the ensemble consists of point-wise in time). Overall, the rate
operator $\mathbf{R2}$ has the most appealing properties from simulation and
fundamental interpretation points of views.

Rate | Asymptotic | Deterministic | Time-independent | Effective ensemble
---|---|---|---|---
operator | jumps? | changes? | post-measurement states? | size
$\mathbf{R1}$ | yes | no | no | $\infty$
$\mathbf{R2}$ | no | no | yes | 3
$\mathbf{R3}$ | yes | no | no | $\infty$
$\mathbf{W}$ | no | yes | no | 2
Table 1:  The basic characterization of the four used rate operators for
stochastic simulations. Quantum jumps can either continue or terminate when
steady state is reached. Depending on the form of the effective Hamiltonian
giving the deterministic evolution, the state of the trajectory can either
change or remain unchanged after renormalization. The states after the
measurements in the continuous-measurement scheme interpretation can be either
time dependent or time independent. The effective ensemble size describes how
many different kinds of state vectors the ensemble consists of point-wise in
time.

###  5.2 Driven master equation

As a second example, we add a time-dependent driving $b(t)$ to the evolution
(44) and consider the decay rates $\gamma_{1}=\gamma_{2}=1$ and
$\gamma_{3}(t)=-\frac{1}{2}\tanh t$. For the driving, we choose

$b(t)=C+\int\limits_{0}^{t}\mathrm{d}s\,\Omega(s;\mu,\sigma)$ (49)

including an integral over the Gaussian function

$\Omega(s;\mu,\sigma)=\frac{1}{\sqrt{2\pi}\sigma}e^{-(s-\mu)^{2}/(2\sigma^{2})}$

and a constant

$C=\int_{-\infty}^{0}\mathrm{d}s\,\Omega(s;\mu,\sigma).$

At time $t=0$, $b(0)=C$ and $b(t)\to 1$ when $t\gg(\mu+\sigma)$.

Figure 3: Simulation results using rate operator $\mathbf{R1}$ [(a), (b)] and
$\mathbf{R1^{\prime}}$ [(c), (d)] for the model (44) with time dependent
driving. Example trajectories are displayed in (a) and (c). Final steady state
distributions for Bloch vector X and Y components are shown in (b) and (d).
The insets in (b) and (d) demonstrate the agreement between the analytical
results (solid line) and simulation results (dots).

The driving does not affect the positivity of the rate operators $\mathbf{R1}$
and $\mathbf{R2}$. However, we use the decomposition (2) and define a new rate
operator

$\displaystyle\mathbf{R1^{\prime}}_{\psi(t)}=\mathbf{R1}_{\psi(t)}+i\frac{b(t)}{2}[\sigma_{z},|{\psi(t)}\rangle\langle{\psi(t)}|],$
(50)

which fully takes into account the driving, so that the deterministic
evolution between the jumps will not depend on it. In other terms, as
announced earlier, we can implement the simulation of time-dependent coherent
driving with pure jump dynamics (provided that the rate operator
$\mathbf{R1^{\prime}}$ is positive).

In Fig. 3 we plot the dynamics for $\mu=1$, $\sigma=1/4$ and for an initial
state $\ket{\psi_{0}}=\cos\frac{\pi}{8}|1\rangle+\sin\frac{\pi}{8}|2\rangle$.
In panel (a) the coherent driving is clearly visible in the deterministic
evolution between the jumps when using the rate operator $\mathbf{R1}$. In
contrast, in panel (c) the sample trajectories do not have any deterministic
evolution between the jumps since the driving is absorbed by the rate operator
$\mathbf{R1^{\prime}}$. In the insets of panels (b) and (d), we verify that we
indeed unravel the dynamics of Eq. (44). The distribution of the states in the
long time limit can be understood by looking at
$\varphi(t)=\arctan\frac{y(t)}{x(t)}$, i.e., the angle between the $x$\- and
$y$-axis of the Bloch sphere, measured from the positive $x$-axis. In the long
time limit, when $b(t)\approx 1$ and all of the trajectories have almost
reached the equator of the Bloch sphere, each jump just shifts the phase
$\varphi(t_{+})=\varphi(t_{-})\pm\frac{\pi}{4}$, where times $t_{\pm}$ are
just after and before the jump, respectively. The initial transient period
randomizes the phases and the resulting smeared distribution for the phase, as
seen in panel (d) of Fig. 3, occurs.

For further confirmation, we plot in Figure 4 $\varphi(t)$ for a single
trajectory.

Figure 4:  The phase evolution (upper panel) and phase jumps (lower panel) for
a single trajectory with time dependent driving. We can see that initially
there is a transient where the phase jumps are not multiples of $\pi/4$. This
leads to the smearing of the ensemble in the long time limit observed in Fig.
3

.

Figure 5:  Time independent driving for $b(t)=1$. In panel a) we plot the
dynamics using rate operator $\mathbf{R1}$ and the driving term is clearly
visible in the real part of the coherences. In panel b) we verify in the inset
the validity with the simulations and plot the distribution of Bloch vectors
in the ensemble at at the final time. In paneld c) we use the rate operator
$\mathbf{R1^{\prime}}$ and the trajectories contain only jumps. In panel d) we
plot the Bloch vectors of the ensemble at the final time. Clearly, the phase
is rational multiple of $\pi/4$ for every member of the ensemble.

We clearly see that after a transient period the phase changes during the
jumps in steps $\frac{\pi}{4}$, which can be understood as follows. Let
$z^{2}=\cos^{2}\theta$ be the squared $z$-component of a unit Bloch vector in
spherical coordinates. The rate operators $\mathbf{R1}$ and
$\mathbf{R1^{\prime}}$ both induce a map
$z^{2}\mapsto\frac{\alpha^{2}z^{2}}{(2-\alpha)^{2}}$ when $t\to\infty$. Hence,
the jumps push the dynamics on the equatorial plane of the Bloch sphere. Since
the coherent driving preserves the $z$-component, the pure states of the
ensemble will be distributed along the equator of the Bloch sphere in the long
time limit. Any such pure state can be written as
$|\psi\rangle=\frac{1}{\sqrt{2}}(|1\rangle+e^{i\varphi}|2\rangle)$. The
eigenstates of the rate operator $\mathbf{R1^{\prime}}_{\psi}$ in the long
time limit are $|\phi_{\psi,\pm}\rangle=\frac{1}{2}(|1\rangle\pm
e^{-i\varphi}\xi|2\rangle)$, where $\xi=\frac{1+ib}{\sqrt{1+b^{2}}}$.
Importantly, $|\xi|=1$ and the relative phase of the eigenstates is
$e^{i\varphi}\xi$. In the long time limit, the effect of the quantum jumps
corresponds to phase jumps $e^{i\phi}\mapsto\pm e^{i\phi}\xi$. In the long
time limit $b(t)=1$. Interestingly, the phase jumps for the rate operator
$\mathbf{R1^{\prime}}$ in the special case $b=1$ are rational multiples of
$\pi$ in the long time limit. Namely, $\xi=e^{i\pi/4}$ which is the 8th root
of unity.

In Fig. 5 we plot the dynamics for time-independent driving $b(t)=1$ and
especially in panel d) we can see how in this case the definite phase
relations are preserved. The phase $\varphi(t)$ aquires values which are
rational multiples of $\pi/4$ for every trajectory of the ensemble.

##  6 Conclusions

We introduced a class of quantum-jump unravelings (denoted as
$\mathbf{R}$-ROQJ) that interpolates between MCWF for CP-divisible dynamics
and the rate-operator unraveling for P-divisible evolutions [40, 30].
$\mathbf{R}$-ROQJ guarantees a well-defined continuous-measurement
interpretation for the class of quantum evolutions defined by the condition
(42), that is strictly stronger than P-divisibility. Furthermore,
$\mathbf{R}$-ROQJ takes full advantage of the freedom we have in dividing the
master equation into a determinisitic and a jump part, by setting a desired
time-independent linear non-Hermitian Hamiltonian, even in the case of an
external driving, or selecting a fixed post-jump basis. In addition to
fundamental insight to open system dynamics, our results also allow us to
simplify numerical simulations for practical purposes. The approach put
forward here will likely be useful also to deal with more general classes of
non-Markovian dynamics, possibly allowing for a greater control of the
reversed jumps [25, 26] via the manipulation of the deterministic and jump
parts of the unraveling.

D.C. was supported by the Polish National Science Centre project
2018/30/A/ST2/00837. J.P. acknowledges support from Magnus Ehrnrooth
Foundation. A.S. acknowledges funding by the FFABR project of MIUR and PSR-2
2020 by UNIMI.

## References

  * [1] H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford Univ. Press, Oxford, 2007).
  * [2] H.J. Carmichael, _An Open System Approach to Quantum Optics, Lectures Notes in Physics_ (Springer, Berlin, 1993)
  * [3] J. Dalibard, Y. Castin, and K. Mølmer, Phys. Rev. Lett. 68, 580 (1992).
  * [4] T. Basche, S. Kummer, and C. Brauchle, Nature 373, 132 (1995).
  * [5] S. Peil and G. Gabrielse, Phys. Rev. Lett. 83, 1287 (1999).
  * [6] F. Jelezko, I. Popa, A. Gruber, C. Tietz, J. Wrachtrup, A. Nizovtsev, and S. Kilin, Appl. Phys. Lett. 81, 2160 (2002).
  * [7] S. Gleyzes, S. Kuhr, C. Guerlin, J. Bernu, S. Deléglise, U.B. Hoff, M. Brune, J.-M. Raimond, and S. Haroche, Nature 446, 297 (2007).
  * [8] R. Vijay, D. H. Slichter, and I. Siddiqi, Phys. Rev. Lett. 106, 110502 (2011).
  * [9] 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, Nature 570, 200 (2019).
  * [10] M. B. Plenio and P. L. Knight, Rev. Mod. Phys. 70, 101 (1998).
  * [11] A.J. Daley, Adv. Phys. 63, 77 (2014).
  * [12] I.Percival, _Quantum State Diffusion_ (Cambridge University Press, Cambridge, England, 2002)
  * [13] A. Barchielli and M. Gregoratti, _Quantum Trajectories and Measurements in Continuous Time: The Diffusive Case_ , Lecture Notes in Physics 782 (Springer, Berlin, 2009).
  * [14] W. T. Strunz, L. Diósi, and N. Gisin, Phys. Rev. Lett. 82, 1801 (1999).
  * [15] T. Yu, L. Diósi, N. Gisin, and W. T. Strunz, Phys. Rev. A 60, 91 (1999).
  * [16] K. Luoma, W.T. Strunz, and J. Piilo, Phys. Rev. Lett. 125, 150403 (2020)
  * [17] A. Barchielli and V.P. Belavkin, J. Phys. A: Math. Gen. 24 1495 (1991).
  * [18] E.-M. Laine, J. Piilo, and H.-P. Breuer, Phys. Rev. A 81 062115 (2010).
  * [19] D. Chrusciński, A. Kossakowski, and Á. Rivas, Phys. Rev. A 83, 052128 (2011).
  * [20] Á. Rivas and S. F. Huelga, _Open Quantum Systems_ (Springer, New York, 2012).
  * [21] Á. Rivas, S. F. Huelga, and M. B. Plenio, Phys. Rev. Lett. 105, 050403 (2010).
  * [22] Á. Rivas, S. F. Huelga, and M. B. Plenio, Rep. Prog. Phys. 77, 094001 (2014).
  * [23] H.-P. Breuer, E.-M. Laine, and J. Piilo, Phys. Rev. Lett. 103, 210401 (2009).
  * [24] H.-P. Breuer, E.-M. Laine, J. Piilo, and B. Vacchini, Rev. Mod. Phys. 88, 021002 (2016).
  * [25] J. Piilo, S. Maniscalco, K. Härkönen, and K.-A. Suominen, Phys. Rev. Lett. 100,180402 (2008).
  * [26] J. Piilo, K. Härkönen, S. Maniscalco, and K.-A. Suominen, Phys. Rev. A 79, 062112 (2009).
  * [27] J. Gambetta and H.M. Wiseman, Phys. Rev. A 68, 062104 (2003).
  * [28] L. Diósi, Phys. Rev. Lett. 100, 080401 (2008).
  * [29] H.M. Wiseman and J.M. Gambetta, Phys. Rev. Lett. 101, 140401 (2008).
  * [30] A. Smirne, M. Caiaffa, and J. Piilo, Phys. Rev. Lett. 124, 190402 (2020).
  * [31] L. Diósi, Phys. Lett. A 114, 451 (1986).
  * [32] L. Diósi, J. Phys. A 21, 2885 (1988).
  * [33] N. Gisin, Helv. Phys. Acta 63 929 (1990).
  * [34] B. Vacchini, A. Smirne, E.-M. Laine, J. Piilo, H.-P. Breuer, New J. Phys. 13, 093004 (2011).
  * [35] D. Chruściński and S. Maniscalco, Phys. Rev. Lett. 112, 120404 (2014).
  * [36] S. Wißmann, H.-P. Breuer, B. Vacchini, Phys. Rev. A 92, 042108 (2015).
  * [37] G. Lindblad, Comm. Math. Phys. 48, 119 (1976).
  * [38] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, J. Math. Phys. 17, 821 (1976).
  * [39] D. Chrusciński, and A. Kossakowski, Phys. Rev. Lett. 104, 070406 (2010).
  * [40] M. Caiaffa, A. Smirne, and A. Bassi, Phys. Rev. A 95 062101 (2017).
  * [41] T. Heinosaari and M. Ziman, _The Mathematical Language of Quantum Theory_ , (Cambridge University Press, Cambridge, 2012).
  * [42] M.J.W. Hall, J.D. Cresser, L. Li, and E. Andersson, Phys. Rev. A 89, 042120 (2014).
  * [43] D. Chruściński and F.A. Wudarski, Phys.Rev. A 91, 012104 (2015).
  * [44] N. Megier, D. Chruscinski, J. Piilo, and W.T. Strunz, Sci. Rep. 7, 6379 (2017).
  * [45] V. Paulsen, Completely Bounded Maps and Operator Algebras (Cambridge University Press, Cambridge, 2003).
  * [46] E. Størmer, Positive Linear Maps of Operator Algebras, Springer Monographs in Mathematics (Springer, New York, 2013).
  * [47] D. Chruściński and F. Mukhamedov, Phys. Rev. A. 100, 052120 (2019).

##  Appendix A Proof of Proposition 1

Consider $\mathcal{H}_{S}=\mathbbm{C}^{2}$ and a master equation as in Eq.(1);
moreover assume that there is a basis
$\left\\{\ket{\varphi_{1}},\ket{\varphi_{2}}\right\\}$ and a linear operator
$\mathbf{C}(t)$ on $\mathbbm{C}^{2}$ such that Eqs.(32) and (38) hold. First
note that Eq.(32), along with the Hermiticity preservation condition, imply
(we neglect from now on the time dependence) $J^{kl}_{ij}=J^{lk}_{ji}$. Using
the latter, it is readily shown that the action of $\mathcal{J}^{\prime}_{t}$
in Eq.(2) with $\mathbf{C}$ defined as in Eq.(33) can be written as

$\mathcal{J}^{\prime}_{t}(\overline{\rho})=\mathcal{J}_{t}(\overline{\rho})+\frac{1}{2}\left(\mathbf{C}\overline{\rho}+\overline{\rho}\mathbf{C}^{\dagger}\right)=\left(\begin{array}[]{cc}\alpha&\beta\\\
\beta&\alpha\end{array}\right)\quad\begin{cases}\alpha=&(2J^{12}_{11}+y)\rho_{12}+J^{22}_{11}\rho_{22}+J^{11}_{22}\rho_{11}\\\
\beta=&(J^{11}_{12}+\frac{y}{2}+J^{12}_{11}+J^{12}_{22})\rho_{11}+(J^{22}_{12}+\frac{y}{2})\rho_{22}\\\
&+\frac{1}{2}(2J^{12}_{12}+2J^{21}_{12}+J^{11}_{22}-J^{11}_{11}+J^{22}_{11}-J^{22}_{22})\rho_{12}\end{cases}$
(A.1)

for any state $\overline{\rho}$ such that

$\rho_{12}=\rho_{21},$ (A.2)

where indeed $\rho_{ij}=\bra{\varphi_{i}}\overline{\rho}\ket{\varphi_{j}}$
(i.e., there are no phases with respect to the selected basis).

Now, the crucial observation is that the matrix in Eq.(A.1) has _fixed_
eigenvectors: for any values of $\alpha$ and $\beta$, they are in fact given
by

$\ket{\varphi_{\pm}}=\frac{1}{\sqrt{2}}\left(\ket{\varphi_{1}}\pm\ket{\varphi_{2}}\right),$
(A.3)

which directly implies the statement to be proven, whenever the rate operator
is positive if referred to pure states satisfying Eq.(A.2):

  * •

as said, we choose an initial state $\ket{\psi_{0}}$ that satisfies Eq.(37);

  * •

but then, due to Eq.(38), the state $\ket{\psi_{0}(t_{1})}$ before the first
jump at time $t_{1}$ (recall that $D(t,s)$ is defined in Eq.(39))

$\frac{D(t_{1},0)\ket{\psi_{0}}}{\left\|D(t_{1},0)\ket{\psi_{0}}\right\|}$

will be such that Eq.(A.2) holds, so that the jump operator
$\mathbf{R}_{\psi_{0}(t_{1})}=\mathcal{J}^{\prime}_{t}(\ket{\psi_{0}(t_{1})}\bra{\psi_{0}(t_{1})})$
fixed by Eq.(2) with $\mathbf{C}(t)$ as in Eq.(33) has the two eigenvectors
$\ket{\varphi_{\pm}}$ in Eq.(A.3);

  * •

thus, after the jump (that is well defined due to the positivity of the rate
operator) the system will be either in $\ket{\varphi_{+}}$ or in
$\ket{\varphi_{-}}$;

  * •

since also the latter are in the form as in Eq.(37) the state before the
second jump,

$\frac{D(t_{2},t_{1})\ket{\varphi_{\pm}}}{\left\|D(t_{2},t_{1})\ket{\varphi_{\pm}}\right\|},$

will still satisfy Eq.(A.2), so that the state after the second jump will be
either $\ket{\varphi_{+}}$ or $\ket{\varphi_{-}}$, and so on.

All in all, we have the 3 possible families of states
$\left\\{\ket{\psi_{0}(t)},\ket{\varphi_{+}(t,s)},\ket{\varphi_{-}(t,s)}\right\\}$
that are obtained from
$\left\\{\ket{\psi_{0}},\ket{\varphi_{+}},\ket{\varphi_{-})}\right\\}$ via the
deterministic evolution in Eq.(39); the instants of occurrence of the jumps
will determine which parts of the deterministic evolution will be actually
involved in each trajectory.

##  Appendix B Proof of Proposition 3

Defining

$\mathbf{K}(t)=\int_{U(N)}\mathcal{L}_{t}^{\dagger}(U^{\dagger})UdU,$ (B.1)

one finds for an arbitrary $V\in U(N)$

$\int_{U(N)}\mathcal{L}_{t}^{\dagger}(VU^{\dagger})UdU=\int_{U(N)}\mathcal{L}_{t}^{\dagger}({U^{\prime}}^{\dagger})U^{\prime}VdU^{\prime}=\mathbf{K}(t)V,$
(B.2)

with $U^{\prime}=UV^{\dagger}$. Hence for an arbitrary system operator $X$ one
has

$\int_{U(N)}\mathcal{L}_{t}^{\dagger}(XU^{\dagger})UdU=\mathbf{K}(t)X.$ (B.3)

Now we use the dissipativity condition, Eq.(42),

$\mathcal{L}^{\dagger}_{t}(Y^{\dagger}Y)\geq\mathcal{L}^{\dagger}_{t}(Y^{\dagger})Y+Y^{\dagger}\mathcal{L}^{\dagger}_{t}(Y),$
(B.4)

for $Y=UX$, with $U\in U(N)$. It leads to

$\displaystyle\mathcal{L}^{\dagger}_{t}(Y^{\dagger}Y)\geq\mathcal{L}^{\dagger}_{t}(X^{\dagger}U^{\dagger})UX+X^{\dagger}U^{\dagger}\mathcal{L}^{\dagger}_{t}(UX)=\mathcal{L}^{\dagger}_{t}(X^{\dagger}U^{\dagger})UX+X^{\dagger}[\mathcal{L}^{\dagger}_{t}(X^{\dagger}U^{\dagger})U]^{\dagger}.$
(B.5)

Averaging over $U(N)$ yields

$\displaystyle\int_{U(N)}\mathcal{L}^{\dagger}_{t}(Y^{\dagger}Y)dU\geq\int_{U(N)}\mathcal{L}^{\dagger}_{t}(X^{\dagger}U^{\dagger})UdU\,X++X^{\dagger}\left(\int_{U(N)}\mathcal{L}^{\dagger}_{t}(X^{\dagger}U^{\dagger})UDU\right)^{\dagger}=\mathbf{K}(t)X^{\dagger}X+X^{\dagger}X\mathbf{K}^{\dagger}(t),$
(B.6)

and hence defining a linear map as

$\mathbf{J}^{\dagger}_{t}(X)=\mathcal{L}_{t}^{\dagger}(X)-(\mathbf{K}(t)X+X\mathbf{K}^{\dagger}(t)),$
(B.7)

we have shown that

$\mathbf{J}^{\dagger}_{t}(X^{\dagger}X)\geq 0,$ (B.8)

which proves that the map $\mathbf{J}_{t}$ is positive.

##  Appendix C Positivity of the rate operators $\mathbf{R1}_{\psi(t)}$,
$\mathbf{R2}_{\psi(t)}$ and $\mathbf{R3}_{\psi(t)}$

Here, we prove the positivity of $\mathbf{R1}_{\psi(t)}$,
$\mathbf{R2}_{\psi(t)}$, and $\mathbf{R3}_{\psi(t)}$ defined in Sec.5, fixed
by Eqs.(46), (47) and (48), assuming only that the $\gamma_{k}(t)$ satisfy
P-divisibility condition, that is, $\gamma_{i}(t)+\gamma_{j}(t)\geq 0$ for
$i\neq j=1,2,3$ and that $\gamma_{3}(t)<0$, so that P-divisibility means that

$\gamma_{k}(t)+\gamma_{3}(t)\geq 0\ ,\ k=1,2.$

Note that the eternal non-Markovian dynamics considered in the main text is a
special case of this class of dynamics.

First, we note that the qubit generator

$\mathcal{L}_{t}(\rho)=\frac{1}{2}\sum_{k=1}^{3}\gamma_{k}(t)(\sigma_{k}\rho\sigma_{k}-\rho),$
(C.1)

may be rewritten as

$\mathcal{L}_{t}(\rho)=\sum_{k=1}^{3}\gamma_{k}(t)(J_{k}(\rho)-\rho),$ (C.2)

where

$J_{k}(\rho)=\frac{1}{2}(\sigma_{k}\rho\sigma_{k}+\rho),$ (C.3)

are CPTP for $k=1,2,3$. We then have the following result.

######  Proposition 4.

The map (cf. Eq. (46))

$\mathbf{J}_{t}(\rho)=\sum_{k=1}^{3}\gamma_{k}(t)J_{k}(\rho)=\mathcal{J}_{t}(\rho)+\frac{\gamma(t)}{2}\rho$
(C.4)

is positive.

The proof is based on the following observation: for a qubit system a linear
map $\Phi:M_{2}(\mathbb{C})\to M_{2}(\mathbb{C})$ is positive if and only the
corresponding Choi matrix $C_{\Phi}$

$C_{\Phi}:=\sum_{i,j=1}^{2}|i\rangle\langle j|\otimes\Phi(|i\rangle\langle
j|)$

satisfies

$C_{\Phi}=A+(\mathbbm{1}\otimes T)B,$

where $A$ and $B$ are positive $4\times 4$ matrices, and ‘$\mathbbm{1}\otimes
T$’ denotes partial transposition. The Choi matrix for the map
$\mathbf{J}_{t}$ reads

$C_{\mathbf{J}_{t}}=\left(\begin{array}[]{cc|cc}{\gamma_{1}}+{\gamma_{2}}+2{\gamma_{3}}&0&0&{\gamma_{1}}+{\gamma_{2}}\\\
0&{\gamma_{1}}+{\gamma_{2}}&{\gamma_{1}}-{\gamma_{2}}&0\\\ \hline\cr
0&{\gamma_{1}}-{\gamma_{2}}&{\gamma_{1}}+{\gamma_{2}}&0\\\
{\gamma_{1}}+{\gamma_{2}}&0&0&{\gamma_{1}}+{\gamma_{2}}+2{\gamma_{3}}\\\
\end{array}\right),$ (C.5)

where we skipped the time dependence. Note, that
$C_{\mathbf{J}_{t}}=A+(\mathbbm{1}\otimes T)B$, where

$A=\left(\begin{array}[]{cc|cc}{\gamma_{1}}+{\gamma_{2}}+2{\gamma_{3}}&0&0&{\gamma_{1}}+{\gamma_{2}}+2\gamma_{3}\\\
0&{\gamma_{1}}+{\gamma_{2}}+2\gamma_{3}&{\gamma_{1}}-{\gamma_{2}}&0\\\
\hline\cr
0&{\gamma_{1}}-{\gamma_{2}}&{\gamma_{1}}+{\gamma_{2}}+2\gamma_{3}&0\\\
{\gamma_{1}}+{\gamma_{2}}+2\gamma_{3}&0&0&{\gamma_{1}}+{\gamma_{2}}+2{\gamma_{3}}\\\
\end{array}\right),$ (C.6)

and

$B=-2\gamma_{3}\left(\begin{array}[]{cc|cc}0&0&0&0\\\ 0&1&1&0\\\ \hline\cr
0&1&1&0\\\ 0&0&0&0\\\ \end{array}\right).$ (C.7)

Since $\gamma_{3}(t)<0$ one has $B\geq 0$. Now, $A\geq 0$ if and only if

${\gamma_{1}}+{\gamma_{2}}+2\gamma_{3}\geq|{\gamma_{1}}-{\gamma_{2}}|.$ (C.8)

Let us consider two scenario:

  1. 1.

if $\gamma_{1}\geq\gamma_{2}$, then (C.8) is equivalent to
$\gamma_{2}+\gamma_{3}\geq 0$,

  2. 2.

conversely, if $\gamma_{2}\geq\gamma_{1}$, then (C.8) is equivalent to
$\gamma_{1}+\gamma_{3}\geq 0$,

which ends the proof of positivity of $A$. Now, since $A$ and $B$ are
positive, the map $\mathbf{J}_{t}$ is positive and hence
$\mathbf{R1}_{\psi(t)}\geq 0$.

Positivity of $\mathbf{R2}_{\psi(t)}$ is evident due to the relation

$\mathbf{R2}_{\psi(t)}=\mathbf{R1}_{\psi(t)}-\gamma_{3}(t)|\psi(t)\rangle\langle\psi(t)|,$
(C.9)

and $\gamma_{3}(t)<0$. Finally, the operator $\mathbf{R3}_{\psi(t)}$
corresponds to the map
$\frac{1}{2}(\sum_{k}\gamma_{k}\sigma_{k}\rho\sigma_{k}-\gamma_{3}\rho)$. Its
Choi matrix reads

$C=\left(\begin{array}[]{cc|cc}0&0&0&-2{\gamma_{3}}\\\
0&{\gamma_{1}}+{\gamma_{2}}&{\gamma_{1}}-{\gamma_{2}}&0\\\ \hline\cr
0&{\gamma_{1}}-{\gamma_{2}}&{\gamma_{1}}+{\gamma_{2}}&0\\\
-2{\gamma_{3}}&0&0&0\\\ \end{array}\right)=A^{\prime}+(\mathbbm{1}\otimes
T)B^{\prime},$ (C.10)

where

$A^{\prime}=\left(\begin{array}[]{cc|cc}0&0&0&0\\\
0&{\gamma_{1}}+{\gamma_{2}}+2\gamma_{3}&{\gamma_{1}}-{\gamma_{2}}&0\\\
\hline\cr
0&{\gamma_{1}}-{\gamma_{2}}&{\gamma_{1}}+{\gamma_{2}}+2\gamma_{3}&0\\\
0&0&0&0\\\ \end{array}\right),$ (C.11)

and $B^{\prime}=B$ defined by (C.7). The proof of positivity of $A^{\prime}$
is very similar to that of $A$.

##  Appendix D Fixed post-jump states for the eternal non-Markovian dynamics

Here, we first show explicitly that the dynamics defined by the master
equation (44) is actually a special case of the two-level system dynamics
identified by Prop.1, so that we can define a $\mathbf{R}$-ROQJ unraveling
with fixed post-jump states via Eq.(33). We will then discuss the positivity
of such an unraveling.

Consider the basis $\left\\{\ket{1},\ket{0}\right\\}$ of eigenvectors of
$\sigma_{z}$; in such a basis the coefficients of $\mathcal{J}_{t}$ according
to Eq.(31) read

$\left(\begin{array}[]{c|cc}J^{11}_{11}(t)&\gamma_{3}(t)/2\\\
J^{11}_{10}(t)&0\\\ J^{11}_{01}(t)&0\\\
J^{11}_{00}(t)&(\gamma_{1}(t)+\gamma_{2}(t))/2\\\ J^{10}_{11}(t)&0\\\
J^{10}_{10}(t)&-\gamma_{3}(t)/2\\\
J^{10}_{01}(t)&(\gamma_{1}(t)-\gamma_{2}(t))/2\\\ J^{10}_{00}(t)&0\\\
J^{01}_{11}(t)&0\\\ J^{01}_{10}(t)&(\gamma_{1}(t)-\gamma_{2}(t))/2\\\
J^{01}_{01}(t)&-\gamma_{3}(t)/2\\\ J^{01}_{00}(t)&0\\\
J^{00}_{11}(t)&(\gamma_{1}(t)+\gamma_{2}(t))/2\\\ J^{00}_{10}(t)&0\\\
J^{00}_{01}(t)&0\\\ J^{00}_{00}(t)&\gamma_{3}(t)/2\end{array}\right)$ (D.1)

and hence (setting $x=0$) the operator $\mathbf{C}(t)$ as in Eq.(33) is

$\displaystyle\mathbf{C}(t)$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}\frac{\gamma_{1}(t)+\gamma_{2}(t)-\gamma_{3}(t)}{2}&y(t)\\\
y(t)&\frac{\gamma_{1}(t)+\gamma_{2}(t)-\gamma_{3}(t)}{2}\end{array}\right)=\frac{\gamma_{1}(t)+\gamma_{2}(t)-\gamma_{3}(t)}{2}\mathbbm{1}+y(t)\sigma_{x};$
(D.4)

while the operator fixing the deterministic evolution according to Eq.(39) is

$D(t,s)=\exp\left(-\frac{1}{2}\int_{s}^{t}d\tau(\gamma_{1}(\tau)+\gamma_{2}(\tau))\mathbbm{1}+y(\tau)\sigma_{x}\right);$
(D.5)

note that for $y(t)=0$ we recover the rate operator $\mathbf{R2}_{\psi(t)}$
defined by Eq.(47). It is easy to see that the assumptions of Prop.1 in
Eq.(32) and (38) hold, so that, starting from an initial state as in Eq.(37),
the $\mathbf{R}$-ROQJ fixed by $\mathbf{C}(t)$ as in Eq.(D.4) consists of
jumps among 3 deterministically evolving states, whenever the positivity of
the rate operator is ensured. In particular, for a continuous family of
$\mathbf{R}$-ROQJ this is the case whenever the dynamics is P-divisible, as
shown below.

First, recall that the P-divisibility condition for the eternal non-Markovian
dynamics is $\gamma_{1}(t),\gamma_{2}(t)\geq|\gamma_{3}(t)|$ for all $t\geq
0$. Moreover, going back to the proof of Prop. 1 in Sec.B, the positivity of
the rate operator is fixed by the eigenvalues

$\alpha(t)\pm\beta(t)$

of $\mathcal{J}^{\prime}_{t}(\ket{\psi(t)}\bra{\psi(t)})$, see the first
equality in Eq.(A.1). Actually, rather than studying the eigenvalues, it is
convenient to look at the trace and determinant of the corresponding $2\times
2$ matrix: for any pure state of the form
$\ket{\psi}=c\ket{\varphi_{1}}\pm\sqrt{1-c^{2}}\ket{\varphi_{2}}$ (i.e., as in
Eq.(37)) we then have the two conditions

$\displaystyle\alpha(t)$ $\displaystyle=$
$\displaystyle\frac{\gamma_{1}(t)+\gamma_{2}(t)}{2}\pm c\sqrt{1-c^{2}}y(t)\geq
0$ $\displaystyle\alpha^{2}(t)-\beta^{2}(t)$ $\displaystyle=$
$\displaystyle\left(\frac{\gamma_{1}(t)+\gamma_{2}(t)}{2}\pm
c\sqrt{1-c^{2}}y(t)\right)^{2}-\left(\frac{y(t)}{2}\pm
c\sqrt{1-c^{2}}(\gamma_{1}(t)-\gamma_{3}(t))\right)^{2}\geq 0.$ (D.6)

Both inequalities hold for any $-1\leq c\leq 1$ if $0\leq
y(t)\leq\gamma_{1}(t)+\frac{1}{2}\gamma_{2}(t)-\frac{1}{2}\gamma_{3}(t)$ for
all $t\geq 0$. The validity of the first inequality directly follows from the
fact that $\gamma_{1}(t),\gamma_{2}(t)\geq 0$ and
$\gamma_{1}(t)+\frac{1}{2}\gamma_{2}(t)-\frac{1}{2}\gamma_{3}(t)\leq\gamma_{1}(t)+\gamma_{2}(t)$
due to P-divisibility and that the maximum value of $|c|\sqrt{1-c^{2}}$ for
$-1\leq c\leq 1$ is $1/2$. For the second inequality, let us consider first
the case where the term $\pm c$ appearing in it is positive, so that the
inequality becomes

$\gamma_{1}(t)+\gamma_{2}(t)-y(t)\geq
2|c|\sqrt{1-c^{2}}\left(\gamma_{1}(t)-\gamma_{3}(t)-y(t)\right)$

as $\gamma_{1}(t)-\gamma_{3}(t)\geq 0$ due to P-divisibility and we used that
$\alpha(t)\geq 0$, as well as $y(t)\geq 0$; but the maximum value with respect
to $c$ of the term at the r.h.s. is
$\max\left\\{0,\gamma_{1}(t)-\gamma_{3}(t)-y(t)\right\\}$, so that the
inequality holds for any $y(t)\leq\gamma_{1}(t)+\gamma_{2}(t)$
($\gamma_{2}(t)\geq-\gamma_{3}(t)$ due to P-divisibility). Instead, if $\pm c$
is negative the second inequality in (D.6) becomes

$\frac{\gamma_{1}(t)+\gamma_{2}(t)}{2}-|c|\sqrt{1-c^{2}}y(t)\geq\left|\frac{y(t)}{2}-|c|\sqrt{1-c^{2}}(\gamma_{1}(t)-\gamma_{3}(t))\right|.$

If $y(t)\geq 2|c|\sqrt{1-c^{2}}(\gamma_{1}(t)-\gamma_{3}(t))$, we have
$\gamma_{1}(t)+\gamma_{2}(t)-y(t)\geq
2|c|\sqrt{1-c^{2}}(y(t)-(\gamma_{1}(t)-\gamma_{3}(t))),$ whose r.h.s. maximum
value with respect to $c$ is
$\max\left\\{0,y(t)-(\gamma_{1}(t)-\gamma_{3}(t))\right\\}$ so that the
inequality holds for any $-1\leq c\leq 1$ if
$y(t)\leq\gamma_{1}(t)+\frac{1}{2}\gamma_{2}(t)-\frac{1}{2}\gamma_{3}(t)$. If
$y(t)<2|c|\sqrt{1-c^{2}}(\gamma_{1}(t)-\gamma_{3}(t))$, we have
$\gamma_{1}(t)+\gamma_{2}(t)+y(t)\geq
2|c|\sqrt{1-c^{2}}(y(t)+(\gamma_{1}(t)-\gamma_{3}(t))),$ whose r.h.s. maximum
is $\max\left\\{0,y(t)+(\gamma_{1}(t)-\gamma_{3}(t))\right\\}$ so that the
inequality is ensured by $\gamma_{2}(t)\geq-\gamma_{3}(t)$, as well as
$\gamma_{1}(t),\gamma_{2}(t),y(t)\geq 0$.

In total, for any $y(t)$ such that $0\leq
y(t)\leq\gamma_{1}(t)+\frac{1}{2}\gamma_{2}(t)-\frac{1}{2}\gamma_{3}(t)$ for
all $t\geq 0$ the rate operator
$\mathbf{R}_{\psi}=\mathcal{J}^{\prime}_{t}(\ket{\psi}\bra{\psi})$ is positive
for any $\ket{\psi}=c\ket{\varphi_{1}}\pm\sqrt{1-c^{2}}\ket{\varphi_{2}}$ with
$-1\leq c\leq 1$ and then if the initial state satisfies the latter condition
(which is preserved by the deterministic evolution in Eq.(D.5)), we will have
a _positive_ rate-operator unravelling with the states
$\left\\{\ket{\psi_{0}(t)},\ket{\varphi_{+}(t,s)},\ket{\varphi_{-}(t,s)}\right\\}$.