# Stochastic entropy production for continuous measurements of an open quantum
system

D. Matos  Department of Physics, King’s College London, Strand London, WC2R
2LS, United Kingdom L. Kantorovich  Department of Physics, King’s College
London, Strand London, WC2R 2LS, United Kingdom I.J. Ford  Department of
Physics and Astronomy, University College London, Gower Street, London, WC1E
6BT, United Kingdom

###### Abstract

We investigate the total stochastic entropy production of a two-level bosonic
open quantum system under protocols of time dependent coupling to a harmonic
environment. These processes are intended to represent the measurement of a
system observable, and consequent selection of an eigenstate, whilst the
system is also subjected to thermalising environmental noise. The entropy
production depends on the evolution of the system variables and their
probability density function, and is expressed through system and
environmental contributions. The continuous stochastic dynamics of the open
system is based on the Markovian approximation to the exact, noise-averaged
stochastic Liouville-von Neumann equation, unravelled through the addition of
stochastic environmental disturbance mimicking a measuring device. Under the
thermalising influence of time independent coupling to the environment, the
mean rate of entropy production vanishes asymptotically, indicating
equilibrium. In contrast, a positive mean production of entropy as the system
responds to time dependent coupling characterises the irreversibility of
quantum measurement, and a comparison of its production for two coupling
protocols, representing connection to and disconnection from the external
measuring device, satisfies a detailed fluctuation theorem.

##  1 Introduction

Most applied and theoretical quantum mechanics research is underpinned by the
theory of open quantum systems. Any realistic quantum system interacts with
its environment, though the interaction in many cases is weak. Open quantum
systems are of interest due to the characteristic dynamical properties they
display, specifically irreversible dissipative behaviour on approach to a
steady state, and decoherence, both of which are not found in closed systems.
These phenomena are of crucial importance in quantum technological
applications such as quantum computing [1] and in theoretical developments in
quantum thermodynamics [2].

In this paper, we concern ourselves with the thermodynamic behaviour of an
open quantum system, more specifically with the stochastic entropy production
associated with the evolution of its reduced density matrix brought about by
changes in interactions with the environment. Such a framework can naturally
describe the consequences arising from a time dependent Hamiltonian coupling
of the system to its environment, but it can also represent the effect of
quantum measurements involving parts of the environment, such as a measuring
device. The field of stochastic thermodynamics began as a generalisation of
the laws of thermodynamics applied to stochastic systems [3], such as the
behaviour of colloidal particles or molecular systems exposed to heat baths
[4]. These developments allow us to compute the entropy production associated
with individual stochastic trajectories of the evolving reduced density
matrix. This stochastic entropy production can satisfy a detailed fluctuation
theorem describing the relationship between the effects of time-reversed
versions of the coupling protocol [5].

In order to apply the tools of stochastic thermodynamics to the quantum
regime, a notion of quantum trajectories must be established. Quantum
trajectories have long been used in the field of quantum optics [6], and they
are typically realised through the unravelling of a deterministic equation of
motion for the variables describing the system. Unravelling is the elaboration
of the deterministic dynamics into a stochastic equation of motion that is
capable of describing randomness in system behaviour, such as fluctuations
around the average evolution. The deterministic dynamics implicitly describes
an average over the range of random behaviour. Randomness might be ascribed to
outcomes of continuous measurements involving the environment [7], such as
photon counting or homodyne/heterodyne detection, though unitary evolution of
a system together with its environment where the state of the environment is
not fully specified can also be represented using a framework of stochastic
unravelling.

To obtain these unravellings, we start with the stochastic Liouville-von
Neumann (SLN) equation for the system’s reduced density matrix [8, 9, 10, 11,
12], a non-Markovian stochastic differential equation derived from the
Feynman-Vernon path integral-based consideration of open quantum system
behaviour [13]. It is important to note that the stochasticity in the SLN is
merely a mathematical feature that requires averaging in order to produce the
exact deterministic behaviour of the reduced density matrix. We then derive
the Markovian limit of the noise-averaged SLN equation in order to make
possible a simple unravelling of the deterministic dynamics using a Kraus
operator representation of physical stochasticity brought about by the
environment. The stochastic differential equation that emerges then describes
physical randomness and can provide a basis for deriving the stochastic
entropy production.

We employ stochastic entropy production to quantify the irreversibility of the
continuous measurement of the quantum state of a two-level bosonic system. We
develop a framework where coupling to an external device causes the system to
select an eigenstate of the measured observable, and subsequent decoupling
returns the system to the initial thermal state. We find that this requires
the system to be coupled, initially with equal strength, to three harmonic
baths representing elements of the environment. Each bath couples to an
observable of the system represented by one of the Pauli matrices. The system
and baths are illustrated in Figure 1. The resulting random exploration of the
Bloch sphere, with bias brought about by the system Hamiltonian, produces
behaviour consistent with a thermal Gibbs state at high temperature. This
situation is then disturbed by increasing one of the environmental coupling
strengths, as a representation of additional interaction associated with a
measuring device, which obliges the system to move towards and dwell in the
vicinity of an appropriate eigenstate according to the Born rule. In our case,
we will be dealing with continuous measurements of the system energy.
Returning the coupling to its initial strength reverses this dynamical
behaviour, such that we can regard the whole sequence as a simple
representation of the dynamics of quantum measurement. Having established the
stochastic dynamics, we can then derive the stochastic entropy production of
measurement using established analysis [5] and show that measurement is
associated, on average, with positive entropy production.

Figure 1: An open system consisting of a bosonic two-level system, interacting
with an environment represented by three sets of harmonic oscillators, each
coupled via a different Pauli matrix.

Note that incorporating stochasticity into the dynamics is crucial to
computing the irreversibility of a process. It can be demonstrated that the
_mean_ evolution of the reduced density matrix towards stationarity entirely
misses the quantum selection of an eigenstate of the measured observable.
Quantifying the irreversibility of a stochastic trajectory, or indeed the
average irreversibility over an ensemble of trajectories, requires the
incorporation of physical randomness into the dynamics.

The plan for the paper is as follows: in Section 2 the essential ideas of
stochastic thermodynamics are introduced, as well as the SLN equation
describing a system coupled to several independent baths. In Section 3 we
obtain a Markovian approximation of the noise-averaged SLN equation by taking
a high temperature limit and assuming the relaxation times of the environment
to be short. We then construct a stochastic unravelling of the equation of
motion by choosing appropriate Kraus operators and approximations. In Section
4 we study process protocols of variation in the strength of one of the bath
couplings; and then our calculations of the stochastic entropy production are
presented, indicating that there is an entropic cost of measurement, and
showing adherence of the entropy production to the detailed fluctuation
theorem.

##  2 Theory

Continuous or weak quantum measurements are central to the work we will be
discussing. They bring about a continuous change to the system being
monitored, induced by the environment to which the system is coupled.
Continuous quantum measurements produce evolution described by a stochastic
differential equation (SDE), with an entropy production which can then be
assessed by applying methods of stochastic thermodynamics. We describe methods
for stochastic entropy production in Section 2.1, followed by the SLN equation
in Section 2.2, which is the starting point for the dynamics we will be
exploring.

###  2.1 Stochastic Thermodynamics

The concept of entropy arose over a century ago from a consideration of the
irreversibility of macroscopic phenomena, but our understanding has evolved
significantly since then. Insights into entropy production were expanded by
applying thermodynamic concepts to small and individual systems in the form of
fluctuation theorems. These were first introduced by Evans _et al._ [14, 15,
16, 17], and then similar ideas appeared in chaos theory [18] and stochastic
modelling [19, 20], and were also developed by Jarzynski [21] and Crooks [22,
23], amongst others. These developments introduce the idea of entropy
production associated with individual stochastic trajectories.

Stochastic entropy production can be defined as the contrast in likelihoods of
forward and reverse sequences of system behaviour [4]. Specifically, it takes
the form of the logarithm of a ratio of the probabilities of a forward
trajectory driven by a _forward_ protocol of driving forces, and the
corresponding reversed trajectory driven by a _reverse_ protocol. These
protocols are defined such that one is the time-reversed version of the other,
e.g., if the system Hamiltonian varies in a specific way for the forward
protocol, the exact opposite time dependence takes place in the reverse
protocol. It may be shown that irreversible behaviour such as relaxation
towards a stationary state is then accompanied by positive mean stochastic
entropy production.

We can write the total stochastic entropy production associated with the
evolution of a set of time-dependent coordinates $\boldsymbol{x}(t)$ for
$0\leq t\leq\tau$ as

$\Delta
s_{tot}(\boldsymbol{x})=\ln\frac{p^{F}(\boldsymbol{x}(0),0)\mathcal{P}^{F}(\boldsymbol{x}(\tau)|\boldsymbol{x}(0))}{p^{F}(\boldsymbol{x}(\tau),\tau)\mathcal{P}^{R}(\boldsymbol{x}^{\dagger}(\tau)|\boldsymbol{x}^{\dagger}(0))},$
(1)

where $\boldsymbol{x}^{\dagger}$ is a set of coordinates evolving to form a
time-reversed trajectory during the reverse protocol. Concretely, the reverse
trajectory is defined as $\boldsymbol{x}^{\dagger}(t)=\boldsymbol{x}(\tau-t)$,
such that the starting and finishing points $\boldsymbol{x}^{\dagger}(0)$ and
$\boldsymbol{x}^{\dagger}(\tau)$ are $\boldsymbol{x}(\tau)$ and
$\boldsymbol{x}(0)$ respectively. It is a detailed reversal of the sequence of
events described by the forward trajectory. For simplicity, our notation
disregards an inversion in coordinates that are odd under time reversal
symmetry when defining the reverse trajectory [5]. $p^{F}(\boldsymbol{x},t)$
is the probability density function (pdf) over the coordinates obtained from
solving the appropriate Fokker-Planck equation for the forward protocol, while
$\mathcal{P}^{F}$ and $\mathcal{P}^{R}$ are the conditional probability
densities for a trajectory from $\boldsymbol{x}(0)$ to $\boldsymbol{x}(\tau)$
under the forward protocol, and for the time reversed trajectory from
$\boldsymbol{x^{\dagger}}(0)$ to $\boldsymbol{x}^{\dagger}(\tau)$ under the
reverse protocol, respectively.

Denoting the pdfs of the total entropy production $\Delta s_{tot}$ in forward
and reverse driving protocols as $P^{F}\left(\Delta s_{tot}\right)$ and
$P^{R}\left(\Delta s_{tot}\right)$, respectively, it is possible in certain
circumstances to relate them as follows [24, 25]:

$e^{\Delta s_{tot}}=\frac{P^{F}\left(\Delta s_{tot}\right)}{P^{R}\left(-\Delta
s_{tot}\right)},$ (2)

which is called the detailed fluctuation theorem. This result states that the
probability of a negative stochastic entropy production for a reverse protocol
is not zero, in spite of the usual demands of the second law, though it is
exponentially smaller than the probability of a positive production of the
same magnitude during the corresponding forward protocol.

####  2.1.1 SDEs and Fokker-Planck equation

We need to construct an appropriate evolution equation for the reduced density
matrix [7]. Kraus operators define mappings of the density matrix between an
initial and final state, and continuous stochastic trajectories of the density
matrix can be generated using a sequence of Kraus operator maps defined for an
infinitesimal time interval $dt$.

The general Kraus representation of the evolution of a density matrix
$\bar{\rho}$ in a time interval $dt$ is given by

$\bar{\rho}(t+dt)=\sum_{k}M_{k}\bar{\rho}(t)M_{k}^{\dagger},$ (3)

where $M_{k}$ are the Kraus operators, labelled by the index $k$. The map
should be trace preserving, which requires the operator sum identity
$\sum_{k}M_{k}^{\dagger}M_{k}=\mathbb{I}$ to be satisfied ($\mathbb{I}$ is the
identity operator). The Kraus operators depend on $dt$, and for continuous
evolution we require the Kraus operators to differ incrementally from the
identity, i.e. $M_{k}\propto\mathbb{I}$, when $dt\rightarrow 0$.

An interpretation of Eq. (3) is that each Kraus operator can implement a
stochastic action on the system by the environment, namely

$\rho(t+dt)=\rho(t)+d\rho=\frac{M_{k}\rho(t)M_{k}^{\dagger}}{{\rm
Tr}(M_{k}\rho(t)M_{k}^{\dagger})},$ (4)

where $\rho$ is a member of an _ensemble_ of density matrices representing the
uncertain current state of the system. The trace of $\rho$ is clearly
preserved. The operation takes place with a conditional probability given by
$p_{k}={\rm Tr}(M_{k}\rho(t)M_{k}^{\dagger})$ such that the average over all
possible transformations of $\rho(t)$ is given by

$\sum_{k}p_{k}\frac{M_{k}\rho(t)M_{k}^{\dagger}}{{\rm
Tr}(M_{k}\rho(t)M_{k}^{\dagger})}=\sum_{k}M_{k}\rho(t)M_{k}^{\dagger},$ (5)

and a further averaging over the ensemble of $\rho(t)$ yields Eq. (3), with
the over-bar therefore denoting an ensemble average. For this interpretation
to be physically acceptable $\rho$ must remain positive definite under the
mapping Eq. (3). This approach can then lead to a Lindblad equation for the
ensemble averaged density matrix $\bar{\rho}$ [26]:

$\displaystyle d\bar{\rho}$ $\displaystyle=$
$\displaystyle-i\left[H_{sys},\bar{\rho}\right]dt+\sum_{k}\left(L_{k}\bar{\rho}L_{k}^{\dagger}-\frac{1}{2}\left\\{L_{k}^{\dagger}L_{k},\bar{\rho}\right\\}\right)dt,$
(6)

with Lindblad operators $L_{k}$ related to the $M_{k}.$ Such a deterministic
Lindblad equation can be unravelled into a stochastic differential equation to
simulate stochastic interactions such as continuous measurements [27]. We show
how this framework can be implemented for our two-level bosonic system in
Section 3.

To calculate stochastic quantum entropy production, forward and reverse Kraus
operators would be needed to construct forward and reverse stochastic
trajectories, respectively. This was first proposed by Crooks [28] by
considering forward and reverse changes in the density matrix with respect to
the invariant equilibrium state of the system. This approach has been used to
calculate the entropy production corresponding to quantum jump unravellings
[29, 30, 31, 32, 33], although the method would not be appropriate for systems
without an equilibrium state [34, 35]. Other recent developments instead
construct the reverse Kraus operators from the time reversal of the forward
operators [34, 35, 36].

However, we need not seek reverse Kraus operators for situations where the
trajectories are continuous and the stochastic evolution is Markovian. We need
only derive a set of Itô SDEs for the forward dynamics (according to
appropriate forward Kraus operators) from which the entropy production
associated with the evolution of the reduced density matrix may be computed
using the approach developed in [5]. We could construct SDEs for each element
of the density matrix, though it is more convenient, and physically more
transparent, to consider the dynamics of quantum expectation values of various
physical operators, as we shall see.

Together with the set of SDEs, we require an associated Fokker-Planck equation
for the pdf of the chosen system variables. Let us therefore consider an Itô
SDE for a vector of dynamical variables $\boldsymbol{x}$ involving a vector of
Wiener increments $\boldsymbol{dW}$, a vector of deterministic rate functions
$\boldsymbol{A}$ and a matrix of noise coefficients $\boldsymbol{B}$:

$d\boldsymbol{x}=\boldsymbol{A}(\boldsymbol{x},t)dt+\boldsymbol{B}(\boldsymbol{x},t)\boldsymbol{dW}.$
(7)

The corresponding Fokker-Planck equation that describes the evolution of the
pdf $p(\boldsymbol{x},t)$ is [37]

$\frac{\partial p(\boldsymbol{x},t)}{\partial
t}=\sum_{i}\frac{\partial}{\partial
x_{i}}\left[-A_{i}(\boldsymbol{x},t)p(\boldsymbol{x},t)+\sum_{j}\frac{\partial}{\partial
x_{j}}\left(D_{ij}(\boldsymbol{x},t)p(\boldsymbol{x},t)\right)\right],$ (8)

where $\boldsymbol{D=\frac{1}{2}\boldsymbol{B}\boldsymbol{B}^{T}}$ is the
diffusion matrix. We shall find that $\boldsymbol{D}$ is diagonal for our
model which simplifies matters.

####  2.1.2 Stochastic entropy production

From the definition in Eq. (1), it is possible to divide the total entropy
production into two contributions associated with the environment and the
system, respectively. This division is somewhat arbitrary and is introduced
largely for ease of interpreting the behaviour. We calculate the total entropy
production using the approach of Ref. [5], first separating the deterministic
part of Eq. (7) into two contributions:

$d\boldsymbol{x}=\boldsymbol{A}^{irr}(\boldsymbol{x},t)dt+\boldsymbol{A}^{rev}(\boldsymbol{x},t)dt+\boldsymbol{B}(\boldsymbol{x},t)\boldsymbol{dW},$
(9)

where $\boldsymbol{A}^{irr}$ and $\boldsymbol{A}^{rev}$ are the irreversible
and reversible contributions to $\boldsymbol{A}$, respectively. Specifically,
if $dx_{i}$ (the $i$-th component of $d\boldsymbol{x}$) is even with respect
to time reversal, then $A_{i}^{rev}dt$ is also even. Since $dt$ is odd,
$A_{i}^{rev}$ has to be odd: $A_{i}^{rev}$ transforms in the opposite way to
$x_{i}$ under the time reversal. By similar reasoning $A_{i}^{irr}$ transforms
in the same way as $x_{i}$. This separation is necessary for the calculation
of the entropy production [5].

Explicit expressions for calculating the two contributions of the entropy
production have been derived. For a system with a diagonal diffusion matrix
$\boldsymbol{D}$, the incremental _environmental_ stochastic entropy
production is [5]:

$\displaystyle d\Delta s_{env}$
$\displaystyle=\sum_{i=1}^{N}\left[\frac{A_{i}^{ir}}{D_{ii}}dx_{i}-\frac{A_{i}^{rev}A_{i}^{irr}}{D_{ii}}dt+\frac{\partial
A_{i}^{irr}}{\partial x_{i}}dt-\frac{\partial A_{i}^{rev}}{\partial
x_{i}}dt-\frac{1}{D_{ii}}\frac{\partial D_{ii}}{\partial x_{i}}dx_{i}\right.$
(10)
$\displaystyle\left.\qquad\qquad\qquad+\frac{A_{i}^{rev}-A_{i}^{irr}}{D_{ii}}\frac{\partial
D_{ii}}{\partial x_{i}}dt-\frac{\partial^{2}D_{ii}}{\partial
x_{i}^{2}}dt+\frac{1}{D_{ii}}\left(\frac{\partial D_{ii}}{\partial
x_{i}}\right)^{2}dt\right],$

where $N$ is the number of dynamical variables contributing to the entropy
production. Note that this expression (10) contains Wiener increments and
hence the environmental entropy production $\Delta s_{env}$ is explicitly a
stochastic variable.

The second contribution, the incremental _system_ stochastic entropy
production, is related to the change in the logarithm of the pdf
$p(\boldsymbol{x},t)$ of the system, obtained from the Fokker-Planck equation:

$\displaystyle d\Delta s_{sys}$ $\displaystyle=-d\left(\ln
p(\boldsymbol{x},t)\right)$ (11) $\displaystyle=-\frac{\partial\ln
p(\boldsymbol{x},t)}{\partial t}dt-\sum_{i}\frac{\partial\ln
p(\boldsymbol{x},t)}{\partial x_{i}}dx_{i}-\sum_{i}D_{ii}\frac{\partial^{2}\ln
p(\boldsymbol{x},t)}{\partial x_{i}^{2}}dt,$

where the second line is obtained using Itô’s Lemma. The total stochastic
entropy production is then given by the sum of the environmental and system
contributions, $d\Delta s_{tot}=d\Delta s_{sys}+d\Delta s_{env}$. From these
individual stochastic contributions, it is possible to obtain the averaged
system and environmental entropy productions by cumulative summing up of their
respective incremental stochastic contributions, followed by the averaging
over many realisations of the noises and hence system evolution.

The averaged total entropy production also has an analytical form. For certain
circumstances, namely sufficient elapsed time for the pdf to become time
independent, and the maintenance of a condition of detailed balance, this
average is given by the Kullback–Leibler divergence between the initial and
final pdfs [5]:

$\langle\langle\Delta s_{tot}\rangle\rangle=\int
d\boldsymbol{x}\,p(\boldsymbol{x},t_{init})\ln\frac{p(\boldsymbol{x},t_{init})}{p(\boldsymbol{x},t_{final})}.$
(12)

The double angled brackets denote averages over all noise histories and
initial coordinates, over the time interval from $t_{init}$ to $t_{final}$
arising from all initial states in the ensemble. Likewise, the averaged
_system_ entropy production also has an analytical form, usually taken to be

$\langle\langle\Delta s_{sys}\rangle\rangle=\Delta S_{G}=-\int
d\boldsymbol{x}\,p(\boldsymbol{x},t_{final})\,\ln
p(\boldsymbol{x},t_{final})+\int
d\boldsymbol{x}\,p(\boldsymbol{x},t_{init})\,\ln p(\boldsymbol{x},t_{init}),$
(13)

namely the difference in the Gibbs entropy of the final and initial
probability distributions. This result again emerges only under certain
conditions, and we investigate additional contributions in Section 3.3.

###  2.2 Stochastic Liouville-von Neumann equation

Having described some of the background concepts of stochastic entropy
production, we now introduce the stochastic dynamics under consideration. As
briefly discussed in the introduction, a non-Markovian SLN equation can
provide an exact treatment of the dynamics of the open system reduced density
matrix, implicitly after averaging over all manifestations of the noises
associated with the environment (see below). The approach involves taking an
average over stochastic solutions of the equation resulting in a deterministic
trajectory from which physical predictions can be obtained [12, 38].

However, in order to unravel the SLN equation in a straightforward way, we
need to start from a Markovian limit of the system behaviour, such that the
system interacts with an environment that does not possess any memory. Another
perspective is to consider the environmental correlation times, which in a
Markovian limit are very short when compared to the characteristic timescales
for the dynamics of the open system. In this section we obtain the
deterministic Markovian equation of motion for the physically meaningful,
ensemble averaged reduced density matrix of our system. In Section 3 that
follows, the unravelling procedure is considered.

Many non-Markovian methods for open quantum systems begin their development
with the full density matrix of the environment and system and proceed to take
a partial trace in order to derive an equation of motion exclusively for the
reduced density matrix of the system. The Feynman-Vernon influence functional
formalism is one such method, where the reduced density matrix of the open
system is represented as a path integral over all environmental modes made up
of an infinite number of harmonic oscillators [13]. Methods which rely on this
include quasi-adiabatic path integrals [39], hierarchical equations of motion
[40, 41], stochastic Schrödinger equations [42] and the Stochastic Liouville-
von Neumann (SLN) equation and its variations [8, 9, 10, 11, 12]. Other
methods that do not rely on the Feynman-Vernon functional include projection
operator methods such as the Nakajima-Zwanzig equation [43, 44, 45].

The SLN equation is a stochastic differential equation containing complex
cross-correlated coloured noises representing the environment; to arrive at
the appropriate deterministic description, an average needs to be taken over
such noises. Note again that the outcome of such a process will describe the
evolution of an average of an ensemble of physical density matrices, according
to the framework of interpretation set out in Section 2.1.1.

The SLN equation describes an open quantum system, with coordinates $q$,
coupled to a large number of harmonic oscillators that represent the
environment. Here we shall consider the system to be coupled to three
independent sets of oscillators; the reason for doing so will become apparent.
The SLN formalism [9] can easily be extended to accommodate this. The
Hamiltonian coupling of the system to the $k^{th}$ set of oscillators is taken
to be linear in the oscillator coordinates $\xi_{ik}$, where $i$ labels the
oscillators, but it can depend on a general coupling operator $f_{k}(q)$ of
the system coordinates $q$. The full Hamiltonian can be written as follows

$H_{total}(q,\\{\xi_{ik}\\},t)=H_{sys}(q,t)+\sum_{k}\left[H_{env;k}(\xi_{ik})-\sum_{i}\xi_{ik}f_{k}(q,t)\right],$
(14)

where $H_{env;k}(\xi_{ik})$ is the Hamiltonian for the $k^{th}$ set of
environmental oscillators, and $H_{sys}$ is the Hamiltonian of the open
system. To derive the SLN equation, the initial density matrix of the full
system is taken to be the tensor product of the reduced density matrix of the
system $\rho_{q}(t_{0})$ and the environmental density matrix
$\rho_{\xi}(t_{0})$, i.e. $\rho_{0}=\rho_{q}(t_{0})\otimes\rho_{\xi}(t_{0})$,
where $t_{0}$ is an initial time. We shall be considering a weak coupling
limit so this partitioned state is an appropriate approximation for the system
we are investigating.

Based on the Hamiltonian of Eq. (14), it is possible to obtain a stochastic
differential equation for the (unphysical) stochastic reduced density matrix
$\rho_{s}(t)$ of the system driven by a pair of complex coloured Gaussian
noises, $\eta_{k}(t)$, and $\nu_{k}(t)$, for each oscillator set $k$ [8, 9,
10, 11]. It should be stressed again that these stochastic noises are purely
mathematical constructs, and do not generate individual physical trajectories
[11, 12]. Physical results only arise when the noise-driven trajectories are
averaged over their many realisations. The SLN equation (with $\hbar=1$)
becomes:

$\frac{d\rho_{s}(t)}{dt}=-i\left[H_{sys},\rho_{s}(t)\right]+i\sum_{k}\left(\eta_{k}(t)\left[f_{k},\rho_{s}(t)\right]+\frac{\nu_{k}(t)}{2}\left\\{f_{k},\rho_{s}(t)\right\\}\right),$
(15)

where the square brackets correspond to a commutator, and the curly brackets
to an anticommutator. Each set of environmental oscillators is connected
independently to the open system and the associated noises satisfy the
following correlations:

$\langle\eta_{k}(t)\eta_{k}(t^{\prime})\rangle=\int_{0}^{\infty}\frac{d\omega}{\pi}J_{k}(\omega)\coth\left(\frac{1}{2}\beta_{k}\omega\right)\cos\left(\omega\left(t-t^{\prime}\right)\right)\equiv
K_{k}^{\eta\eta}(t-t^{\prime}),$ (16)
$\langle\eta_{k}(t)\nu_{k}(t^{\prime})\rangle=-2i\Theta(t-t^{\prime})\int_{0}^{\infty}\frac{d\omega}{\pi}J_{k}(\omega)\sin\left(\omega\left(t-t^{\prime}\right)\right)\equiv
K_{k}^{\eta\nu}(t-t^{\prime}),$ (17)

where $J_{k}(\omega)$ is the spectral density of the $k^{th}$ oscillator set,
$\beta_{k}=1/k_{B}T_{k}$ is its inverse temperature (with $k_{B}$ being set to
1 hereafter), $\Theta(t-t^{\prime})$ the Heaviside function, and the angled
brackets denote an average over the environmental noises; all other
correlation functions are zero. There exist several ways to construct the
coloured noises, with each scheme affecting the convergence of results in
different ways [38].

###  2.3 Noise-averaged SLN equation

We need to average Eq. (15) with respect to realisations of the environmental
noises (to be denoted by single angle brackets) to obtain an exact and
deterministic evolution equation for the physical reduced density matrix of
the open system. Taking the average of both sides of Eq. (15) and denoting the
physical density matrix averaged over the noises as
$\bar{\rho}=\langle\rho_{s}\rangle$, we write:

$\frac{d\bar{\rho}(t)}{dt}=-i\left[H_{sys},\bar{\rho}(t)\right]+i\sum_{k}\left(\left[f_{k},\langle\eta_{k}(t)\rho_{s}(t)\rangle\right]+\frac{1}{2}\left\\{f_{k},\langle\nu_{k}(t)\rho_{s}(t)\rangle\right\\}\right).$
(18)

In order to calculate the averages of the stochastic density matrix multiplied
by the noises that appear in the right hand side, the Furutsu-Novikov theorem
may be used [46, 47]. For a set of noises $\zeta_{i}$ with the correlation

$\langle\zeta_{i}(t)\zeta_{j}(t^{\prime})\rangle=F_{ij}(t,t^{\prime}),$ (19)

the Furutsu-Novikov theorem states that

$\langle\zeta_{i}(t)A[\zeta]\rangle=\sum_{j}\int_{0}^{t}dt^{\prime}F_{ij}(t,t^{\prime})\left\langle\frac{\delta
A[\phi]}{\delta\zeta_{j}(t^{\prime})}\right\rangle,$ (20)

where $A[\zeta]$ is a functional of the noises and $\delta
A[\zeta]/\delta\zeta_{j}$ is its functional derivative with respect to
$\zeta_{j}$. The averaged product of a noise and a noise-dependent functional
may therefore be transformed into an integral. We use Eqs. (20) in (18) to
obtain the required explicit expressions for the averages of the noises
multiplied by $\rho_{s}$:

$\displaystyle\frac{d\bar{\rho}(t)}{dt}$
$\displaystyle=-i\left[H(t),\bar{\rho}(t)\right]+i\sum_{k}\left[f_{k},\int_{0}^{t}dt^{\prime}K_{k}^{\eta\eta}(t-t^{\prime})\left\langle\frac{\delta\rho_{s}(t)}{\delta\eta_{k}(t^{\prime})}\right\rangle+\int_{0}^{t}dt^{\prime}K_{k}^{\eta\nu}(t-t^{\prime})\left\langle\frac{\delta\rho_{s}(t)}{\delta\nu_{k}(t^{\prime})}\right\rangle\right]$
(21)
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad+\frac{i}{2}\sum_{k}\left\\{f_{k},\int_{0}^{t}dt^{\prime}K_{k}^{\eta\nu}(t^{\prime}-t)\left\langle\frac{\delta\rho_{s}(t)}{\delta\eta_{k}(t^{\prime})}\right\rangle\right\\}.$

From Eq. (17) we observe the presence of the Heaviside function
$\Theta(t-t^{\prime})$, which vanishes for $t^{\prime}>t.$ The integrals in
Eq. (21) have an upper limit $t$, which implies that
$K_{k}^{\eta\nu}(t^{\prime}-t)=0$ in the last term since $t^{\prime}$ is
always smaller than $t$ due to the integration limits. Hence, the last term is
zero and can be dropped; only the commutator terms remain. We then need to
calculate the functional derivatives, and after some manipulation, see
Appendix A.1, we obtain

$\displaystyle\frac{d\bar{\rho}(t)}{dt}$
$\displaystyle=-i\left[H_{sys},\bar{\rho}(t)\right]-\int_{0}^{t}dt^{\prime}K^{\eta\eta}(t-t^{\prime})\sum_{k}\left[f_{k},\left\langle
U_{+}(t,t^{\prime})f_{k}(t^{\prime})U_{+}(t^{\prime},t)\rho_{s}(t)-\rho_{s}(t)U_{-}(t,t^{\prime})f_{k}(t^{\prime})U_{-}(t^{\prime},t)\right\rangle\right]$
(22)
$\displaystyle\qquad\qquad\qquad+\frac{1}{2}\int_{0}^{t}dt^{\prime}K_{k}^{\eta\nu}(t-t^{\prime})\sum_{k}\left[f_{k},\left\langle
U_{+}(t,t^{\prime})f_{k}(t^{\prime})U_{+}(t^{\prime},t)\rho_{s}(t)+\rho_{s}(t)U_{-}(t,t^{\prime})f_{k}(t^{\prime})U_{-}(t^{\prime},t)\right\rangle\right],$

where $U_{+}$ and $U_{-}$ are forward and reverse propagators, respectively
(see Appendix A.1). Eq. (22) is the most general form of the averaged SLN
equation, as no approximations have been made so far. Note its non-Markovian
nature. Also note that it describes the evolution of a reduced density matrix
averaged over physical fluctuations, along the lines of Eq. (6). However, this
equation still contains the stochastic density matrix $\rho_{s}$ and hence is
not a self-contained equation for the physical ensemble-averaged density
matrix $\bar{\rho}$. To be able to work with an equation depending exclusively
on the latter, some simplifications need to be made; this will be done in the
next Section.

###  2.4 Markovian limit of the SLN equation

Eq. (22) in its exact non-Markovian form is difficult to work with as it still
contains the stochastic density matrix. A self-contained equation for the
physical density matrix $\bar{\rho}(t)$ can be obtained, however, in the
Markovian limit by approximating certain properties of the environment. We
begin by assuming that all the sets of oscillators coupled to the system have
the same temperature $T_{k}\equiv T$, and the same spectral density. This is
not a necessary assumption or approximation in our analysis, but it will
simplify our notations. Next we take the limit in which the temperature of the
environment $T$ is very large such that $T$ is much greater than any
characteristic energy quantum of its harmonic oscillators, $\beta\omega\ll 1$,
and we will also adopt the same Ohmic spectral density
$J(\omega)=\alpha\omega$ for all oscillator sets, where $\alpha$ is a
proportionality constant. In this case we can use the first order expansion of
$\coth\left(\frac{\beta\omega}{2}\right)\approx\frac{2}{\beta\omega}$,
following [48], and obtain the following results for the correlation functions
of Eqs. (16) and (17):

$K^{\eta\eta}(t-t^{\prime})\approx\int_{0}^{\infty}\frac{d\omega}{\pi}\alpha\omega\frac{2}{\beta\omega}\cos\left(\omega(t-t^{\prime})\right)=\frac{2\alpha}{\beta}\delta(t-t^{\prime}),$
(23)

and

$\displaystyle K^{\eta\nu}(t-t^{\prime})$
$\displaystyle=2i\alpha\Theta(t-t^{\prime})\frac{\partial}{\partial
t}\left[\frac{1}{\pi}\int_{0}^{\infty}d\omega\cos\left(\omega(t-t^{\prime})\right)\right]=-2i\alpha\Theta(t-t^{\prime})\frac{\partial}{\partial
t^{\prime}}\left[\frac{1}{\pi}\int_{0}^{\infty}d\omega\cos\left(\omega(t-t^{\prime})\right)\right]$
(24) $\displaystyle=-2i\alpha\Theta(t-t^{\prime})\frac{\partial}{\partial
t^{\prime}}\delta(t-t^{\prime}).$

Note that the correlation functions are the same for all the oscillator sets,
hence we have dropped the index $k$ here. Eqs. (23) and (24) appear in [49],
although a different route to obtain them was taken, specifically that of
introducing a cutoff in the spectral density and allowing it to be much larger
than the dynamical timescales of the system. Inserting Eqs. (23) and (24) into
Eq. (22) leads to the Markovian limit of the averaged SLN equation that
contains only $\bar{\rho}(t)$. The equation takes a similar form to that found
in other work [48]:

$\displaystyle\frac{d\bar{\rho}(t)}{dt}$
$\displaystyle=-i\left[H_{sys},\bar{\rho}(t)\right]-\sum_{k}\left[f_{k},\frac{\alpha}{\beta}\left[f_{k},\bar{\rho}(t)\right]-\frac{\alpha}{2}\left\\{\left[H_{sys},f_{k}\right],\bar{\rho}(t)\right\\}+\frac{i\alpha}{2}\left\\{\frac{\partial
f_{k}}{\partial t},\bar{\rho}(t)\right\\}\right]$
$\displaystyle=-i\left[H_{sys},\bar{\rho}(t)\right]+\sum_{k}\left(-\frac{\alpha}{\beta}\left[f_{k},\left[f_{k},\bar{\rho}(t)\right]\right]+\frac{\alpha}{2}\left[f_{k},\left\\{\left[H_{sys},f_{k}\right],\bar{\rho}(t)\right\\}\right]-\frac{i\alpha}{2}\left[f_{k},\left\\{\frac{\partial
f_{k}}{\partial t},\bar{\rho}(t)\right\\}\right]\right),$ (25)

but unlike the result published in [48], this equation contains a term
involving the time dependence of the system coupling operators $f_{k}$.

##  3 System setup and equations of motion

In the last section we obtained the Markovian limit of the evolution equation
for the ensemble averaged reduced density matrix. Here we shall unravel this
equation, namely add noise terms, such that individual trajectories correspond
to different physical realisations of the quantum state diffusion of the
system under the influence of an underspecified environment.

###  3.1 Setup of the system and the unravelling procedure

####  3.1.1 Two-level system and choice of environments

We consider a two-level bosonic system with Hamiltonian

$H_{sys}=\epsilon\sigma_{z},$ (26)

where $\sigma_{z}$ is the usual Pauli matrix and $\epsilon$ has the dimension
of energy. We take the system Hamiltonian to be constant in time. We require
the environmental interaction with the open system to be consistent with a
thermal state subject to our earlier assumptions of weak coupling and high
temperature. This means that the insertion of
$\bar{\rho}\propto\mathbb{I}-\beta H_{sys}$ should cause the right hand side
of Eq. (25), with time independent $f_{k}$, to vanish:

$\displaystyle\frac{d\bar{\rho}}{dt}$
$\displaystyle\propto-i\left[H_{sys},\mathbb{I}-\beta
H_{sys}\right]+\frac{\alpha}{\beta}\sum_{k}\left(-\left[f_{k},\left[f_{k},\mathbb{I}-\beta
H_{sys}\right]\right]+\frac{1}{2}\left[f_{k},\left\\{\left[\beta
H_{sys},f_{k}\right],\mathbb{I}-\beta H_{sys}\right\\}\right]\right)$
$\displaystyle=-\frac{\alpha\beta}{2}\sum_{k}\left[f_{k},\left\\{\left[H_{sys},f_{k}\right],H_{sys}\right\\}\right].$
(27)

The anticommutator vanishes exactly for our choice of Hamiltonian where
$H_{sys}^{2}\propto\mathbb{I}$ and the implication is that any choice of the
operators $f_{k}$, through which the system interacts with the environment, is
consistent with the appropriate stationary thermal state. Since the $f_{k}$
are hermitian according to the structure of the Hamiltonian in Eq. (14), it is
natural, in our case, to employ three coupling operators proportional to the
Pauli matrices. Furthermore, it is natural for the initial strength of the
coupling to be the same for each operator, in order that the system-
environment interaction should be isotropic. Note that the same framework, if
needed, could then be employed for system Hamiltonians proportional to
$\sigma_{x}$ or $\sigma_{y}$ as well as the choice in Eq. (26).

This coupling scheme has the added advantage that it allows us to model the
act of measurement of the system energy in a straightforward way, merely by
elevating the strength of the coupling to the operator $\sigma_{z}$. This
represents an additional interaction between the system and an external
measuring device, and in the next section we shall see that the unravelled
system dynamics under such an interaction can indeed correspond to the
stochastic selection of one of the energy eigenstates. Specifically, the
Hamiltonian, Eq. (26), has two eigenvectors $|\pm\epsilon\rangle$ with
eigenvalues $\pm\epsilon$, respectively, and the environmental coupling to
$\sigma_{z}$ drives the density matrix stochastically towards the form
$|+\epsilon\rangle\langle+\epsilon|$ or $|-\epsilon\rangle\langle-\epsilon|$.
However, the scheme of coupling to the three Pauli matrices means that these
eigenstates are not fixed points of the dynamics: the tendency is only for the
system to dwell in the vicinity of one or other of the eigenstates, to an
extent that depends on the strength of environmental coupling to $\sigma_{z}$.

To summarise, we shall couple three independent sets of harmonic oscillators
to our open system using operators:

$f_{x}=\gamma_{0}\sigma_{x}\,,\,\,\,f_{y}=\gamma_{0}\sigma_{y}\,\,\,\text{and}\,\,\,f_{z}=\gamma(t)\sigma_{z},$
(28)

where for simplicity we replace labels $k$ by the Cartesian indices $x$, $y$
and $z$. $\gamma_{0}$ is the coupling coefficient between the system and the
environment, and $\delta\gamma(t)=\gamma(t)-\gamma_{0}$ is the coupling
coefficient between the system and the device measuring $H_{sys}$. If
$\gamma(t)$ is increased from an initial value of $\gamma_{0}$, we can model
measurement of the system energy while under the influence of a thermal bath.
We can calculate the stochastic entropy production related to this process.
Furthermore, we can model the detachment of the measuring device by a protocol
where $\gamma(t)$ returns to $\gamma_{0}$ allowing the system to re-
thermalise, and again compute the associated entropy production.

Note also that, with this choice of the coupling, the last term in Eq. (25)
containing the time derivative of the $f_{k}$ vanishes exactly; this is
obvious for constant $f_{x}$ and $f_{y}$, while for $f_{z}$ it follows from
the fact that $\partial f_{k}/\partial t$ is proportional to $\sigma_{z}$.

####  3.1.2 Unravelled stochastic equation

To obtain a form of Eq. (25) that allows for a straightforward stochastic
unravelling, we would need to write it in a Lindblad form with positive
coefficients [7]. It would be possible to do so by extending the Hilbert space
and then adding stochasticity [50, 51], though that approach will not be taken
here. Instead we construct Lindblad operators that match the dynamics of Eq.
(25) in accordance with our high temperature approximation $\beta H_{sys}\ll
1$.

We shall now show that the Markovian averaged (deterministic) SLN Eq. (25) is
equivalent, to lowest order in $\beta\epsilon$, to the Lindblad equation when
an appropriate choice of the Lindblad operators is made. Indeed, consider the
Lindblad equation (6) using operators

$L_{k}=\lambda\left(f_{k}-\frac{1}{4}\left[\beta H_{sys},f_{k}\right]\right),$
(29)

where $\lambda=\sqrt{\frac{2\alpha}{\beta}}$. By inserting these operators
into Eq. (6) we obtain

$\displaystyle d\bar{\rho}$
$\displaystyle=-i\left[H_{sys}-\frac{i\alpha}{4}\sum_{k}\left[H_{sys},f_{k}^{2}\right],\bar{\rho}\right]dt+\sum_{k}\left(-\frac{\lambda^{2}}{2}\left[f_{k},\left[f_{k},\bar{\rho}\right]\right]+\frac{\lambda^{2}}{4}\beta\left[f_{k},\left\\{\left[H_{sys},f_{k}\right],\bar{\rho}\right\\}\right]\right)dt$
(30)
$\displaystyle\qquad\qquad\qquad\qquad\qquad-\frac{\lambda^{2}}{16}\sum_{k}\beta^{2}\left(\left[H_{sys},f_{k}\right]\bar{\rho}\left[H_{sys},f_{k}\right]-\frac{1}{2}\left\\{\left[H_{sys},f_{k}\right]^{2},\bar{\rho}\right\\}\right)dt.$

We note that $\left[H_{sys},f_{k}^{2}\right]$ vanishes when
$f_{k}^{2}\propto\mathbb{I}$, which holds in our case. The final contribution
on the right hand side is a factor $\beta H_{sys}$ smaller in magnitude than
the penultimate term, and therefore can be neglected. With these
considerations we recover the Markovian averaged SLN Eq. (25) confirming the
suitability of the Lindblad operators in Eq. (29).

Starting from Eqs. (6) and (29) we can now construct an unravelling consistent
with the quantum state diffusion (QSD) approach [52, 53, 54, 55, 56]. Using
the set of Kraus operators [26, 57, 58]

$M_{\pm
k}\propto\left(\mathbb{I}-iH_{sys}dt-\frac{1}{2}L_{k}^{\dagger}L_{k}dt\pm
L_{k}\sqrt{dt}\right),$ (31)

consistent with the development in Section 2.1.1, and substituting Eq. (31)
into Eq. (4) we obtain (see Refs. [26, 58]):

$d\rho=-i\left[H_{sys},\rho\right]dt+\sum_{k}\left[\left(L_{k}\rho
L_{k}^{\dagger}-\frac{1}{2}\left\\{L_{k}^{\dagger}L_{k},\rho\right\\}\right)dt+\left(\rho
L_{k}^{\dagger}+L_{k}\rho-{\rm
Tr}\left[\rho\left(L_{k}+L_{k}^{\dagger}\right)\right]\rho\right)dW_{k}\right],$
(32)

where we use $\rho$ to denote the physical reduced density matrix of the
system that evolves stochastically, and where environmental noise is
represented by a set of independent Wiener increments $dW_{k}$. Notice that
upon averaging over the Wiener noises this evolution equation corresponds to
Eq. (6) and hence the density matrix $\bar{\rho}$ corresponds to the noise
averaged density matrix $\rho$. Eq. (32) finally provides us with the dynamics
and the means by which to determine the irreversibility of system behaviour in
various cases.

####  3.1.3 Equations of motion for components of the coherence vector

With the particular environmental couplings defined in Eq. (28) ($k=x,y,z$),
and with $\gamma_{0}=1$, we obtain Lindblad operators:

$L_{x}=\lambda\left(\sigma_{x}-i\frac{\beta\epsilon}{2}\sigma_{y}\right),\hskip
18.00003ptL_{y}=\lambda\left(\sigma_{y}+i\frac{\beta\epsilon}{2}\sigma_{x}\right),\hskip
18.00003ptL_{z}=\lambda\gamma\sigma_{z}.$ (33)

Instead of simulating the dynamics of elements of the density matrix directly,
we will evolve the coherence vector defined by the three components
$r_{i}=\text{Tr}(\sigma_{i}\rho)$. This is consistent with a representation of
the density matrix in the form

$\rho=\frac{1}{2}\left(\mathbb{I}+\boldsymbol{r}\cdot\boldsymbol{\sigma}\right).$
(34)

We obtain from Eq. (32) the following equations:

$dr_{x}=-2\epsilon
r_{y}dt-2\lambda^{2}\left(1+\gamma^{2}\right)r_{x}dt+\lambda\left(\beta\epsilon
r_{z}+2\left(1-r_{x}^{2}\right)\right)dW_{x}-2\lambda
r_{x}r_{y}dW_{y}-2\gamma\lambda r_{x}r_{z}dW_{z}$ (35)

$dr_{y}=2\epsilon
r_{x}dt-2\lambda^{2}\left(1+\gamma^{2}\right)r_{y}dt-2\lambda
r_{x}r_{y}dW_{x}+\lambda\left(\beta\epsilon
r_{z}+2\left(1-r_{y}^{2}\right)\right)dW_{y}-2\gamma\lambda r_{y}r_{z}dW_{z}$
(36)

$dr_{z}=-4\lambda^{2}\left(\beta\epsilon+r_{z}\right)dt-\lambda
r_{x}\left(\beta\epsilon+2r_{z}\right)dW_{x}-\lambda
r_{y}\left(\beta\epsilon+2r_{z}\right)dW_{y}+2\gamma\lambda\left(1-r_{z}^{2}\right)dW_{z}.$
(37)

Notice how important it is to consider the _stochastic_ evolution of the
coherence vector rather than the averaged behaviour. For the latter situation,
described by $d\langle r_{x}\rangle=-2\epsilon\langle r_{y}\rangle
dt-2\lambda^{2}\left(1+\gamma^{2}\right)\langle r_{x}\rangle dt$, etc, any
initial state results in $\langle r_{x}\rangle\rightarrow 0$, $\langle
r_{y}\rangle\rightarrow 0$ and $\langle r_{z}\rangle\rightarrow-\beta\epsilon$
at long times and hence is not disturbed by the evolution of $\gamma$ away
from unity. What does change in these circumstances is the pdf of the system
variables, and hence higher moments of the components of the coherence vector.
The noise-averaged Lindblad equation for $\bar{\rho}$ cannot capture the
stochastic selection of an eigenstate under measurement.

####  3.1.4 Equations of motion in cylindrical coordinates

Next we derive the equations of motion for the convenient set of coordinates
$(r^{2},r_{z},\phi)$, where

$r^{2}=r_{x}^{2}+r_{y}^{2}+r_{z}^{2}\,\,\,\text{and}\,\,\,\phi=\arctan\left(\frac{r_{y}}{r_{x}}\right).$
(38)

Using Itô’s Lemma, see Appendix A.2, we can obtain:

$\displaystyle dr^{2}$
$\displaystyle=4\lambda^{2}\left(r^{2}-1\right)\left(\gamma^{2}r_{z}^{2}-\gamma^{2}+r^{2}-r_{z}^{2}-2\right)dt+4\lambda\sqrt{r^{2}-r_{z}^{2}}\cos\phi\left(1-r^{2}\right)dW_{x}$
(39)
$\displaystyle+4\lambda\sqrt{r^{2}-r_{z}^{2}}\sin\phi\left(1-r^{2}\right)dW_{y}+4\lambda\gamma
r_{z}\left(1-r^{2}\right)dW_{z}\,,$

$\displaystyle dr_{z}$
$\displaystyle=-4\lambda^{2}\left(\beta\epsilon+r_{z}\right)dt-\lambda\left(\beta\epsilon+2r_{z}\right)\sqrt{r^{2}-r_{z}^{2}}\cos\phi\,dW_{x}$
$\displaystyle-\lambda\left(\beta\epsilon+2r_{z}\right)\sqrt{r^{2}-r_{z}^{2}}\sin\phi\,dW_{y}+2\gamma\lambda\left(1-r_{z}^{2}\right)\,dW_{z}\,,$
(40)

and

$d\phi=2\epsilon dt-\lambda\frac{\left(\beta\epsilon
r_{z}+2\right)\sin\phi}{\sqrt{r^{2}-r_{z}^{2}}}dW_{x}+\lambda\frac{\left(\beta\epsilon
r_{z}+2\right)\cos\phi}{\sqrt{r^{2}-r_{z}^{2}}}dW_{y}\,.$ (41)

Due to the coordinate singularities, suitable care needs to be taken with
regard to initial state choices and the dynamics themselves.

####  3.1.5 Initial state simplification

We see from Eq. (39) that $r=1$ is a fixed point of the dynamics. Defining a
pure state by the condition ${\rm Tr}(\rho^{2})=1$, it can be shown that $r=1$
corresponds to a pure state by inserting $\rho$ in its coherence vector
representation 34. In order to reduce the computational difficulty of solving
the equations, we set $r=1$ and consider

$\displaystyle dr_{z}$
$\displaystyle=-4\lambda^{2}\left(\beta\epsilon+r_{z}\right)dt-\lambda\left(\beta\epsilon+2r_{z}\right)\sqrt{1-r_{z}^{2}}\cos\phi\,dW_{x}$
(42)
$\displaystyle-\lambda\left(\beta\epsilon+2r_{z}\right)\sqrt{1-r_{z}^{2}}\sin\phi\,dW_{y}+2\gamma\lambda\left(1-r_{z}^{2}\right)\,dW_{z}$

$d\phi=2\epsilon dt-\lambda\frac{\left(\beta\epsilon
r_{z}+2\right)\sin\phi}{\sqrt{1-r_{z}^{2}}}dW_{x}+\lambda\frac{\left(\beta\epsilon
r_{z}+2\right)\cos\phi}{\sqrt{1-r_{z}^{2}}}dW_{y}.$ (43)

Hence, the coherence vector describing our system will always lie on the
surface of the Bloch sphere, and the pdf will depend only on two variables:
$r_{z}$ and $\phi$. These equations form the basis of our simulations with
corresponding results to be presented in Section 4.

###  3.2 Fokker-Planck equation

The definition of the entropy production in Eq. (1) requires an evolving pdf
for the system coordinates over time, and hence we need to derive and solve
the associated Fokker-Planck equation. To do so, we need expressions for the
vector $\boldsymbol{A}$ and matrices $\boldsymbol{B}$ and
$\boldsymbol{D}=\frac{1}{2}\boldsymbol{B}\boldsymbol{B}^{T}$, obtained by
comparing their definitions in Eq. (7) with the equations of motion (42) and
(43). We get

$\boldsymbol{A}=\left(\begin{matrix}A_{z}\\\
A_{\phi}\end{matrix}\right)=\left(\begin{matrix}-4\lambda^{2}(\beta\epsilon+r_{z})\\\
2\epsilon\end{matrix}\right),$ (44)

$\boldsymbol{B}=\lambda\left(\begin{matrix}-(\beta\epsilon+2r_{z})\sqrt{1-r_{z}^{2}}\cos\phi&-(\beta\epsilon+2r_{z})\sqrt{1-r_{z}^{2}}\sin\phi&2\gamma(1-r_{z}^{2})\\\
-(\beta\epsilon r_{z}+2)\frac{\sin\phi}{\sqrt{1-r_{z}^{2}}}&(\beta\epsilon
r_{z}+2)\frac{\cos\phi}{\sqrt{1-r_{z}^{2}}}&0\end{matrix}\right),$ (45)

and hence

$\boldsymbol{D}=\left(\begin{array}[]{cc}D_{zz}&0\\\
0&D_{\phi\phi}\end{array}\right)$ (46)

with

$D_{zz}=2\lambda^{2}\left(1-r_{z}^{2}\right)\left(\left(\frac{\beta\epsilon}{2}\right)^{2}+\beta\epsilon
r_{z}+(1-\gamma^{2})r_{z}^{2}+\gamma^{2}\right)$ (47)

and

$D_{\phi\phi}=2\lambda^{2}\frac{\left(\frac{\beta\epsilon
r_{z}}{2}\right)^{2}+\beta\epsilon r_{z}+1}{1-r_{z}^{2}}\,.$ (48)

There exists no $\phi$ dependence in either $\boldsymbol{A}$ or
$\boldsymbol{D}$, and the latter matrix is indeed diagonal, as anticipated.
With these expressions, the Fokker-Planck equation for the pdf
$p(r_{z},\phi,t)$ becomes:

$\displaystyle\frac{\partial}{\partial t}p(r_{z},\phi,t)$
$\displaystyle=\frac{\partial}{\partial
r_{z}}\left[-A_{z}p(r_{z},\phi,t)+\frac{\partial}{\partial
r_{z}}\left(D_{zz}p(r_{z},\phi,t)\right)\right]+\frac{\partial}{\partial\phi}\left[-A_{\phi}p(r_{z},\phi,t)+\frac{\partial}{\partial\phi}\left(D_{\phi\phi}p(r_{z},\phi,t)\right)\right]$
(49) $\displaystyle=\frac{\partial}{\partial
r_{z}}\left[\frac{\lambda^{2}}{2}\left(-2\beta^{2}\epsilon^{2}r_{z}+12\beta\epsilon\left(1-r_{z}^{2}\right)+16(1-\gamma^{2})r_{z}\left(1-r_{z}^{2}\right)\right)p(r_{z},\phi,t)\right.$
$\displaystyle\qquad\qquad\left.+2\lambda^{2}\left(1-r_{z}^{2}\right)\left(\left(\frac{\beta\epsilon}{2}\right)^{2}+\beta\epsilon
r_{z}+(1-\gamma^{2})r_{z}^{2}+\gamma^{2}\right)\frac{\partial}{\partial
r_{z}}p(r_{z},\phi,t)\right]$
$\displaystyle\qquad\qquad+\frac{\partial}{\partial\phi}\left[-2\epsilon
p(r_{z},\phi,t)+2\lambda^{2}\frac{\left(\frac{\beta\epsilon
r_{z}}{2}\right)^{2}+\beta\epsilon
r_{z}+1}{1-r_{z}^{2}}\frac{\partial}{\partial\phi}p(r_{z},\phi,t)\right].$

The Fokker-Planck equation can be simplified further if we assume the initial
pdf to be independent of $\phi$. As there is no explicit $\phi$ dependence in
Eq. (49), if the initial pdf is $\phi$ independent, it will always remain so.
This leads us to the simplification $p(r_{z},\phi,t)\rightarrow p(r_{z},t)$,
and a more succinct equation:

$\displaystyle\frac{\partial}{\partial t}p(r_{z},t)$
$\displaystyle=\frac{\partial}{\partial
r_{z}}\left[\frac{\lambda^{2}}{2}\left(-2\beta^{2}\epsilon^{2}r_{z}+12\beta\epsilon\left(1-r_{z}^{2}\right)+16(1-\gamma^{2})r_{z}\left(1-r_{z}^{2}\right)\right)p(r_{z},t)\right.$
(50)
$\displaystyle\left.\qquad\qquad+2\lambda^{2}\left(1-r_{z}^{2}\right)\left(\left(\frac{\beta\epsilon}{2}\right)^{2}+\beta\epsilon
r_{z}+(1-\gamma^{2})r_{z}^{2}+\gamma^{2}\right)\frac{\partial}{\partial
r_{z}}p(r_{z},t)\right].$

The stationary solution $p_{st}(r_{z})$ of the Fokker-Planck equation
corresponding to setting $\partial p(r_{z},t)/\partial t=0$ can be obtained
analytically using tools such as Mathematica. For $\gamma\neq 1$, we have:

$\displaystyle p_{st}(r_{z})_{\gamma\neq 1}$
$\displaystyle=\frac{1}{N_{\gamma}}(1-r_{z})^{a}(1+r_{z})^{b}\left(\beta^{2}\epsilon^{2}+4\beta\epsilon
r_{z}+4\gamma^{2}+4(1-\gamma^{2})\right)^{c}$ (51)
$\displaystyle\qquad\qquad\qquad\times\exp\left[\frac{8\beta\epsilon\left(\beta^{2}\epsilon^{2}(1+2\gamma^{2})-4\right){\rm
arctanh}\left[\frac{-\beta\epsilon-2r_{z}(1-\gamma^{2})}{\gamma\sqrt{4\gamma^{2}-4+\beta^{2}\epsilon^{2}}}\right]}{(4-\beta^{2}\epsilon^{2})^{2}\gamma\sqrt{4\gamma^{2}-4+\beta^{2}\epsilon^{2}}}\right],$

while for $\gamma=1$ we have

$p_{st}(r_{z})_{\gamma=1}=\frac{1}{N_{1}}(1-r_{z})^{a}(1+r_{z})^{b}(\beta^{2}\epsilon^{2}+4\beta\epsilon
r_{z}+4)^{d},$ (52)

where $N_{\gamma}$ is a $\gamma$ dependent normalisation factor, and the
exponents are given by

$a=-\left(\frac{\beta\epsilon}{2+\beta\epsilon}\right)^{2},\,\,\,b=-\left(\frac{\beta\epsilon}{2-\beta\epsilon}\right)^{2},\,\,\,c=-1+\frac{4(-4+3\beta^{2}\epsilon^{2})}{(-4+\beta^{2}\epsilon^{2})^{2}},\,\,\,d=-1+\frac{8(-4+3\beta^{2}\epsilon^{2})}{(-4+\beta^{2}\epsilon^{2})^{2}},$
(53)

where we note that $c\neq d$. The stationary pdfs in Eqs. (51) and (52) are
shown in Figure 2 for a number of values of $\gamma$.

Figure 2: Stationary pdf solutions of the Fokker-Planck Eq. (50) for different
values of $\gamma$. Parameters used are $\beta=0.1$ and $\epsilon=1$.

We can see that the coupling strength $\gamma$ has a very significant effect
on the stationary probability density of the system coordinate $r_{z}$.
Increasing $\gamma$ concentrates the pdf closer to the boundaries and in an
asymmetric manner due to the system Hamiltonian bias $\epsilon$. By evolving
$\gamma$ from unity to a larger value, we can capture the effect of a
measurement of the system energy, as this would cause individual trajectories
to dwell in the vicinity of the eigenstates of $H_{sys}$ at $r_{z}=\pm 1$, or
equivalently at system energies $\text{Tr}\left(H_{sys}\rho\right)=\epsilon
r_{z}$ of $\pm\epsilon$.

###  3.3 Average system entropy production and boundary contributions

Using its definition in Eq. (11), the incremental stochastic system entropy
production over each time step $dt$ can be calculated from the time-dependent
pdf $p(r_{z},t)$ obtained from the Fokker-Planck equation (50). We then
average over many realisations of the behaviour. However, we have noticed in
practice that combining the average system entropy production with the average
environmental entropy production Eq. (10), obtained numerically for time
independent environmental coupling, does not yield a vanishing total entropy
production that we would expect for a system in a stationary state with zero
probability current. We have therefore investigated this situation further.
The incremental system entropy production averaged over $r_{z}$ and the noises
for our system may be written as

$\displaystyle d\langle\langle\Delta s_{sys}\rangle\rangle$
$\displaystyle=\int_{-1}^{1}d\langle\Delta s_{sys}\rangle\,p(r_{z},t)dr_{z},$
(54)

where $d\langle\Delta s_{sys}\rangle$ is the noise averaged increment in
system entropy production, which is a function of $r_{z}$ and $t$. By
integrating by parts twice and retaining the boundary terms, see Appendix A.3
for details, we obtain:

$\displaystyle d\langle\langle\Delta s_{sys}\rangle\rangle$
$\displaystyle=dS_{G}-\left[J_{z}\ln
p(r_{z},t)\right]_{-1}^{1}dt-\left[D_{zz}\frac{\partial p(r_{z},t)}{\partial
r_{z}}\right]_{-1}^{1}dt,$ (55)

where $dS_{G}$ is the increment in the Gibbs entropy of the system, shown in
Eq. (13), and $J_{z}$ is the $r_{z}$-component of the probability current
defined in Appendix A.3.

In Eq. (13) we noted that $d\langle\langle\Delta s_{sys}\rangle\rangle=dS_{G}$
held only under certain conditions, now shown to be where the boundary terms
in Eq. (55) vanish. For example, the current $J_{z}$ and often the pdf and its
derivatives vanish at the boundaries of the parameter space, and diffusion
coefficients can also go to zero. For our system, although $D_{zz}=0$ at
$r_{z}=\pm 1$, it can be shown that the gradient $\partial
p_{st}(r_{z},t)/\partial r_{z}$ is singular at the boundaries for the
stationary solution of the Fokker-Planck equation, and that by implication
this applies to nonstationary situations as well. There is therefore a need to
correct the usual result $d\langle\langle\Delta s_{sys}\rangle\rangle=dS_{G}$.

In contrast to this analysis, we observe a mean system entropy production
consistent with zero when our system is in a stationary state, which suggests
that the numerically obtained time dependent pdf is not sufficiently accurate
near the boundaries, or the sampling of these regions by the generated
trajectories is too limited, to produce the appropriate additional
contributions. Our approach, therefore, is to estimate the boundary terms
analytically for the stationary pdf and to add them to the numerical results.

Using the stationary pdfs given in Eqs. (52) and (51), together with the
diffusion coefficient $D_{zz}$ given in Eq. (46), it may be shown that the
boundary correction with $J_{z}$ in Eq. (55) vanishes at equilibrium, and so
does $dS_{G}.$ This allows us to write the mean stationary system entropy
production increments for general $\gamma$ as

$d\langle\langle\Delta s_{sys}\rangle\rangle_{st}=-\left[D_{zz}\frac{\partial
p_{st}(r_{z})_{\gamma}}{\partial r_{z}}\right]_{-1}^{1}dt,$ (56)

which does not vanish. In fact, it is possible to obtain a simple analytical
expression for $\gamma=1$ to lowest order in $\beta\epsilon$:

$d\langle\langle\Delta
s_{sys}\rangle\rangle_{st}^{\gamma=1}=-\lambda^{2}(\beta\epsilon)^{2}dt.$ (57)

This shows that this boundary correction can be obtained exactly analytically
for $\gamma=1$, though it becomes harder for $\gamma\neq 1$. In the latter
case we have a more complicated form of $p_{st}(r_{z})$, and hence it is
necessary to expand it analytically in terms of small $\beta\epsilon$ before
calculating the boundary correction.

###  3.4 Average environmental entropy production

With the choice of the couplings and Hamiltonian, we can calculate the
environmental entropy production from Eq. (11). We need to determine
$\boldsymbol{A}^{irr}$ and $\boldsymbol{A}^{rev}$ from Eqs. (42) and (43).
This is done by understanding how the components of the coherence vector and
the various contributions to $\boldsymbol{A}$ transform under time reversal.
For a system without spin degrees of freedom, the time reversal operator is
given by $\Theta=K$ , where $K$ is the operation of complex conjugation [59]
and the time reversal operation on the Pauli matrices produces

$\Theta\sigma_{z}\Theta^{-1}=\sigma_{z}\qquad\Theta\sigma_{y}\Theta^{-1}=-\sigma_{y}\qquad\Theta\sigma_{x}\Theta^{-1}=\sigma_{x}.$
(58)

Thus $r_{x}$ and $r_{z}$ are even under time reversal, and $r_{y}$ is odd, and
hence $\phi$ is also odd, and we can separate the coefficients of the $dt$
terms in Eqs. (42) and (43) into their $\boldsymbol{A}^{irr}$ and
$\boldsymbol{A}^{rev}$components. We can write

$\boldsymbol{A}^{irr}=\left(\begin{matrix}-4\lambda^{2}\left(\beta\epsilon+r_{z}\right)\\\
0\end{matrix}\right),\hskip
18.00003pt\boldsymbol{A}^{rev}=\left(\begin{matrix}0\\\
2\epsilon\end{matrix}\right)\,,$ (59)

and obtain an expression for the environmental entropy production using Eq.
(10), since we already have the form of $\boldsymbol{D}$ from Eq. (46).

##  4 Simulations and Results

In this section we describe two protocols designed to represent the dynamics
of connecting the system to and disconnecting it from a measuring device,
respectively, and we compute the associated stochastic entropy production for
each.

###  4.1 The computational procedure

In order to demonstrate adherence to the detailed fluctuation theorem we
consider two simple protocols that are a time reversal of each other:

  * •

Connecting a measuring device (protocol M): begin at $t=0$ with the system
thermalised for $\gamma=1$ and the initial pdf defined as in Figure 2; for
$t>0$ perform simulations of Eqs. (42) and (43) using $\gamma=2$;

  * •

Disconnecting a measuring device (protocol $\overline{{\rm M}}$): begin at
$t=0$ with the system thermalised for $\gamma=2$ and the initial pdf defined
as in Figure 2; then, for $t>0$ perform simulations of Eqs. (42) and (43)
using $\gamma=1$.

The thermalisation of the system is defined by the pdf given by the stationary
Fokker-Planck equation as plotted in Figure 2. For both protocols, the initial
value of $r_{z}$ for a trajectory is randomly sampled from the appropriate pdf
$p_{st}(r_{z})$ corresponding to either $\gamma=1$ or $\gamma=2$, while the
value of $\phi$ is randomly sampled from a uniform distribution within the
range $[0,2\pi)$. Whenever we specify a particular value of $\gamma$, we will
be referring to the value used for the dynamics (i.e., for $t>0$), unless we
explicitly mention the thermalisation condition ($t=0$). For both protocols
several parameters of the system are kept constant: the environmental inverse
temperature $\beta=0.1$, the Hamiltonian parameter $\epsilon=1$, the
proportionality constant in the spectral density $\alpha=0.01$, and
$\lambda=\sqrt{0.2}$.

In order to calculate the system entropy production, we solve the Fokker-
Planck equation (50) for $r_{z}$ and observe the behaviour of the system going
from a stationary state at a particular value of $\gamma$ to another
stationary state defined by a different $\gamma$. The protocols will drive the
system between two states obtained from the Fokker-Planck equation at
different $\gamma$ values. The evolution from one stationary pdf to the other
occurs over a timescale of order $t=1$, for both protocols.

In summary, the computational procedure is as follows:

  1. 1.

Obtain the solution $p(r_{z},t)$ of the Fokker-Planck equation (50) using the
appropriate stationary pdf as the initial condition (at $t=0$) for the
trajectories; in both cases the boundary conditions at $r_{z}=\pm 1$ are
chosen corresponding to a zero probability current $J_{z}$ defined in Appendix
A.3. Note that each pdf is obtained on a grid of $r_{z}$ values.

  2. 2.

Start a loop over the noises (independently for each protocol):

    1. (a)

generate noises $dW_{x}$, $dW_{y}$ and $dW_{z}$ for each time step $dt$;

    2. (b)

run the stochastic dynamics for $r_{z}(t)$ and $\phi(t)$ for both protocols
using Eqs. (42) and (43) based on initial values sampled from the
corresponding $\gamma$-dependent stationary pdf;

    3. (c)

for each time increment $dt$ and using obtained values of $r_{z}(t)$ and
$\phi(t)$, calculate the incremental contributions to the environmental,
$d\Delta s_{env}$, and system, $d\Delta
s_{sys}=-d\left(\ln\,p(r_{z},t)\right)$, stochastic entropy productions via
respectively, Eq. (10) and the first line of Eq. (11); in the latter case the
increment of the logarithm of the pdf is calculated as a difference between
its values at two consecutive times. The values of the pdf $p(r_{z},t)$ for
these times at the required value of $r_{z}$ are obtained by linearly
interpolating each pdf between the two nearest $r_{z}$ values on the grid;

    4. (d)

the incremental total stochastic entropy production, at each time step $dt$,
is calculated as a sum, $d\Delta s_{tot}=d\Delta s_{env}+d\Delta s_{sys}$;
accumulate these contributions over an entire trajectory;

    5. (e)

go back to step 2(a) to run another trajectory; the calculation is repeated
the required number of times using different manifestations of the noises. The
stochastic entropy production is averaged over the trajectories. Each protocol
will have its corresponding independent ensemble of entropy productions.

  3. 3.

Once the averaged values of the entropy production increments are obtained
from a large ensemble of trajectories, the boundary term from Eq. (55) is
added to the ensemble average $d\left\langle\left\langle\Delta
s_{tot}\right\rangle\right\rangle$ as

$d\langle\langle\Delta s_{tot}\rangle\rangle=d\langle\langle\Delta
s_{env}\rangle\rangle+d\langle\langle\Delta
s_{sys}\rangle\rangle-\left[D_{zz}\frac{\partial p_{st}(r_{z})}{\partial
r_{z}}\right]_{-1}^{1}.$ (60)

for the specific value of $\gamma$. This yields the final increment of the
mean total entropy production $d\langle\langle\Delta s_{tot}\rangle\rangle$
for the given time step $dt$, separately for the two protocols.

###  4.2 Dynamics of $r_{z}$ and $\phi$

The thermal state $\bar{\rho}_{eq}$ of the system is approximately
proportional to $e^{-\beta H_{sys}}$ for weak system-environment coupling. In
the derivation of Eq. (32) we have assumed $\beta H_{sys}\ll 1,$ so we can
write $\bar{\rho}_{eq}\approx\frac{1}{2}(\mathbb{I}-\beta H_{sys})$ and hence
obtain the average of $r_{z}$ as $\bar{r}_{z}={\rm
Tr}\left(\sigma_{z}\bar{\rho}_{eq}\right)=-\beta\epsilon\approx-0.1$ for our
parameters. We find that our numerics supports this: by running one million
realisations of the dynamics for protocol M (representing the connection of
the measuring device), we observe that the averaged results of $r_{z}$ do
indeed remain constant throughout the simulation. Furthermore, the dynamics of
$\phi$ maintains a pdf that is constant in $\phi$ which matches the
assumptions used to derive Eq. (50). While the averaged system behaviour
remains unchanged, the stochastic trajectories differ significantly due to
measurements, shown in Figure 3. Note that each individual trajectory, if run
for a sufficiently long time, would dwell in the vicinity of either of the
boundaries $r_{z}=\pm 1$, jumping between the two. However, once the system
reaches stationarity (equilibrium), the ensemble of trajectories reaches the
stationary distribution $p_{st}(r_{z})_{\gamma}$ associated with the
corresponding value of $\gamma$. One would expect this distribution to match
the distribution following from solving the Fokker-Planck equation from
Section 3.2. As can be seen from Figure 3 (right panel), this is indeed the
case in our simulations.

Figure 3: Left panel: a set of 10 random trajectories that were initialised at
$r_{z}=0$ for protocol M ($\gamma=2$), displaying the drifting behaviour of
the system towards the eigenstates driven by the measuring device. Right
panel: a normalised histogram for protocol M based on one million individual
trajectories at $t=2,$ and compared with $p_{st}(r_{z})_{\gamma=2}$ of Eq.
(51) obtained by solving the Fokker-Planck equation.

This demonstrates the stochastic behaviour of the system and the effect of the
measurement interactions that send $r_{z}$ to the vicinity of the eigenstates,
while maintaining the same thermal averaged state throughout the evolution.
This exemplifies how crucial it is to represent the stochastic behaviour of
quantum measurements, and to recognise that all physical changes in the system
are stochastic in nature. We then show in Figure 4 how the system explores the
surface of the Bloch sphere biased by the Hamiltonian and under the influence
of the environment for a couple of individual trajectories for protocols M and
$\overline{{\rm M}}$ (connecting and disconnecting the measurement device).
The trajectory with $\gamma=1$ explores the Bloch sphere fairly uniformly as
seen in Figure 2, but for $\gamma=2$, the trajectories tend to dwell near the
boundaries at $r_{z}=\pm 1$.

Figure 4: Example individual trajectories for the two protocols ($\gamma=2$,
blue; $\gamma=1$, red) , exploring the Bloch sphere. Parameters used are:
$\beta=0.1$, $\epsilon=1$, $t_{max}=2$, $dt=10^{-5}$, and $\alpha=0.01$.
Initial states for the protocols were randomly sampled from the stationary
solutions to the Fokker-Planck equation, for $\gamma=1$ and $\gamma=2$,
respectively.

###  4.3 Entropy production

The total stochastic entropy production is obtained by summing up the
stochastic system entropy production in Eq. (11), the stochastic environmental
entropy production in Eq. (10), and (for the ensemble average) the boundary
correction term from Eq. (55). The results for protocol M with $\gamma=2$ are
shown in Figure 5.

Figure 5: Stochastic entropy production for protocol M with $\gamma=2$. The
noise and $r_{z}$ averaged total entropy production is shown with (solid
black) and without (dotted black) the boundary correction. The red curve
corresponds to a single stochastic realisation. The dashed blue line denotes
the asymptotic analytical average total stochastic entropy production for the
process from Eq. (12). The averages result from one million realisations for
$\beta=0.1$, $\epsilon=1$, $t_{max}=2$, $dt=10^{-5}$, and $\alpha=0.01$.

As is seen from the Figure, without the boundary correction the mean entropy
production continues to increase which is inconsistent with a stationary
system. The need for boundary corrections is probably a consequence of our
choice of system variables: it arises from the singular pdf density at the
boundaries.

Figure 6: The black curve corresponds to noise- and $r_{z}$-averaged total
stochastic entropy production for protocol $\overline{{\rm M}}$, with
$\gamma=1$, including the boundary correction. The red curve is a single
stochastic realisation, and the dashed blue curve is the asymptotic analytical
averaged total stochastic entropy production from Eq. (12). The averages are
obtained from one million realisations for $\beta=0.1$, $\epsilon=1$,
$t_{max}=2$, $dt=10^{-5}$, and $\alpha=0.01$. The boundary correction for
$\gamma=1$ is significantly smaller than that for $\gamma=2$.

The entropy production for protocol $\overline{{\rm M}}$ with $\gamma=1$ is
shown in Figure 6. As with protocol M, these results also lead to a ceiling in
the average total stochastic entropy production once the system equilibrates.
We see that the process of disconnecting the measurement device also leads to
a positive mean entropy production; in faimagesct, the mean entropy production
is higher than for the connection process. Note that entropy production tends
to a constant asymptotically even for an individual trajectory. Both Figs. 5
and 6 display perfect alignment between the analytically averaged total
entropy production for late times, and its numerical counterpart. This
validates the approximations we have made, and the form we have used for the
boundary correction terms.

###  4.4 Entropy production pdf and detailed fluctuation theorem

As a check on the accuracy of the entropy production results, we verify that
the detailed fluctuation theorem Eq. (2) is satisfied. The detailed
fluctuation relation concerns processes described by forward and reverse
protocols that are time reversals of one another, and where the final pdf
under the forward protocol is the same as the initial pdf under the reverse
protocol. However, as we deal here with _distributions_ of the entropy
production calculated over _many_ trajectories, the strict reversal condition
of the trajectories in protocols M and $\overline{\text{M}}$ is not necessary
and we can compare $e^{\Delta s_{tot}}$ directly with the ratio of entropy
pdfs $P^{M}\left(\Delta s_{tot}\right)/P^{\overline{M}}\left(-\Delta
s_{tot}\right)$ . Figure 7 displays the asymptotic distributions of total
stochastic entropy production for protocols M and $\overline{{\rm M}}$. These
distributions should satisfy the detailed fluctuation relation, and we indeed
find a very good adherence except for deviations at the extremes of the range
due to insufficient sampling.

Figure 7: Histograms displaying the asymptotic values of the total stochastic
entropy production $\Delta s_{tot}$ for protocol M (blue; $\gamma=2$) and
protocol $\overline{{\rm M}}$ (green; $\gamma=1$) based on one million
individual trajectories (realisations). Black lines show both sides of Eq.
(2), revealing the adherence to the detailed fluctuation theorem (2).

##  5 Discussion and Conclusions

In this paper, we have considered the stochastic entropy production induced by
continuous measurements of a simple open quantum system. We begin with the
stochastic Liouville-von Neumann (SLN) equation for the system dynamics, in
which the environments are represented by coloured noises, taking the
appropriate stochastic average with respect to these noises and retaining the
non-Markovian generality of the dynamics. Trajectories of the reduced density
matrix generated by the SLN equation are not physical, only the averaged
evolution is of physical significance. In contrast, we require a stochastic
dynamical unravelling that can be interpreted as the outcome of a physical
weak continuous measurement of the system in conjunction with thermalising
influences of the environment. To do so, we take the Markovian limit of the
behaviour of the environment in the SLN by assuming a high temperature with a
specific choice of a two-level bosonic Hamiltonian
$H_{sys}=\epsilon\sigma_{z}$. These steps enable us to perform a Markovian
stochastic unravelling described by an Itô process, which then allows us to
derive an appropriate Fokker-Planck equation and understand how the pdf of the
reduced density matrix changes with the strength of the continuous
measurements while coupled to a thermalising bath. We have set up a system
that interacts with three independent harmonic baths through the three Pauli
matrices, allowing for a quasi-isotropic stochastic exploration of the Bloch
sphere consistent with residence in a thermal state. We take the view that a
system under constant interaction with an underspecified environment is best
represented by an ensemble of reduced density matrices that conveys its
uncertainty.

The measurement of system energy is then realised by increasing $\gamma$, the
strength of environmental coupling to $\sigma_{z}$, causing the system to
dwell in the vicinity of one of the two energy eigenstates. We can then obtain
the stochastic entropy production associated with each stochastic trajectory
of the system. Note that if the pdf of the reduced density matrix is
stationary for a particular set of environmental couplings, raising $\gamma$
will cause the pdf to increase in the vicinity of the eigenstates of the
system Hamiltonian but this will not change the _average_ energy of an
ensemble of systems. A model of the dynamics and thermodynamics of measurement
requires the system to be represented by a _member_ of an ensemble: the
average behaviour will not suffice.

To calculate the entropy production associated with measurement, we use an
analysis of the Markovian system dynamics developed in [5]. This contrasts
with the calculation of the stochastic entropy production for quantum systems
based on forward and reverse trajectories constructed using forward and
reverse Kraus operators [31, 36, 28]. We also compute system entropy
production by considering the evolution of the pdf of the reduced density
matrix of the system.

We find that the numerically calculated entropy production contains some
subtle numerical artefacts. These are corrected by calculating the mean
analytical system entropy production in a stationary state, allowing us to
obtain boundary correction terms that eliminate the artefacts.

We find that the increase and decrease in coupling strength $\gamma$,
corresponding to the operation of a measuring device, are both accompanied by
a positive mean stochastic entropy production. This is the entropic cost of
quantum measurement. Furthermore, we show that the stochastic entropy
production associated with a quantum measurement is finite as the system
achieves a new stationary state. Finally, we have shown that the processes of
attachment and detachment of a measuring device generate distributions of
stochastic entropy production that satisfy a detailed fluctuation relation. We
aim to extend this approach next to more complicated systems, to include non-
Markovian dynamics and to relax the requirement for a high temperature
approximation.

## Acknowledgments

D. Matos is supported by the EPSRC Centre for Doctoral Training in Cross-
Disciplinary Approaches to Non-Equilibrium Systems (CANES, Grant No.
EP/L015854/1). Calculations in this paper were performed using the King’s
College HPC cluster Gravity.

##  Appendix A Appendices

###  A.1 Averaged SLN equation

Here we shall consider the calculation of the functional derivatives of
$\rho_{s}$ with respect to the noises needed for the derivation of Eq. (22).
The formal solution of the SLN equation is given by

$\rho_{s}(t)=U_{+}(t,0)\rho_{s}(0)U_{-}(0,t),$ (61)

where $U_{+}$ and $U_{-}$ are the appropriate forward and backward propagators
defined as

$\displaystyle U_{+}(t,0)$ $\displaystyle=$
$\displaystyle\widehat{\mathcal{T}}_{+}\exp\left[-i\int_{0}^{t}dt^{\prime}H_{+}(t^{\prime})\right],$
(62) $\displaystyle U_{-}(0,t)$ $\displaystyle=$
$\displaystyle\widehat{\mathcal{T}}_{-}\exp\left[i\int_{0}^{t}dt^{\prime}H_{-}(t^{\prime})\right],$
(63)

where $\widehat{\mathcal{T}}_{+}$ and $\widehat{\mathcal{T}}_{-}$ are the
forward and backward time ordering operators, respectively, and

$H_{\pm}(t)=H_{sys}(t)-\sum_{k}\left(\eta_{k}(t)\pm\frac{\nu_{k}(t)}{2}\right)f_{k}$
(64)

are the Hamiltonian operators. Note that they are not adjoints of each other,
see Ref. [11] for details. Hence the functional derivative of the density
matrix can be written as

$\frac{\delta\rho_{s}(t)}{\delta\zeta(t^{\prime})}=\frac{\delta
U_{+}(t,0)}{\delta\zeta(t^{\prime})}\rho_{s}(0)U_{-}(0,t)+U_{+}(t,0)\rho_{s}(0)\frac{\delta
U_{-}(0,t)}{\delta\zeta(t^{\prime})}.$ (65)

The derivatives of the forward and backward propagators with respect to
arbitrary perturbations of the noises are given by:

$\displaystyle\frac{\delta U_{+}(t,t^{\prime})}{\delta\zeta(\tau^{\prime})}$
$\displaystyle=-i\int_{t^{\prime}}^{t}U_{+}(t,\tau)\frac{\delta
V_{+}(\tau)}{\delta\zeta(\tau^{\prime})}U_{+}(\tau,t^{\prime})d\tau$ (66)
$\displaystyle\frac{\delta U_{-}(t^{\prime},t)}{\delta\zeta(\tau^{\prime})}$
$\displaystyle=i\int_{t^{\prime}}^{t}U_{-}(t^{\prime},\tau)\frac{\delta
V_{-}(\tau)}{\delta\zeta(\tau^{\prime})}U_{-}(\tau,t)d\tau,$ (67)

where $V_{\pm}(\tau)=-\left(\eta(\tau)\pm\frac{\nu(\tau)}{2}\right)f$, for
convenience dropping the suffix $k.$ We can now define the variation of the
propagators with respect to the noises $\eta(t)$ and $\nu(t)$. Starting with
$\eta(t)$:

$\displaystyle\frac{\delta U_{+}(t,t^{\prime})}{\delta\eta(\tau)}$
$\displaystyle=-i\int_{t^{\prime}}^{t}U_{+}(t,\tau^{\prime})\left[-\delta(\tau-\tau^{\prime})f(\tau^{\prime})\right]U_{+}(\tau^{\prime},t^{\prime})d\tau^{\prime}$
(68) $\displaystyle=iU_{+}(t,\tau)f(\tau)U_{+}(\tau,t^{\prime})\,,$
$\displaystyle\frac{\delta U_{-}(t^{\prime},t)}{\delta\eta(\tau)}$
$\displaystyle=i\int_{t^{\prime}}^{t}U_{-}(t^{\prime},\tau^{\prime})\left[-\delta(\tau-\tau^{\prime})f(\tau^{\prime})\right]U_{-}(\tau^{\prime},t)d\tau^{\prime}$
$\displaystyle=-iU_{-}(t^{\prime},\tau)f(\tau)U_{-}(\tau,t).$

For $\nu(t)$:

$\displaystyle\frac{\delta U_{+}(t,t^{\prime})}{\delta\nu(\tau)}$
$\displaystyle=-i\int_{t^{\prime}}^{t}U_{+}(t,\tau^{\prime})\left[-\frac{1}{2}\delta(\tau-\tau^{\prime})f(\tau^{\prime})\right]U_{+}(\tau^{\prime},t^{\prime})d\tau^{\prime}$
(69) $\displaystyle=\frac{i}{2}U_{+}(t,\tau)f(\tau)U_{+}(\tau,t^{\prime})\,,$
$\displaystyle\frac{\delta U_{-}(t^{\prime},t)}{\delta\nu(\tau)}$
$\displaystyle=i\int_{t^{\prime}}^{t}U_{-}(t^{\prime},\tau^{\prime})\left[\frac{1}{2}\delta(\tau-\tau^{\prime})f(\tau^{\prime})\right]U_{-}(\tau^{\prime},t)d\tau^{\prime}$
$\displaystyle=\frac{i}{2}U_{-}(t^{\prime},\tau)f(\tau)U_{-}(\tau,t).$

This then allows us to write $\left\langle\eta(t)\rho_{s}(t)\right\rangle$ as

$\displaystyle\langle\eta(t)\rho_{s}(t)\rangle$
$\displaystyle=\int_{0}^{t}dt^{\prime}K^{\eta\eta}(t-t^{\prime})\left\langle\frac{\delta\rho_{s}(t)}{\delta\eta(t^{\prime})}\right\rangle+\int_{0}^{t}dt^{\prime}K^{\eta\nu}(t-t^{\prime})\left\langle\frac{\delta\rho_{s}(t)}{\delta\nu(t^{\prime})}\right\rangle$
$\displaystyle=\int_{0}^{t}dt^{\prime}K^{\eta\eta}(t-t^{\prime})\left\langle
iU_{+}(t,t^{\prime})f(t^{\prime})U_{+}(t^{\prime},0)\rho_{s}(0)U_{-}(0,t)-U_{+}(t,0)\rho_{s}(0)iU_{-}(0,t^{\prime})f(t^{\prime})U_{-}(t^{\prime},t)\right\rangle$
(70)
$\displaystyle\qquad+\int_{0}^{t}dt^{\prime}K^{\eta\nu}(t-t^{\prime})\left\langle\frac{i}{2}U_{+}(t,t^{\prime})f(t^{\prime})U_{+}(t^{\prime},0)\rho_{s}(0)U_{-}(0,t)+U_{+}(t,0)\rho_{s}(0)\frac{i}{2}U_{-}(0,t^{\prime})f(t^{\prime})U_{-}(t^{\prime},t)\right\rangle.$

Using this expression, we arrive at Eq. (22) given in the text.

###  A.2 Cylindrical coordinate equations of motion

In this section, we derive the equations of motion in cylindrical coordinates
$(r^{2},r_{z},\phi)$ for general $\gamma$. With these coordinates, the $x$ and
$y$ components are given by

$r_{x}=\sqrt{r^{2}-r_{z}^{2}}\cos\phi\hskip
45.00008ptr_{y}=\sqrt{r^{2}-r_{z}^{2}}\sin\phi.$ (71)

From these expressions and using Itô’s Lemma, we can obtain the equations of
motion for this set of coordinates:

$\displaystyle dr^{2}$ $\displaystyle=\frac{\partial r^{2}}{\partial
r_{z}}dr_{z}+\frac{1}{2}\frac{\partial^{2}r^{2}}{\partial
r_{z}^{2}}\left(dr_{z}\right)^{2}+\frac{\partial r^{2}}{\partial
r_{y}}dr_{y}+\frac{1}{2}\frac{\partial r^{2}}{\partial
r_{y}^{2}}\left(dr_{y}\right)^{2}+\frac{\partial r^{2}}{\partial
r_{x}}dr_{x}+\frac{1}{2}\frac{\partial^{2}r^{2}}{\partial
r_{x}^{2}}\left(dr_{x}\right)^{2}$ (72) $\displaystyle d\phi$
$\displaystyle=\frac{\partial\phi}{\partial
r_{x}}dr_{x}+\frac{1}{2}\frac{\partial^{2}\phi}{\partial
r_{x}^{2}}\left(dr_{x}\right)^{2}+\frac{\partial\phi}{\partial
r_{y}}dr_{y}+\frac{1}{2}\frac{\partial^{2}\phi}{\partial
r_{y}^{2}}\left(dr_{y}\right)^{2}+\frac{\partial^{2}\phi}{\partial
r_{x}\partial r_{y}}dr_{x}dr_{y}.$ (73)

The required derivatives are:

$\displaystyle\frac{\partial r^{2}}{\partial r_{x}}=2r_{x}\qquad\frac{\partial
r^{2}}{\partial r_{y}}=2r_{y}\qquad\frac{\partial r^{2}}{\partial
r_{z}}=2r_{z}$ (74) $\displaystyle\frac{\partial^{2}r^{2}}{\partial
r_{x}^{2}}=\frac{\partial^{2}r^{2}}{\partial
r_{y}^{2}}=\frac{\partial^{2}r^{2}}{\partial r_{z}^{2}}=2$ (75)
$\displaystyle\frac{\partial\phi}{\partial
r_{x}}=-\frac{r_{y}}{r_{x}^{2}+r_{y}^{2}}\qquad\frac{\partial\phi}{\partial
r_{y}}=\frac{r_{x}}{r_{x}^{2}+r_{y}^{2}}$ (76)
$\displaystyle\frac{\partial^{2}\phi}{\partial
r_{x}^{2}}=\frac{2r_{x}r_{y}}{\left(r_{x}^{2}+r_{y}^{2}\right)^{2}}\qquad\frac{\partial^{2}\phi}{\partial
r_{y}^{2}}=-\frac{2r_{x}r_{y}}{\left(r_{x}^{2}+r_{y}^{2}\right)^{2}}\qquad\frac{\partial^{2}\phi}{\partial
r_{x}r_{y}}=\frac{r_{y}^{2}-r_{x}^{2}}{\left(r_{x}^{2}+r_{y}^{2}\right)^{2}},$
(77)

which allows us to simplify Eqs. (72) and (73):

$\displaystyle dr^{2}$
$\displaystyle=2r_{z}dr_{z}+\left(dr_{z}\right)^{2}+2r_{y}dr_{y}+\left(dr_{y}\right)^{2}+2r_{x}dr_{x}+\left(dr_{x}\right)^{2}$
(78) $\displaystyle d\phi$
$\displaystyle=-\frac{r_{y}}{r_{x}^{2}+r_{y}^{2}}dr_{x}+\frac{r_{x}r_{y}}{\left(r_{x}^{2}+r_{y}^{2}\right)^{2}}\left(dr_{x}\right)^{2}+\frac{r_{x}}{r_{x}^{2}+r_{y}^{2}}dr_{y}-\frac{r_{x}r_{y}}{\left(r_{x}^{2}+r_{y}^{2}\right)^{2}}\left(dr_{y}\right)^{2}+\frac{r_{y}^{2}-r_{x}^{2}}{\left(r_{x}^{2}+r_{y}^{2}\right)^{2}}dr_{x}dr_{y}.$
(79)

We can calculate $(dr_{x})^{2}$, $(dr_{y})^{2}$, $(dr_{z})^{2}$ and
$dr_{x}dr_{y}$ using Eqs. (35)-(37) and the substitution
$dW_{i}dW_{j}=\delta_{ij}dt$, keeping linear terms in $dt$:

$\displaystyle(dr_{x})^{2}$ $\displaystyle=\lambda^{2}\left(\beta\epsilon
r_{z}+2(1-r_{x}^{2})\right)^{2}dt+4\lambda^{2}r_{x}^{2}r_{y}^{2}dt+4\lambda^{2}\gamma^{2}r_{x}^{2}r_{z}^{2}dt$
(80) $\displaystyle(dr_{y})^{2}$
$\displaystyle=4\lambda^{2}r_{x}^{2}r_{y}^{2}dt+\lambda^{2}\left(\beta\epsilon
r_{z}+2(1-r_{y}^{2})\right)^{2}dt+4\lambda^{2}\gamma^{2}r_{y}^{2}r_{z}^{2}dt$
(81) $\displaystyle(dr_{z})^{2}$
$\displaystyle=\lambda^{2}r_{x}^{2}\left(\beta\epsilon+2r_{z}\right)^{2}dt+\lambda^{2}r_{y}^{2}\left(\beta\epsilon+2r_{z}\right)^{2}dt+4\lambda^{2}\gamma^{2}\left(1-r_{z}^{2}\right)^{2}dt$
(82) $\displaystyle dr_{x}dr_{y}$
$\displaystyle=-2\lambda^{2}r_{x}r_{y}\left(\beta\epsilon
r_{z}+2(1-r_{x}^{2})\right)dt-2\lambda^{2}r_{x}r_{y}\left(\beta\epsilon
r_{z}+2(1-r_{y}^{2})\right)dt+4\lambda^{2}\gamma^{2}r_{x}r_{y}r_{z}^{2}dt.$
(83)

After some algebra and dropping quadratic terms in $\beta$, we obtain

$\displaystyle(dr_{x})^{2}$
$\displaystyle=4\lambda^{2}\left[r_{x}^{2}\left(r_{x}^{2}+r_{y}^{2}+\gamma^{2}r_{z}^{2}\right)+\beta\epsilon
r_{z}\left(1-r_{x}^{2}\right)+1-2r_{x}^{2}\right]dt$
$\displaystyle(dr_{y})^{2}$
$\displaystyle=4\lambda^{2}\left[r_{y}^{2}\left(r_{x}^{2}+r_{y}^{2}+\gamma^{2}r_{z}^{2}\right)+\beta\epsilon
r_{z}\left(1-r_{y}^{2}\right)+1-2r_{y}^{2}\right]dt$ (84)
$\displaystyle(dr_{z})^{2}$ $\displaystyle=4\lambda^{2}\left[\beta\epsilon
r_{z}\left(r_{x}^{2}+r_{y}^{2}\right)+r_{z}^{2}\left(r_{x}^{2}+r_{y}^{2}\right)-2\gamma^{2}r_{z}^{2}+\gamma^{2}\left(1+r_{z}^{4}\right)\right]dt.$

The sum of Eqs. (84) is given by

$\displaystyle(dr_{x})^{2}+(dr_{y})^{2}+(dr_{z})^{2}$
$\displaystyle=4\lambda^{2}\left[-2\gamma^{2}r_{z}^{2}+2\beta\epsilon
r_{z}+\gamma^{2}\left(1+r_{z}^{4}\right)\right.$ (85)
$\displaystyle\left.\qquad\qquad\qquad+\left(r^{2}-r_{z}^{2}\right)\left(r^{2}+\gamma^{2}r_{z}^{2}-2\right)+2\right]dt.$

Moving on to the $2r_{i}dr_{i}$ terms:

$\displaystyle 2r_{x}dr_{x}$ $\displaystyle=-4\epsilon
r_{x}r_{y}dt-4\lambda^{2}\left(1+\gamma^{2}\right)r_{x}^{2}dt+2\lambda
r_{x}\left(\beta\epsilon r_{z}+2\left(1-r_{x}^{2}\right)\right)dW_{x}-4\lambda
r_{x}^{2}r_{y}dW_{y}-4\gamma\lambda r_{x}^{2}r_{z}dW_{z}$ $\displaystyle
2r_{y}dr_{y}$ $\displaystyle=4\epsilon
r_{x}r_{y}dt-4\lambda^{2}\left(1+\gamma^{2}\right)r_{y}^{2}dt-4\lambda
r_{x}r_{y}^{2}dW_{x}+2\lambda r_{y}\left(\beta\epsilon
r_{z}+2\left(1-r_{y}^{2}\right)\right)dW_{y}-4\gamma\lambda
r_{y}^{2}r_{z}dW_{z}$ (86) $\displaystyle 2r_{z}dr_{z}$
$\displaystyle=-8\lambda^{2}r_{z}^{2}dt-8\lambda^{2}\beta\epsilon
r_{z}dt-2\lambda r_{x}r_{z}\left(\beta\epsilon+2r_{z}\right)dW_{x}-2\lambda
r_{y}r_{z}\left(\beta\epsilon+2r_{z}\right)dW_{y}+4\gamma\lambda
r_{z}\left(1-r_{z}^{2}\right)dW_{z}.$

Adding up Eqs. (86) leads to

$\displaystyle 2r_{x}dr_{x}+2r_{y}dr_{y}+2r_{z}dr_{z}$
$\displaystyle=-4\lambda^{2}\left(2r_{z}^{2}+(1+\gamma^{2})(r^{2}-r_{z}^{2})\right)dt-8\lambda^{2}\beta\epsilon
r_{z}dt+4\lambda r_{x}\left(1-r^{2}\right)dW_{x}$ (87)
$\displaystyle\qquad+4\lambda r_{y}\left(1-r^{2}\right)dW_{y}+4\lambda\gamma
r_{z}\left(1-r^{2}\right)dW_{z}.$

By adding Eqs. (85) and (87) we obtain the final expression for $dr^{2}$:

$\displaystyle\begin{aligned}
dr^{2}&=4\lambda^{2}\left(r^{2}-1\right)\left(\gamma^{2}r_{z}^{2}-\gamma^{2}+r^{2}-r_{z}^{2}-2\right)dt+4\lambda\sqrt{1-\frac{r_{z}^{2}}{r^{2}}}r\cos\phi\left(1-r^{2}\right)dW_{x}\\\
&+4\lambda\sqrt{1-\frac{r_{z}^{2}}{r^{2}}}r\sin\phi\left(1-r^{2}\right)dW_{y}+4\lambda\gamma
r_{z}\left(1-r^{2}\right)dW_{z}\,.\end{aligned}$ (88)

Moving on to calculating $d\phi$, the contributions linear in $dr_{i}$ are:

$-\frac{r_{y}}{r_{x}^{2}+r_{y}^{2}}dr_{x}+\frac{r_{x}}{r_{x}^{2}+r_{y}^{2}}dr_{y}=\frac{1}{r_{x}^{2}+r_{y}^{2}}\left[2\epsilon\left(r_{x}^{2}+r_{y}^{2}\right)dt-\lambda
r_{y}\left(\beta\epsilon r_{z}+2\right)dW_{x}+\lambda r_{x}\left(\beta\epsilon
r_{z}+2\right)dW_{y}\right].$ (89)

The terms quadratic in $dr_{i}$ lead to

$\frac{r_{x}r_{y}}{\left(r_{x}^{2}+r_{y}^{2}\right)^{2}}\left(dr_{x}\right)^{2}-\frac{r_{x}r_{y}}{\left(r_{x}^{2}+r_{y}^{2}\right)^{2}}\left(dr_{y}\right)^{2}+\frac{r_{y}^{2}-r_{x}^{2}}{\left(r_{x}^{2}+r_{y}^{2}\right)^{2}}dr_{x}dr_{y}=0,$
(90)

allowing us to write the final expression for $d\phi$ as

$d\phi=2\epsilon dt-\lambda\frac{\left(\beta\epsilon
r_{z}+2\right)\sin\phi}{\sqrt{r^{2}-r_{z}^{2}}}dW_{x}+\lambda\frac{\left(\beta\epsilon
r_{z}+2\right)\cos\phi}{\sqrt{r^{2}-r_{z}^{2}}}dW_{y}.$ (91)

###  A.3 Averaged system entropy production and boundary terms

We derive an analytical expression for $d\langle\langle\Delta
s_{sys}\rangle\rangle$, the noise- and coordinate-averaged, incremental system
entropy production, written as

$\displaystyle d\langle\langle\Delta s_{sys}\rangle\rangle$
$\displaystyle=\int_{-1}^{1}d\langle\Delta s_{sys}\rangle\,p(r_{z},t)dr_{z},$
(92)

where $d\langle\Delta s_{sys}\rangle$ is given by the noise averaged form of
Eq. (11). We proceed as follows:

$\displaystyle d\langle\langle\Delta s_{sys}\rangle\rangle$
$\displaystyle=\int_{-1}^{1}\left(-\frac{\partial\ln p(r_{z},t)}{\partial
t}dt-\frac{\partial\ln p(r_{z},t)}{\partial r_{z}}d\langle r_{z}\rangle-
D_{zz}\frac{\partial^{2}\ln p(r_{z},t)}{\partial
r_{z}^{2}}dt\right)p(r_{z},t)dr_{z}$
$\displaystyle=\int_{-1}^{1}\left(-\frac{\partial\ln p(r_{z},t)}{\partial
t}dt-A_{z}\frac{\partial\ln p(r_{z},t)}{\partial r_{z}}dt-
D_{zz}\frac{\partial^{2}\ln p(r_{z},t)}{\partial
r_{z}^{2}}dt\right)p(r_{z},t)dr_{z}.$ (93)

Integrating the last term by parts, and introducing the $r_{z}$ component of
the probability current,

$J_{z}=A_{z}p(r_{z},t)-\frac{\partial\left(D_{zz}p(r_{z},t)\right)}{\partial
r_{z}}\,,$ (94)

that appears in the Fokker-Planck equation (8) as

$\frac{\partial p(r_{z},t)}{\partial t}=-\frac{\partial J_{z}}{\partial
r_{z}}\,,$ (95)

we can write:

$\displaystyle d\langle\langle\Delta s_{sys}\rangle\rangle$
$\displaystyle=\int_{-1}^{1}\left(\left(-\frac{\partial\ln
p(r_{z},t)}{\partial t}dt-A_{z}\frac{\partial\ln p(r_{z},t)}{\partial
r_{z}}dt\right)p(r_{z},t)+\frac{\partial\left(D_{zz}p(r_{z},t)\right)}{\partial
r_{z}}\frac{\partial\ln p(r_{z},t)}{\partial r_{z}}dt\right)dr_{z}$
$\displaystyle\qquad\qquad-\left[D_{zz}\frac{\partial\ln p(r_{z},t)}{\partial
r_{z}}p(r_{z},t)\right]_{-1}^{1}dt$ (96)
$\displaystyle=\int_{-1}^{1}\left(-\frac{\partial\ln p(r_{z},t)}{\partial
t}p(r_{z},t)dt-J_{z}\frac{\partial\ln p(r_{z},t)}{\partial
r_{z}}dt\right)dr_{z}-\left[D_{zz}\frac{\partial p(r_{z},t)}{\partial
r_{z}}\right]_{-1}^{1}dt.$

Next we integrate the second term inside the integral by parts and replace the
$r_{z}$ derivative of the probability current with the left hand side of the
Fokker-Planck equation (95):

$\displaystyle d\langle\langle\Delta s_{sys}\rangle\rangle$
$\displaystyle=\int_{-1}^{1}\left(-\frac{\partial\ln p(r_{z},t)}{\partial
t}p(r_{z},t)dt+\ln p(r_{z},t)\,\frac{\partial J_{z}}{\partial
r_{z}}dt\right)dr_{z}-\left[J_{z}\ln
p(r_{z},t)\right]_{-1}^{1}dt-\left[D_{zz}\frac{\partial p(r_{z},t)}{\partial
r_{z}}\right]_{-1}^{1}dt$ $\displaystyle=-\int_{-1}^{1}\left(\frac{\partial\ln
p(r_{z},t)}{\partial t}p(r_{z},t)dt+\ln p(r_{z},t)\,\frac{\partial
p(r_{z},t)}{\partial t}dt\right)dr_{z}-\left[J_{z}\ln
p(r_{z},t)\right]_{-1}^{1}dt-\left[D_{zz}\frac{\partial p(r_{z},t)}{\partial
r_{z}}\right]_{-1}^{1}dt$
$\displaystyle=-\left(\frac{d}{dt}\int_{-1}^{1}p(r_{z},t)\ln
p(r_{z},t)dr_{z}\right)dt-\left[J_{z}\ln
p(r_{z},t)\right]_{-1}^{1}dt-\left[D_{zz}\frac{\partial p(r_{z},t)}{\partial
r_{z}}\right]_{-1}^{1}dt$ $\displaystyle=dS_{G}(t)-\left[J_{z}\ln
p(r_{z},t)\right]_{-1}^{1}dt-\left[D_{zz}\frac{\partial p(r_{z},t)}{\partial
r_{z}}\right]_{-1}^{1}dt,$ (97)

where $dS_{G}$ is the increment of the Gibbs entropy, see Eq. (13).

## References

  * [1] Peter W Shor.  Scheme for reducing decoherence in quantum computer memory.  Physical review A, 52(4):R2493, 1995.
  * [2] Ulrich Weiss.  Quantum dissipative systems, volume 13.  World scientific, 2012.
  * [3] Ken Sekimoto.  Langevin equation and thermodynamics.  Progress of Theoretical Physics Supplement, 130:17–27, 1998.
  * [4] Udo Seifert.  Entropy production along a stochastic trajectory and an integral fluctuation theorem.  Physical review letters, 95(4):040602, 2005.
  * [5] Richard E Spinney and Ian J Ford.  Entropy production in full phase space for continuous stochastic dynamics.  Physical Review E, 85(5):051113, 2012.
  * [6] Howard Carmichael.  An open systems approach to quantum optics: lectures presented at the Université Libre de Bruxelles, October 28 to November 4, 1991, volume 18.  Springer Science & Business Media, 2009.
  * [7] Heinz-Peter Breuer, Francesco Petruccione, et al.  The theory of open quantum systems.  Oxford University Press on Demand, 2002.
  * [8] Jürgen T Stockburger and Hermann Grabert.  Non-markovian quantum state diffusion.  Chemical Physics, 268(1-3):249–256, 2001.
  * [9] Jürgen T Stockburger and Hermann Grabert.  Exact c-number representation of non-markovian quantum dissipation.  Physical review letters, 88(17):170407, 2002.
  * [10] Jürgen T Stockburger.  Simulating spin-boson dynamics with stochastic liouville–von neumann equations.  Chemical physics, 296(2-3):159–169, 2004.
  * [11] GMG McCaul, CD Lorenz, and L Kantorovich.  Partition-free approach to open quantum systems in harmonic environments: An exact stochastic liouville equation.  Physical Review B, 95(12):125124, 2017.
  * [12] MA Lane, D Matos, IJ Ford, and L Kantorovich.  Exactly thermalised quantum dynamics of the spin-boson model coupled to a dissipative environment.  arXiv preprint arXiv:2002.07700, 2020.
  * [13] Richard Phillips Feynman and FL Vernon Jr.  The theory of a general quantum system interacting with a linear dissipative system.  Annals of physics, 24:118–173, 1963.
  * [14] Denis J Evans, Ezechiel Godert David Cohen, and Gary P Morriss.  Probability of second law violations in shearing steady states.  Physical review letters, 71(15):2401, 1993.
  * [15] Denis J Evans and Debra J Searles.  Steady states, invariant measures, and response theory.  Physical Review E, 52(6):5839, 1995.
  * [16] Denis J Evans and Debra J Searles.  The fluctuation theorem.  Advances in Physics, 51(7):1529–1585, 2002.
  * [17] DM Carberry, James Cowie Reid, GM Wang, Edith M Sevick, Debra J Searles, and Denis J Evans.  Fluctuations and irreversibility: An experimental demonstration of a second-law-like theorem using a colloidal particle held in an optical trap.  Physical review letters, 92(14):140601, 2004.
  * [18] Giovanni Gallavotti and Ezechiel Godert David Cohen.  Dynamical ensembles in nonequilibrium statistical mechanics.  Physical review letters, 74(14):2694, 1995.
  * [19] Jorge Kurchan.  Fluctuation theorem for stochastic dynamics.  Journal of Physics A: Mathematical and General, 31(16):3719, 1998\.
  * [20] Joel L Lebowitz and Herbert Spohn.  A gallavotti–cohen-type symmetry in the large deviation functional for stochastic dynamics.  Journal of Statistical Physics, 95(1):333–365, 1999.
  * [21] Christopher Jarzynski.  Nonequilibrium equality for free energy differences.  Physical Review Letters, 78(14):2690, 1997.
  * [22] Gavin E Crooks.  Nonequilibrium measurements of free energy differences for microscopically reversible markovian systems.  Journal of Statistical Physics, 90(5):1481–1487, 1998.
  * [23] Gavin E Crooks.  Entropy production fluctuation theorem and the nonequilibrium work relation for free energy differences.  Physical Review E, 60(3):2721, 1999.
  * [24] Rosemary J Harris and Gunther M Schütz.  Fluctuation theorems for stochastic dynamics.  Journal of Statistical Mechanics: Theory and Experiment, 2007(07):P07020, 2007.
  * [25] Ian J Ford, Zachary PL Laker, and Henry J Charlesworth.  Stochastic entropy production arising from nonstationary thermal transport.  Physical review E, 92(4):042108, 2015.
  * [26] Kurt Jacobs.  Quantum measurement theory and its applications.  Cambridge University Press, 2014.
  * [27] Kurt Jacobs and Daniel A Steck.  A straightforward introduction to continuous quantum measurement.  Contemporary Physics, 47(5):279–303, 2006.
  * [28] Gavin E Crooks.  Quantum operation time reversal.  Physical Review A, 77(3):034101, 2008.
  * [29] Jordan M Horowitz and Juan MR Parrondo.  Entropy production along nonequilibrium quantum jump trajectories.  New Journal of Physics, 15(8):085028, 2013.
  * [30] B Leggio, A Napoli, A Messina, and H-P Breuer.  Entropy production and information fluctuations along quantum trajectories.  Physical Review A, 88(4):042111, 2013.
  * [31] Cyril Elouard, David A Herrera-Martí, Maxime Clusel, and Alexia Auffeves.  The role of quantum measurement in stochastic thermodynamics.  npj Quantum Information, 3(1):1–10, 2017.
  * [32] Cyril Elouard, NK Bernardes, ARR Carvalho, MF Santos, and A Auffèves.  Probing quantum fluctuation theorems in engineered reservoirs.  New Journal of Physics, 19(10):103011, 2017.
  * [33] Juliette Monsel, Cyril Elouard, and Alexia Auffeves.  An autonomous quantum machine to measure the thermodynamic arrow of time.  npj Quantum Information, 4(1):1–9, 2018.
  * [34] Justin Dressel, Areeya Chantasri, Andrew N Jordan, and Alexander N Korotkov.  Arrow of time for continuous quantum measurement.  Physical review letters, 119(22):220507, 2017.
  * [35] Sreenath K Manikandan and Andrew N Jordan.  Time reversal symmetry of generalized quantum measurements with past and future boundary conditions.  Quantum studies: mathematics and foundations, 6(2):241–268, 2019\.
  * [36] Sreenath K Manikandan, Cyril Elouard, and Andrew N Jordan.  Fluctuation theorems for continuous quantum measurements and absolute irreversibility.  Physical Review A, 99(2):022117, 2019.
  * [37] Hannes Risken.  Fokker-planck equation.  In The Fokker-Planck Equation, pages 63–95. Springer, 1996.
  * [38] Daniel Matos, Matthew A Lane, Ian J Ford, and Lev Kantorovich.  Efficient choice of colored noise in the stochastic dynamics of open quantum systems.  Physical Review E, 102(6):062134, 2020.
  * [39] Nancy Makri and Dmitrii E Makarov.  Tensor propagator for iterative quantum time evolution of reduced density matrices. i. theory.  The Journal of chemical physics, 102(11):4600–4610, 1995.
  * [40] Yun-an Yan, Fan Yang, Yu Liu, and Jiushu Shao.  Hierarchical approach based on stochastic decoupling to dissipative systems.  Chemical physics letters, 395(4-6):216–221, 2004.
  * [41] D Suess, A Eisfeld, and WT Strunz.  Hierarchy of stochastic pure states for open quantum system dynamics.  Physical review letters, 113(15):150403, 2014.
  * [42] Peter P Orth, Adilet Imambekov, and Karyn Le Hur.  Nonperturbative stochastic method for driven spin-boson model.  Physical Review B, 87(1):014305, 2013.
  * [43] Sadao Nakajima.  On quantum theory of transport phenomena: Steady diffusion.  Progress of Theoretical Physics, 20(6):948–959, 1958.
  * [44] Robert Zwanzig.  Ensemble method in the theory of irreversibility.  The Journal of Chemical Physics, 33(5):1338–1341, 1960.
  * [45] Hazime Mori.  Transport, collective motion, and brownian motion.  Progress of theoretical physics, 33(3):423–455, 1965.
  * [46] Koichi Furutsu.  On the statistical theory of electromagnetic waves in a fluctuating medium.  United States Department of Commerce, 1964.
  * [47] Evgenii A Novikov.  Functionals and the random-force method in turbulence theory.  Sov. Phys. JETP, 20(5):1290–1294, 1965.
  * [48] Yun-An Yan and Jiushu Shao.  Stochastic description of quantum brownian dynamics.  Frontiers of Physics, 11(4):1–24, 2016.
  * [49] Amir O Caldeira and Anthony J Leggett.  Path integral approach to quantum brownian motion.  Physica A: Statistical mechanics and its Applications, 121(3):587–616, 1983.
  * [50] Heinz-Peter Breuer, Bernd Kappler, and Francesco Petruccione.  Stochastic wave-function method for non-markovian quantum master equations.  Physical Review A, 59(2):1633, 1999.
  * [51] Heinz-Peter Breuer.  Genuine quantum trajectories for non-markovian processes.  Physical Review A, 70(1):012106, 2004.
  * [52] Nicolas Gisin and Ian C Percival.  The quantum-state diffusion model applied to open systems.  Journal of Physics A: Mathematical and General, 25(21):5677, 1992\.
  * [53] Nicolas Gisin and Ian C Percival.  Quantum state diffusion, localization and quantum dispersion entropy.  Journal of Physics A: Mathematical and General, 26(9):2233, 1993\.
  * [54] Walter T Strunz.  Linear quantum state diffusion for non-markovian open quantum systems.  Physics Letters A, 224(1-2):25–30, 1996.
  * [55] Todd A Brun.  Continuous measurements, quantum trajectories, and decoherent histories.  Physical Review A, 61(4):042107, 2000.
  * [56] Ian Percival.  Quantum state diffusion.  Cambridge University Press, 1998.
  * [57] Claudia L. Clarke.  Irreversibility Measures in a Quantum Setting.  PhD thesis, UCL (University College London), 2021.
  * [58] Claudia L. Clarke and Ian J. Ford.  In preparation.  2022\.
  * [59] Robert G Sachs.  The physics of time reversal.  University of Chicago Press, 1987.