Stochastic entropy production for continuous measurements of an open quantum system

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–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 x (t) for 0 ≤ t ≤ τ as

$$\begin{eqnarray}&&{\rm{\Delta }}{s}_{{tot}}({\boldsymbol{x}})=\mathrm{ln}\displaystyle \frac{{p}^{F}({\boldsymbol{x}}(0),0){{ \mathcal P }}^{F}({\boldsymbol{x}}(\tau )\left|{\boldsymbol{x}}(0))\right.}{{p}^{F}({\boldsymbol{x}}(\tau ),\tau ){{ \mathcal P }}^{R}({{\boldsymbol{x}}}^{\dagger }(\tau )\left|{{\boldsymbol{x}}}^{\dagger }(0))\right.},\end{eqnarray} \tag{ 1 }$$

where x ^† is a set of coordinates evolving to form a time-reversed trajectory during the reverse protocol. Concretely, the reverse trajectory is defined as x ^†(t) = ε x (τ − t), where the elements of ε are +1 and −1 for even and odd variables, respectively. In this study we will only be dealing with even variables, such that the starting and finishing points of the reversed trajectory x ^†(0) and x ^†(τ) are x (τ) and 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], since no odd variables are involved here. p^F ( 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 x (0) to x (τ) under the forward protocol, and for the time reversed trajectory from x ^†(0) to x ^†(τ) under the reverse protocol, respectively.

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

$$\begin{eqnarray}&&{e}^{{\rm{\Delta }}{s}_{{tot}}}=\displaystyle \frac{{P}^{F}\left({\rm{\Delta }}{s}_{{tot}}\right)}{{P}^{R}\left(-{\rm{\Delta }}{s}_{{tot}}\right)},\end{eqnarray} \tag{ 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. Physical interpretation of stochastic entropy production for a quantum system

2.1.2. 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 each 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

$$\begin{eqnarray}&&\bar{\rho }(t+{dt})=\displaystyle \sum _{k}{M}_{k}\bar{\rho }(t){M}_{k}^{\dagger },\end{eqnarray} \tag{ 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 → 0.

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

$$\begin{eqnarray}&&\rho (t+{dt})=\rho (t)+d\rho =\displaystyle \frac{{M}_{k}\rho (t){M}_{k}^{\dagger }}{\mathrm{Tr}({M}_{k}\rho (t){M}_{k}^{\dagger })},\end{eqnarray} \tag{ 4 }$$

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

$$\begin{eqnarray}&&\displaystyle \sum _{k}{p}_{k}\displaystyle \frac{{M}_{k}\rho (t){M}_{k}^{\dagger }}{\mathrm{Tr}({M}_{k}\rho (t){M}_{k}^{\dagger })}=\displaystyle \sum _{k}{M}_{k}\rho (t){M}_{k}^{\dagger },\end{eqnarray} \tag{ 5 }$$

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

$$\begin{eqnarray}&&d\bar{\rho }=-i\left[{H}_{{sys}},\bar{\rho }\right]{dt}+\displaystyle \sum _{k}\left({L}_{k}\bar{\rho }{L}_{k}^{\dagger }-\displaystyle \frac{1}{2}\left\{{L}_{k}^{\dagger }{L}_{k},\bar{\rho }\right\}\right){dt},\end{eqnarray} \tag{ 6 }$$

with ℏ = 1, system Hamiltonian H_sys and 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–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–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 x involving a vector of Wiener increments d W, a vector of deterministic rate functions A and a matrix of noise coefficients B :

$$\begin{eqnarray}&&d{\boldsymbol{x}}={\boldsymbol{A}}({\boldsymbol{x}},t){dt}+{\boldsymbol{B}}({\boldsymbol{x}},t){\boldsymbol{d}}W.\end{eqnarray} \tag{ 7 }$$

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

$$\begin{eqnarray}&&\displaystyle \frac{\partial p({\boldsymbol{x}},t)}{\partial t}=\displaystyle \sum _{i}\displaystyle \frac{\partial }{\partial {x}_{i}}\left[-{A}_{i}({\boldsymbol{x}},t)p({\boldsymbol{x}},t)+\displaystyle \sum _{j}\displaystyle \frac{\partial }{\partial {x}_{j}}\left({D}_{{ij}}({\boldsymbol{x}},t)p({\boldsymbol{x}},t)\right)\right],\end{eqnarray} \tag{ 8 }$$

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

2.1.3. Stochastic entropy production

From the definition in equation (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 [5], first separating the deterministic part of equation (7) into two contributions:

$$\begin{eqnarray}&&d{\boldsymbol{x}}={{\boldsymbol{A}}}^{{irr}}({\boldsymbol{x}},t){dt}+{{\boldsymbol{A}}}^{{rev}}({\boldsymbol{x}},t){dt}+{\boldsymbol{B}}({\boldsymbol{x}},t){\boldsymbol{d}}W,\end{eqnarray} \tag{ 9 }$$

where A ^irr and A ^rev are the irreversible and reversible contributions to A , respectively. Specifically, if dx_i (the i-th component of d 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 D , the incremental environmental stochastic entropy production is [5]:

$$\begin{eqnarray}\begin{array}{rcl}d{\rm{\Delta }}{s}_{{env}} & = & \displaystyle \sum _{i=1}^{N}\left[\displaystyle \frac{{A}_{i}^{{ir}}}{{D}_{{ii}}}{{dx}}_{i}-\displaystyle \frac{{A}_{i}^{{rev}}{A}_{i}^{{irr}}}{{D}_{{ii}}}{dt}+\displaystyle \frac{\partial {A}_{i}^{{irr}}}{\partial {x}_{i}}{dt}-\displaystyle \frac{\partial {A}_{i}^{{rev}}}{\partial {x}_{i}}{dt}-\displaystyle \frac{1}{{D}_{{ii}}}\displaystyle \frac{\partial {D}_{{ii}}}{\partial {x}_{i}}{{dx}}_{i}\right.\\ & & \left.+\displaystyle \frac{{A}_{i}^{{rev}}-{A}_{i}^{{irr}}}{{D}_{{ii}}}\displaystyle \frac{\partial {D}_{{ii}}}{\partial {x}_{i}}{dt}-\displaystyle \frac{{\partial }^{2}{D}_{{ii}}}{\partial {x}_{i}^{2}}{dt}+\displaystyle \frac{1}{{D}_{{ii}}}{\left(\displaystyle \frac{\partial {D}_{{ii}}}{\partial {x}_{i}}\right)}^{2}{dt}\right],\end{array}\end{eqnarray} \tag{ 10 }$$

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 Δ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( x , t) of the system, obtained from the Fokker-Planck equation:

$$\begin{eqnarray}\begin{array}{rcl}d{\rm{\Delta }}{s}_{{sys}} & = & -d\left(\mathrm{ln}p({\boldsymbol{x}},t)\right)\\ & = & -\displaystyle \frac{\partial \mathrm{ln}p({\boldsymbol{x}},t)}{\partial t}{dt}-\displaystyle \sum _{i}\displaystyle \frac{\partial \mathrm{ln}p({\boldsymbol{x}},t)}{\partial {x}_{i}}{{dx}}_{i}-\displaystyle \sum _{i}{D}_{{ii}}\displaystyle \frac{{\partial }^{2}\mathrm{ln}p({\boldsymbol{x}},t)}{\partial {x}_{i}^{2}}{dt},\end{array}\end{eqnarray} \tag{ 11 }$$

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Δs_tot = dΔs_sys + dΔ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]:

$$\begin{eqnarray}&&\langle \langle {\rm{\Delta }}{s}_{{tot}}\rangle \rangle =\int d{\boldsymbol{x}}\,p({\boldsymbol{x}},{t}_{{init}})\mathrm{ln}\displaystyle \frac{p({\boldsymbol{x}},{t}_{{init}})}{p({\boldsymbol{x}},{t}_{{final}})}.\end{eqnarray} \tag{ 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

$$\begin{eqnarray}&&\langle \langle {\rm{\Delta }}{s}_{{sys}}\rangle \rangle ={\rm{\Delta }}{S}_{G}=-\int d{\boldsymbol{x}}\,p({\boldsymbol{x}},{t}_{{final}})\mathrm{ln}p({\boldsymbol{x}},{t}_{{final}})+\int d{\boldsymbol{x}}\,p({\boldsymbol{x}},{t}_{{init}})\mathrm{ln}p({\boldsymbol{x}},{t}_{{init}}),\end{eqnarray} \tag{ 13 }$$

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

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–12]. Other methods that do not rely on the Feynman-Vernon functional include projection operator methods such as the Nakajima-Zwanzig equation [43–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.2.

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 ξ_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

$$\begin{eqnarray}&&{H}_{{total}}(q,\{{\xi }_{{ik}}\},t)={H}_{{sys}}(q,t)+\displaystyle \sum _{k}\left[{H}_{{env};k}({\xi }_{{ik}})-\displaystyle \sum _{i}{\xi }_{{ik}}{f}_{k}(q,t)\right],\end{eqnarray} \tag{ 14 }$$

where H_env;k (ξ_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 ρ_q (t₀) and the environmental density matrix ρ_ξ (t₀), i.e. ρ₀ = ρ_q (t₀) ⨂ ρ_ξ (t₀), where t₀ 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 equation (14), it is possible to obtain a stochastic differential equation for the (unphysical) stochastic reduced density matrix ρ_s (t) of the system driven by a pair of complex coloured Gaussian noises, η_k (t), and ν_k (t), for each oscillator set k [8–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 ℏ = 1) becomes:

$$\begin{eqnarray}&&\displaystyle \frac{d{\rho }_{s}(t)}{{dt}}=-i\left[{H}_{{sys}},{\rho }_{s}(t)\right]+i\displaystyle \sum _{k}\left({\eta }_{k}(t)\left[{f}_{k},{\rho }_{s}(t)\right]+\displaystyle \frac{{\nu }_{k}(t)}{2}\left\{{f}_{k},{\rho }_{s}(t)\right\}\right),\end{eqnarray} \tag{ 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:

$$\begin{eqnarray}&&\langle {\eta }_{k}(t){\eta }_{k}(t^{\prime} )\rangle ={\int }_{0}^{\infty }\displaystyle \frac{d\omega }{\pi }{J}_{k}(\omega )\coth \left(\displaystyle \frac{1}{2}{\beta }_{k}\omega \right)\cos \left(\omega \left(t-t^{\prime} \right)\right)\equiv {K}_{k}^{\eta \eta }(t-t^{\prime} ),\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\langle {\eta }_{k}(t){\nu }_{k}(t^{\prime} )\rangle =-2i{\rm{\Theta }}(t-t^{\prime} ){\int }_{0}^{\infty }\displaystyle \frac{d\omega }{\pi }{J}_{k}(\omega )\sin \left(\omega \left(t-t^{\prime} \right)\right)\equiv {K}_{k}^{\eta \nu }(t-t^{\prime} ),\end{eqnarray} \tag{ 17 }$$

where J_k (ω) is the spectral density of the k^th oscillator set, β_k = 1/k_B T_k is its inverse temperature (with k_B being set to 1 hereafter), ${\rm{\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].

We need to average equation (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 equation (15) and denoting the physical density matrix averaged over the noises as $\bar{\rho }=\langle {\rho }_{s}\rangle$, we write:

$$\begin{eqnarray}&&\displaystyle \frac{d\bar{\rho }(t)}{{dt}}=-i\left[{H}_{{sys}},\bar{\rho }(t)\right]+i\displaystyle \sum _{k}\left(\left[{f}_{k},\langle {\eta }_{k}(t){\rho }_{s}(t)\rangle \right]+\displaystyle \frac{1}{2}\left\{{f}_{k},\langle {\nu }_{k}(t){\rho }_{s}(t)\rangle \right\}\right).\end{eqnarray} \tag{ 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 ζ_i with the correlation

$$\begin{eqnarray}&&\langle {\zeta }_{i}(t){\zeta }_{j}(t^{\prime} )\rangle ={F}_{{ij}}(t,t^{\prime} ),\end{eqnarray} \tag{ 19 }$$

the Furutsu-Novikov theorem states that

$$\begin{eqnarray}&&\langle {\zeta }_{i}(t)A[\zeta ]\rangle =\displaystyle \sum _{j}{\int }_{0}^{t}{dt}^{\prime} {F}_{{ij}}(t,t^{\prime} )\left\langle \displaystyle \frac{\delta A[\zeta ]}{\delta {\zeta }_{j}(t^{\prime} )}\right\rangle ,\end{eqnarray} \tag{ 20 }$$

where A[ζ] is a functional of the noises and δ A[ζ]/δ ζ_j is its functional derivative with respect to ζ_j . The averaged product of a noise and a noise-dependent functional may therefore be transformed into an integral. We use equations (20) in (18) to obtain the required explicit expressions for the averages of the noises multiplied by ρ_s :

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d\bar{\rho }(t)}{{dt}} & = & -i\left[{H}_{{sys}}(t),\bar{\rho }(t)\right]+i\displaystyle \sum _{k}\left[{f}_{k},{\int }_{0}^{t}{dt}^{\prime} {K}_{k}^{\eta \eta }(t-t^{\prime} )\left\langle \displaystyle \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 \displaystyle \frac{\delta {\rho }_{s}(t)}{\delta {\nu }_{k}(t^{\prime} )}\right\rangle \right]\\ & & +\,\displaystyle \frac{i}{2}\displaystyle \sum _{k}\left\{{f}_{k},{\int }_{0}^{t}{dt}^{\prime} {K}_{k}^{\eta \nu }(t^{\prime} -t)\left\langle \displaystyle \frac{\delta {\rho }_{s}(t)}{\delta {\eta }_{k}(t^{\prime} )}\right\rangle \right\}.\end{array}\end{eqnarray} \tag{ 21 }$$

From equation (17) we observe the presence of the Heaviside function ${\rm{\Theta }}(t-t^{\prime} )$, which vanishes for $t^{\prime} \gt t.$ The integrals in equation (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

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d\bar{\rho }(t)}{{dt}} & = & -i\left[{H}_{{sys}},\bar{\rho }(t)\right]-{\int }_{0}^{t}{dt}^{\prime} {K}^{\eta \eta }(t-t^{\prime} )\displaystyle \sum _{k}\left[{f}_{k},\left\langle {U}_{+}(t,t^{\prime} ){f}_{k}(t^{\prime} ){U}_{+}(t^{\prime} ,t){\rho }_{s}(t)\right.\right.\\ & & \left.\left.-{\rho }_{s}(t){U}_{-}(t,t^{\prime} ){f}_{k}(t^{\prime} ){U}_{-}(t^{\prime} ,t)\right\rangle \right]\\ & & +\displaystyle \frac{1}{2}{\int }_{0}^{t}{dt}^{\prime} {K}_{k}^{\eta \nu }(t-t^{\prime} )\displaystyle \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],\end{array}\end{eqnarray} \tag{ 22 }$$

where U₊ and U₋ are forward and reverse propagators, respectively (see appendix A.1). Equation (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 equation (6). However, this equation still contains the stochastic density matrix ρ_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.

Equation (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 ≡ 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, β ω ≪ 1, and we will also adopt the same Ohmic spectral density J(ω) = α ω for all oscillator sets, where α is a proportionality constant. In this case we can use the first order expansion of $\coth \left(\tfrac{\beta \omega }{2}\right)\approx \tfrac{2}{\beta \omega }$, following [48], and obtain the following results for the correlation functions of equations (16) and (17):

$$\begin{eqnarray}&&{K}^{\eta \eta }(t-t^{\prime} )\approx {\int }_{0}^{\infty }\displaystyle \frac{d\omega }{\pi }\alpha \omega \displaystyle \frac{2}{\beta \omega }\cos \left(\omega (t-t^{\prime} )\right)=\displaystyle \frac{2\alpha }{\beta }\delta (t-t^{\prime} ),\end{eqnarray} \tag{ 23 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}{K}^{\eta \nu }(t-t^{\prime} ) & = & 2i\alpha {\rm{\Theta }}(t-t^{\prime} )\displaystyle \frac{\partial }{\partial t}\left[\displaystyle \frac{1}{\pi }{\int }_{0}^{\infty }d\omega \cos \left(\omega (t-t^{\prime} )\right)\right]\\ & = & -2i\alpha {\rm{\Theta }}(t-t^{\prime} )\displaystyle \frac{\partial }{\partial t^{\prime} }\left[\displaystyle \frac{1}{\pi }{\int }_{0}^{\infty }d\omega \cos \left(\omega (t-t^{\prime} )\right)\right]\\ & = & -2i\alpha {\rm{\Theta }}(t-t^{\prime} )\displaystyle \frac{\partial }{\partial t^{\prime} }\delta (t-t^{\prime} ).\end{array}\end{eqnarray} \tag{ 24 }$$

Note that the correlation functions are the same for all the oscillator sets, hence we have dropped the index k here. Equations (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 equations (23) and (24) into equation (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]:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d\bar{\rho }(t)}{{dt}} & = & -i\left[{H}_{{sys}},\bar{\rho }(t)\right]-\displaystyle \sum _{k}\left[{f}_{k},\displaystyle \frac{\alpha }{\beta }\left[{f}_{k},\bar{\rho }(t)\right]-\displaystyle \frac{\alpha }{2}\left\{\left[{H}_{{sys}},{f}_{k}\right],\bar{\rho }(t)\right\}+\displaystyle \frac{i\alpha }{2}\left\{\displaystyle \frac{\partial {f}_{k}}{\partial t},\bar{\rho }(t)\right\}\right]\\ & = & -i\left[{H}_{{sys}},\bar{\rho }(t)\right]+\displaystyle \sum _{k}\left(-\displaystyle \frac{\alpha }{\beta }\left[{f}_{k},\left[{f}_{k},\bar{\rho }(t)\right]\right]+\displaystyle \frac{\alpha }{2}\left[{f}_{k},\left\{\left[{H}_{{sys}},{f}_{k}\right],\bar{\rho }(t)\right\}\right]-\displaystyle \frac{i\alpha }{2}\left[{f}_{k},\left\{\displaystyle \frac{\partial {f}_{k}}{\partial t},\bar{\rho }(t)\right\}\right]\right),\end{array}\end{eqnarray} \tag{ 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.1.1. Two-level system and choice of environments

We consider a two-level bosonic system with Hamiltonian

$$\begin{eqnarray}&&{H}_{{sys}}=\epsilon {\sigma }_{z},\end{eqnarray} \tag{ 26 }$$

where σ_z is the usual Pauli matrix and has the dimension of energy. A more complicated Hamiltonian could be used but would make our core message less clear, namely that in order to perform a measurement of an observable, the coupling of the system to the environment by way of that observable needs to be elevated. We can therefore proceed with a simple framework of environmental couplings through the Pauli matrices, as will be seen shortly. 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 equation (25), with time independent f_k , to vanish:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d\bar{\rho }}{{dt}} & \propto & -i\left[{H}_{{sys}},{\mathbb{I}}-\beta {H}_{{sys}}\right]+\displaystyle \frac{\alpha }{\beta }\displaystyle \sum _{k}\left(-\left[{f}_{k},\left[{f}_{k},{\mathbb{I}}-\beta {H}_{{sys}}\right]\right]+\displaystyle \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}\displaystyle \sum _{k}\left[{f}_{k},\left\{\left[{H}_{{sys}},{f}_{k}\right],{H}_{{sys}}\right\}\right].\end{array}\end{eqnarray} \tag{ 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 equation (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 σ_x or σ_y as well as the choice in equation (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 σ_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, equation (26), has two eigenvectors ∣ ± 〉 with eigenvalues ±, respectively, and the environmental coupling to σ_z drives the reduced density matrix stochastically towards the form ∣ + 〉〈 + ∣ or ∣ − 〉〈 − ∣. 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 σ_z .

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

$$\begin{eqnarray}&&{f}_{x}={\gamma }_{0}{\sigma }_{x},\,{f}_{y}={\gamma }_{0}{\sigma }_{y}\,\,\mathrm{and}\,\,{f}_{z}=\gamma (t){\sigma }_{z},\end{eqnarray} \tag{ 28 }$$

where for simplicity we replace labels k by the Cartesian indices x, y and z. γ₀ is the coupling coefficient between the system and the environment, and δ γ(t) = γ(t) − γ₀ is the coupling coefficient between the system and the device measuring H_sys . If γ(t) is increased from an initial value of γ₀, 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 γ(t) returns to γ₀ 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 equation (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 ∂f_k /∂t is proportional to σ_z .

3.1.2. Unravelled stochastic equation

To obtain a form of equation (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 [55, 56], though that approach will not be taken here. Instead we construct Lindblad operators that match the dynamics of equation (25) in accordance with our high temperature approximation β H_sys ≪ 1.

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

$$\begin{eqnarray}&&{L}_{k}=\lambda \left({f}_{k}-\displaystyle \frac{1}{4}\left[\beta {H}_{{sys}},{f}_{k}\right]\right),\end{eqnarray} \tag{ 29 }$$

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

$$\begin{eqnarray}\begin{array}{rcl}d\bar{\rho } & = & -i\left[{H}_{{sys}}-\displaystyle \frac{i\alpha }{4}\displaystyle \sum _{k}\left[{H}_{{sys}},{f}_{k}^{2}\right],\bar{\rho }\right]{dt}+\displaystyle \sum _{k}\left(-\displaystyle \frac{{\lambda }^{2}}{2}\left[{f}_{k},\left[{f}_{k},\bar{\rho }\right]\right]+\displaystyle \frac{{\lambda }^{2}}{4}\beta \left[{f}_{k},\left\{\left[{H}_{{sys}},{f}_{k}\right],\bar{\rho }\right\}\right]\right){dt}\\ & & -\,\displaystyle \frac{{\lambda }^{2}}{16}\displaystyle \sum _{k}{\beta }^{2}\left(\left[{H}_{{sys}},{f}_{k}\right]\bar{\rho }\left[{H}_{{sys}},{f}_{k}\right]-\displaystyle \frac{1}{2}\left\{{\left[{H}_{{sys}},{f}_{k}\right]}^{2},\bar{\rho }\right\}\right){dt}.\end{array}\end{eqnarray} \tag{ 30 }$$

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 β H_sys smaller in magnitude than the penultimate term, and therefore can be neglected. With these considerations we recover the Markovian averaged SLN equation (25) confirming the suitability of the Lindblad operators in equation (29).

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

$$\begin{eqnarray}&&{M}_{\pm k}\propto \left({\mathbb{I}}-{{iH}}_{{sys}}{dt}-\displaystyle \frac{1}{2}{L}_{k}^{\dagger }{L}_{k}{dt}\pm {L}_{k}\sqrt{{dt}}\right),\end{eqnarray} \tag{ 31 }$$

consistent with the development in section 2.1.2, and substituting equation (31) into equation (4) we obtain (see [26, 58]):

$$\begin{eqnarray}&&d\rho =-i\left[{H}_{{sys}},\rho \right]{dt}+\displaystyle \sum _{k}\left[\left({L}_{k}\rho {L}_{k}^{\dagger }-\displaystyle \frac{1}{2}\left\{{L}_{k}^{\dagger }{L}_{k},\rho \right\}\right){dt}+\left(\rho {L}_{k}^{\dagger }+{L}_{k}\rho -\mathrm{Tr}\left[\rho \left({L}_{k}+{L}_{k}^{\dagger }\right)\right]\rho \right){{dW}}_{k}\right],\end{eqnarray} \tag{ 32 }$$

where we use ρ 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 equation (6) and hence $\bar{\rho }$ corresponds to the noise averaged reduced density matrix ρ. Equation (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 equation (28) (k = x, y, z), and with γ₀ = 1, we obtain Lindblad operators:

$$\begin{eqnarray}&&{L}_{x}=\lambda \left({\sigma }_{x}-i\displaystyle \frac{\beta \epsilon }{2}{\sigma }_{y}\right),\quad {L}_{y}=\lambda \left({\sigma }_{y}+i\displaystyle \frac{\beta \epsilon }{2}{\sigma }_{x}\right),\quad {L}_{z}=\lambda \gamma {\sigma }_{z}.\end{eqnarray} \tag{ 33 }$$

Instead of simulating the dynamics of elements of the density matrix directly, we will evolve the coherence or Bloch vector defined by the three components ${r}_{i}=\mathrm{Tr}({\sigma }_{i}\rho )$. This is consistent with a representation of the reduced density matrix in the form

$$\begin{eqnarray}&&\rho =\displaystyle \frac{1}{2}\left({\mathbb{I}}+{\boldsymbol{r}}\cdot {\boldsymbol{\sigma }}\right).\end{eqnarray} \tag{ 34 }$$

We obtain from equation (32) the following equations:

$$\begin{eqnarray}&&{{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}\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{{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}\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&{{dr}}_{z}=-4{\lambda }^{2}\left(\beta \epsilon +{r}_{z}\right){dt}-\lambda {r}_{x}\left(\beta \epsilon +2{r}_{z}\right){{dW}}_{x}-\lambda {r}_{y}\left(\beta \epsilon +2{r}_{z}\right){{dW}}_{y}+2\gamma \lambda \left(1-{r}_{z}^{2}\right){{dW}}_{z}.\end{eqnarray} \tag{ 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 〈r_x 〉 → 0, 〈r_y 〉 → 0 and 〈r_z 〉 → − β at long times and hence is not disturbed by the evolution of γ 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², r_z , ϕ), where

$$\begin{eqnarray}&&{r}^{2}={r}_{x}^{2}+{r}_{y}^{2}+{r}_{z}^{2}\,\,\mathrm{and}\,\,\phi =\arctan \left(\displaystyle \frac{{r}_{y}}{{r}_{x}}\right).\end{eqnarray} \tag{ 38 }$$

Using Itô's Lemma, see appendix A.2, we can obtain:

$$\begin{eqnarray}\begin{array}{rcl}{{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{{r}^{2}-{r}_{z}^{2}}\cos \phi \left(1-{r}^{2}\right){{dW}}_{x}\\ & & +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},\end{array}\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}\begin{array}{rcl}{{dr}}_{z} & = & -4{\lambda }^{2}\left(\beta \epsilon +{r}_{z}\right){dt}-\lambda \left(\beta \epsilon +2{r}_{z}\right)\sqrt{{r}^{2}-{r}_{z}^{2}}\cos \phi \,{{dW}}_{x}\\ & & -\lambda \left(\beta \epsilon +2{r}_{z}\right)\sqrt{{r}^{2}-{r}_{z}^{2}}\sin \phi \,{{dW}}_{y}+2\gamma \lambda \left(1-{r}_{z}^{2}\right){{dW}}_{z},\end{array}\end{eqnarray} \tag{ 40 }$$

and

$$\begin{eqnarray}&&d\phi =2\epsilon {dt}-\lambda \displaystyle \frac{\left(\beta \epsilon {r}_{z}+2\right)\sin \phi }{\sqrt{{r}^{2}-{r}_{z}^{2}}}{{dW}}_{x}+\lambda \displaystyle \frac{\left(\beta \epsilon {r}_{z}+2\right)\cos \phi }{\sqrt{{r}^{2}-{r}_{z}^{2}}}{{dW}}_{y}.\end{eqnarray} \tag{ 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 equation (39) that r = 1 is a fixed point of the dynamics. Defining a pure state by the condition $\mathrm{Tr}({\rho }^{2})=1$, it can be shown that r = 1 corresponds to a pure state by inserting ρ in its coherence vector representation 34. In order to reduce the computational difficulty of solving the equations, we set r = 1 and consider

$$\begin{eqnarray}\begin{array}{rcl}{{dr}}_{z} & = & -4{\lambda }^{2}\left(\beta \epsilon +{r}_{z}\right){dt}-\lambda \left(\beta \epsilon +2{r}_{z}\right)\sqrt{1-{r}_{z}^{2}}\cos \phi \,{{dW}}_{x}\\ & & -\lambda \left(\beta \epsilon +2{r}_{z}\right)\sqrt{1-{r}_{z}^{2}}\sin \phi \,{{dW}}_{y}+2\gamma \lambda \left(1-{r}_{z}^{2}\right){{dW}}_{z}\end{array}\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&d\phi =2\epsilon {dt}-\lambda \displaystyle \frac{\left(\beta \epsilon {r}_{z}+2\right)\sin \phi }{\sqrt{1-{r}_{z}^{2}}}{{dW}}_{x}+\lambda \displaystyle \frac{\left(\beta \epsilon {r}_{z}+2\right)\cos \phi }{\sqrt{1-{r}_{z}^{2}}}{{dW}}_{y}.\end{eqnarray} \tag{ 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 ϕ. These equations form the basis of our simulations with corresponding results to be presented in section 4.

The definition of the entropy production in equation (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 A and matrices B and ${\boldsymbol{D}}=\tfrac{1}{2}{\boldsymbol{B}}{{\boldsymbol{B}}}^{T}$, obtained by comparing their definitions in equation (7) with the equations of motion (42) and (43). We get

$$\begin{eqnarray}&&{\boldsymbol{A}}=\left(\begin{array}{c}{A}_{z}\\ {A}_{\phi }\end{array}\right)=\left(\begin{array}{c}-4{\lambda }^{2}(\beta \epsilon +{r}_{z})\\ 2\epsilon \end{array}\right),\end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray}{\boldsymbol{B}}=\lambda \left(\begin{array}{ccc}-(\beta \epsilon +2{r}_{z})\sqrt{1-{r}_{z}^{2}}\cos \phi & -(\beta \epsilon +2{r}_{z})\sqrt{1-{r}_{z}^{2}}\sin \phi & 2\gamma (1-{r}_{z}^{2})\\ -(\beta \epsilon {r}_{z}+2)\displaystyle \frac{\sin \phi }{\sqrt{1-{r}_{z}^{2}}} & (\beta \epsilon {r}_{z}+2)\displaystyle \frac{\cos \phi }{\sqrt{1-{r}_{z}^{2}}} & 0\end{array}\right),\end{eqnarray} \tag{ 45 }$$

and hence

$$\begin{eqnarray}{\boldsymbol{D}}=\left(\begin{array}{cc}{D}_{{zz}} & 0\\ 0 & {D}_{\phi \phi }\end{array}\right),\end{eqnarray} \tag{ 46 }$$

with

$$\begin{eqnarray}&&{D}_{{zz}}=2{\lambda }^{2}\left(1-{r}_{z}^{2}\right)\left({\left(\displaystyle \frac{\beta \epsilon }{2}\right)}^{2}+\beta \epsilon {r}_{z}+(1-{\gamma }^{2}){r}_{z}^{2}+{\gamma }^{2}\right),\end{eqnarray} \tag{ 47 }$$

and

$$\begin{eqnarray}&&{D}_{\phi \phi }=2{\lambda }^{2}\displaystyle \frac{{\left(\tfrac{\beta \epsilon {r}_{z}}{2}\right)}^{2}+\beta \epsilon {r}_{z}+1}{1-{r}_{z}^{2}}.\end{eqnarray} \tag{ 48 }$$

There exists no ϕ dependence in either A or D , and the latter matrix is indeed diagonal, as anticipated. With these expressions, the Fokker-Planck equation for the pdf p(r_z , ϕ, t) becomes:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial }{\partial t}p({r}_{z},\phi ,t) & = & \displaystyle \frac{\partial }{\partial {r}_{z}}\left[-{A}_{z}p({r}_{z},\phi ,t)+\displaystyle \frac{\partial }{\partial {r}_{z}}\left({D}_{{zz}}p({r}_{z},\phi ,t)\right)\right]\\ & & +\displaystyle \frac{\partial }{\partial \phi }\left[-{A}_{\phi }p({r}_{z},\phi ,t)+\displaystyle \frac{\partial }{\partial \phi }\left({D}_{\phi \phi }p({r}_{z},\phi ,t)\right)\right]\\ & = & \displaystyle \frac{\partial }{\partial {r}_{z}}\left[\displaystyle \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.\\ & & \left.+2{\lambda }^{2}\left(1-{r}_{z}^{2}\right)\left({\left(\displaystyle \frac{\beta \epsilon }{2}\right)}^{2}+\beta \epsilon {r}_{z}+(1-{\gamma }^{2}){r}_{z}^{2}+{\gamma }^{2}\right)\displaystyle \frac{\partial }{\partial {r}_{z}}p({r}_{z},\phi ,t)\right]\\ & & +\displaystyle \frac{\partial }{\partial \phi }\left[-2\epsilon p({r}_{z},\phi ,t)+2{\lambda }^{2}\displaystyle \frac{{\left(\tfrac{\beta \epsilon {r}_{z}}{2}\right)}^{2}+\beta \epsilon {r}_{z}+1}{1-{r}_{z}^{2}}\displaystyle \frac{\partial }{\partial \phi }p({r}_{z},\phi ,t)\right].\end{array}\end{eqnarray} \tag{ 49 }$$

The Fokker-Planck equation can be simplified further if we assume the initial pdf to be independent of ϕ. There is no explicit ϕ dependence in equation (49), so an initial pdf that is ϕ independent will always remain so. This leads us to the simplification p(r_z , ϕ, t) → p(r_z , t), and a more succinct equation:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial }{\partial t}p({r}_{z},t) & = & \displaystyle \frac{\partial }{\partial {r}_{z}}\left[\displaystyle \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.\\ & & \left.+2{\lambda }^{2}\left(1-{r}_{z}^{2}\right)\left({\left(\displaystyle \frac{\beta \epsilon }{2}\right)}^{2}+\beta \epsilon {r}_{z}+(1-{\gamma }^{2}){r}_{z}^{2}+{\gamma }^{2}\right)\displaystyle \frac{\partial }{\partial {r}_{z}}p({r}_{z},t)\right].\end{array}\end{eqnarray} \tag{ 50 }$$

The stationary solution p_st (r_z ) of this Fokker-Planck equation, corresponding to setting ∂p(r_z , t)/∂t = 0, can be obtained analytically using tools such as Mathematica. For γ ≠ 1, we have:

$$\begin{eqnarray}\begin{array}{rcl}{p}_{{st}}{\left({r}_{z}\right)}_{\gamma \ne 1} & = & \displaystyle \frac{1}{{N}_{\gamma }}{\left(1-{r}_{z}\right)}^{a}{\left(1+{r}_{z}\right)}^{b}{\left({\beta }^{2}{\epsilon }^{2}+4\beta \epsilon {r}_{z}+4{\gamma }^{2}+4(1-{\gamma }^{2})\right)}^{c}\\ & & \times \,\exp \left[\displaystyle \frac{8\beta \epsilon \left({\beta }^{2}{\epsilon }^{2}(1+2{\gamma }^{2})-4\right)\mathrm{arctanh}\left[\tfrac{-\beta \epsilon -2{r}_{z}(1-{\gamma }^{2})}{\gamma \sqrt{4{\gamma }^{2}-4+{\beta }^{2}{\epsilon }^{2}}}\right]}{{\left(4-{\beta }^{2}{\epsilon }^{2}\right)}^{2}\gamma \sqrt{4{\gamma }^{2}-4+{\beta }^{2}{\epsilon }^{2}}}\right],\end{array}\end{eqnarray} \tag{ 51 }$$

while for γ = 1 we have

$$\begin{eqnarray}&&{p}_{{st}}{\left({r}_{z}\right)}_{\gamma =1}=\displaystyle \frac{1}{{N}_{1}}{\left(1-{r}_{z}\right)}^{a}{\left(1+{r}_{z}\right)}^{b}{\left({\beta }^{2}{\epsilon }^{2}+4\beta \epsilon {r}_{z}+4\right)}^{d},\end{eqnarray} \tag{ 52 }$$

where N_γ is a γ-dependent normalisation factor, and the exponents are given by

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

where we note that c ≠ d. The stationary pdfs in equations (51) and (52) are shown in figure 2 for a number of values of γ.

Zoom In Zoom Out Reset image size

We can see that the coupling strength γ has a very significant effect on the stationary probability density of the system coordinate r_z . Increasing γ concentrates the pdf closer to the boundaries and in an asymmetric manner due to the system Hamiltonian bias . By evolving γ 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 = ± 1, or equivalently at system energies $\mathrm{Tr}\left({H}_{{sys}}\rho \right)=\epsilon {r}_{z}$ equal to ±.

Using its definition in equation (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 noticed in practice that combining the average stochastic system entropy production with the average stochastic environmental entropy production equation (10), obtained numerically for time independent environmental coupling, does not yield a vanishing average total stochastic 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

$$\begin{eqnarray}&&d\langle \langle {\rm{\Delta }}{s}_{{sys}}\rangle \rangle ={\int }_{-1}^{1}d\langle {\rm{\Delta }}{s}_{{sys}}\rangle \,p({r}_{z},t){{dr}}_{z},\end{eqnarray} \tag{ 54 }$$

where d〈Δs_sys 〉 is the noise averaged increment in system entropy production, which is a function of r_z and t. (Hereafter, we occasionally drop 'stochastic' for simplicity of terminology). Integrating by parts twice and retaining the boundary terms, see appendix A.3 for details, we obtain:

$$\begin{eqnarray}&&d\langle \langle {\rm{\Delta }}{s}_{{sys}}\rangle \rangle ={{dS}}_{G}-{\left[{J}_{z}\mathrm{ln}p({r}_{z},t)\right]}_{-1}^{1}{dt}-{\left[{D}_{{zz}}\displaystyle \frac{\partial p({r}_{z},t)}{\partial {r}_{z}}\right]}_{-1}^{1}{dt},\end{eqnarray} \tag{ 55 }$$

where dS_G is the increment in the Gibbs entropy of the system, shown in equation (13), and J_z is the r_z -component of the probability current defined in appendix A.3.

In equation (13) we noted that d〈〈Δs_sys 〉〉 = dS_G held only under certain conditions, now shown to be where the boundary terms in equation (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, eliminating the boundary terms For our system, although D_zz = 0 at r_z = ± 1, it can be shown that the gradient ∂p_st (r_z , t)/∂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〈〈Δs_sys 〉〉 = 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 equations (52) and (51), together with the diffusion coefficient D_zz given in equation (46), it may be shown that the boundary correction involving J_z in equation (55) vanishes at equilibrium, and so does dS_G . This allows us to write the mean stationary system entropy production increment for general γ as

$$\begin{eqnarray}&&d\langle \langle {\rm{\Delta }}{s}_{{sys}}\rangle {\rangle }_{{st}}=-{\left[{D}_{{zz}}\displaystyle \frac{\partial {p}_{{st}}{\left({r}_{z}\right)}_{\gamma }}{\partial {r}_{z}}\right]}_{-1}^{1}{dt},\end{eqnarray} \tag{ 56 }$$

which does not vanish. In fact, for γ = 1 and to lowest order in β , it is possible to obtain a simple analytical expression:

$$\begin{eqnarray}&&d\langle \langle {\rm{\Delta }}{s}_{{sys}}\rangle {\rangle }_{{st}}^{\gamma =1}=-{\lambda }^{2}{\left(\beta \epsilon \right)}^{2}{dt}.\end{eqnarray} \tag{ 57 }$$

This shows that the exact boundary correction can be obtained analytically for γ = 1, though it becomes harder for γ ≠ 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 β before calculating the boundary correction.

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

$$\begin{eqnarray}&&{\rm{\Theta }}{\sigma }_{z}{{\rm{\Theta }}}^{-1}={\sigma }_{z}\qquad {\rm{\Theta }}{\sigma }_{y}{{\rm{\Theta }}}^{-1}=-{\sigma }_{y}\qquad {\rm{\Theta }}{\sigma }_{x}{{\rm{\Theta }}}^{-1}={\sigma }_{x}.\end{eqnarray} \tag{ 58 }$$

Thus r_x and r_z are even under time reversal, and r_y is odd, and hence ϕ is also odd, and we can separate the coefficients of the dt terms in equations (42) and (43) into their A ^irr and A ^rev components. We can write

$$\begin{eqnarray}&&{{\boldsymbol{A}}}^{{irr}}=\left(\begin{array}{c}-4{\lambda }^{2}\left(\beta \epsilon +{r}_{z}\right)\\ 0\end{array}\right),\quad {{\boldsymbol{A}}}^{{rev}}=\left(\begin{array}{c}0\\ 2\epsilon \end{array}\right),\end{eqnarray} \tag{ 59 }$$

and obtain an expression for the environmental entropy production using equation (10), since we already have the form of D from equation (46).

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 γ = 1 and the initial pdf defined as in figure 2; for t > 0 perform simulations of equations (42) and (43) using γ = 2;
• Disconnecting a measuring device (protocol $\overline{{\rm{M}}}$): begin at t = 0 with the system thermalised for γ = 2 and the initial pdf defined as in figure 2; then, for t > 0 perform simulations of equations (42) and (43) using γ = 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 γ = 1 or γ = 2, while the value of ϕ is randomly sampled from a uniform pdf in the range [0, 2π). Whenever we specify a particular value of γ, 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 β = 0.1, the Hamiltonian parameter = 1, the proportionality constant in the spectral density α = 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 pdf going from a stationary state at a particular value of γ to another stationary state defined by a different γ. The protocols will drive the system between two stationary states obtained from the Fokker-Planck equation at different γ 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.  
  Obtain the solution p(r_z , t) of the Fokker-Planck equation (50) using the appropriate stationary pdf to select the initial condition (at t = 0) for the trajectories; in both cases the boundary conditions at r_z = ± 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.  
  Start a loop over the noises (independently for each protocol):

  • (a)  
    generate noises dW_x , dW_y and dW_z for each time step dt;
  • (b)  
    run the stochastic dynamics for r_z (t) and ϕ(t) for both protocols using equations (42) and (43) based on initial values sampled from the corresponding γ-dependent stationary pdf;
  • (c)  
    for each time increment dt and using obtained values of r_z (t) and ϕ(t), calculate the incremental contributions to the environmental, dΔs_env , and system, $d{\rm{\Delta }}{s}_{{sys}}=-d\left(\mathrm{ln}\,p({r}_{z},t)\right)$, stochastic entropy productions via respectively, equation (10) and the first line of equation (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;
  • (d)  
    the incremental total stochastic entropy production, at each time step dt, is calculated as a sum, dΔs_tot = dΔs_env + dΔs_sys ; accumulate these contributions over an entire trajectory;
  • (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.  
  Once the averaged values of the entropy production increments are obtained from a large ensemble of trajectories, the boundary term from equation (55) is added to the ensemble average $d\left\langle \left\langle {\rm{\Delta }}{s}_{{tot}}\right\rangle \right\rangle$ as

  $$\begin{eqnarray}&&d\langle \langle {\rm{\Delta }}{s}_{{tot}}\rangle \rangle =d\langle \langle {\rm{\Delta }}{s}_{{env}}\rangle \rangle +d\langle \langle {\rm{\Delta }}{s}_{{sys}}\rangle \rangle -{\left[{D}_{{zz}}\displaystyle \frac{\partial {p}_{{st}}({r}_{z})}{\partial {r}_{z}}\right]}_{-1}^{1}.\end{eqnarray} \tag{ 60 }$$

  for the specific value of γ. This yields the increment of the mean total stochastic entropy production d〈〈Δs_tot 〉〉 for the given time step dt, separately for the two protocols.

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 equation (32) we have assumed β H_sys ≪ 1, so we can write ${\bar{\rho }}_{{eq}}\approx \tfrac{1}{2}({\mathbb{I}}-\beta {H}_{{sys}})$ and hence obtain the average of r_z as ${\bar{r}}_{z}=\mathrm{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 value of r_z does indeed remain constant throughout the simulation. Furthermore, the dynamics of ϕ maintains a pdf that is constant in ϕ which matches the assumptions used to derive equation (50). While the averaged system behaviour remains unchanged, the stochastic trajectories differ significantly, as 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 = ± 1, jumping between the two. Once the system reaches stationarity (equilibrium), the ensemble of trajectories adopts the stationary density ${p}_{{st}}{\left({r}_{z}\right)}_{\gamma }$ associated with the corresponding value of γ. One would expect this to match the density obtained 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.

Zoom In Zoom Out Reset image size

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 explicitly, and to recognise that all physical changes in the system are stochastic in nature. We 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 γ = 1 explores the Bloch sphere fairly broadly as seen in figure 2, but for γ = 2, the trajectory tends to dwell near the boundaries at r_z = ± 1.

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

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

The entropy production for protocol $\overline{{\rm{M}}}$ with γ = 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 fact, the mean entropy production in this case is higher than for the connection process. Note that entropy production tends to a constant asymptotically even for an individual trajectory. Both figures 5 and 6 display excellent 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.

Zoom In Zoom Out Reset image size

As a check on the accuracy of the entropy production results, we verify that the detailed fluctuation theorem equation (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. Figure 7 displays the asymptotic pdfs of total stochastic entropy production for protocols M and $\overline{{\rm{M}}}$. These distributions should satisfy the detailed fluctuation relation. We compare ${e}^{{\rm{\Delta }}{s}_{{tot}}}$ with the ratio of numerical entropy pdfs ${P}^{M}\left({\rm{\Delta }}{s}_{{tot}}\right)/{P}^{\overline{M}}\left(-{\rm{\Delta }}{s}_{{tot}}\right)$ and indeed find a very good adherence excluding deviations at the extremes of the range due to insufficient sampling.

Zoom In Zoom Out Reset image size

Here we shall consider the calculation of the functional derivatives of ρ_s with respect to the noises needed for the derivation of equation (22). The formal solution of the SLN equation is given by

$$\begin{eqnarray}&&{\rho }_{s}(t)={U}_{+}(t,0){\rho }_{s}(0){U}_{-}(0,t),\end{eqnarray} \tag{ 61 }$$

where U₊ and U₋ are the appropriate forward and backward propagators defined as

$$\begin{eqnarray}&&{U}_{+}(t,0)={\widehat{{ \mathcal T }}}_{+}\exp \left[-i{\int }_{0}^{t}{dt}^{\prime} {H}_{+}(t^{\prime} )\right]\end{eqnarray} \tag{ 62 }$$

$$\begin{eqnarray}&&{U}_{-}(0,t)={\widehat{{ \mathcal T }}}_{-}\exp \left[i{\int }_{0}^{t}{dt}^{\prime} {H}_{-}(t^{\prime} )\right],\end{eqnarray} \tag{ 63 }$$

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

$$\begin{eqnarray}&&{H}_{\pm }(t)={H}_{{sys}}(t)-\displaystyle \sum _{k}\left({\eta }_{k}(t)\pm \displaystyle \frac{{\nu }_{k}(t)}{2}\right){f}_{k}\end{eqnarray} \tag{ 64 }$$

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

$$\begin{eqnarray}&&\displaystyle \frac{\delta {\rho }_{s}(t)}{\delta \zeta (t^{\prime} )}=\displaystyle \frac{\delta {U}_{+}(t,0)}{\delta \zeta (t^{\prime} )}{\rho }_{s}(0){U}_{-}(0,t)+{U}_{+}(t,0){\rho }_{s}(0)\displaystyle \frac{\delta {U}_{-}(0,t)}{\delta \zeta (t^{\prime} )}.\end{eqnarray} \tag{ 65 }$$

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

$$\begin{eqnarray}&&\displaystyle \frac{\delta {U}_{+}(t,t^{\prime} )}{\delta \zeta (\tau ^{\prime} )}=-i{\int }_{t^{\prime} }^{t}{U}_{+}(t,\tau )\displaystyle \frac{\delta {V}_{+}(\tau )}{\delta \zeta (\tau ^{\prime} )}{U}_{+}(\tau ,t^{\prime} )d\tau \end{eqnarray} \tag{ 66 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\delta {U}_{-}(t^{\prime} ,t)}{\delta \zeta (\tau ^{\prime} )}=i{\int }_{t^{\prime} }^{t}{U}_{-}(t^{\prime} ,\tau )\displaystyle \frac{\delta {V}_{-}(\tau )}{\delta \zeta (\tau ^{\prime} )}{U}_{-}(\tau ,t)d\tau ,\end{eqnarray} \tag{ 67 }$$

where ${V}_{\pm }(\tau )=-\left(\eta (\tau )\pm \tfrac{\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 η(t) and ν(t). Starting with η(t):

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\delta {U}_{+}(t,t^{\prime} )}{\delta \eta (\tau )} & = & -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} \\ & = & {{iU}}_{+}(t,\tau )f(\tau ){U}_{+}(\tau ,t^{\prime} ),\\ \displaystyle \frac{\delta {U}_{-}(t^{\prime} ,t)}{\delta \eta (\tau )} & = & 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} \\ & = & -{{iU}}_{-}(t^{\prime} ,\tau )f(\tau ){U}_{-}(\tau ,t).\end{array}\end{eqnarray} \tag{ 68 }$$

For ν(t):

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\delta {U}_{+}(t,t^{\prime} )}{\delta \nu (\tau )} & = & -i{\int }_{t^{\prime} }^{t}{U}_{+}(t,\tau ^{\prime} )\left[-\displaystyle \frac{1}{2}\delta (\tau -\tau ^{\prime} )f(\tau ^{\prime} )\right]{U}_{+}(\tau ^{\prime} ,t^{\prime} )d\tau ^{\prime} \\ & = & \displaystyle \frac{i}{2}{U}_{+}(t,\tau )f(\tau ){U}_{+}(\tau ,t^{\prime} ),\\ \displaystyle \frac{\delta {U}_{-}(t^{\prime} ,t)}{\delta \nu (\tau )} & = & i{\int }_{t^{\prime} }^{t}{U}_{-}(t^{\prime} ,\tau ^{\prime} )\left[\displaystyle \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).\end{array}\end{eqnarray} \tag{ 69 }$$

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

$$\begin{eqnarray}\begin{array}{rcl}\langle \eta (t){\rho }_{s}(t)\rangle & = & {\int }_{0}^{t}{dt}^{\prime} {K}^{\eta \eta }(t-t^{\prime} )\left\langle \displaystyle \frac{\delta {\rho }_{s}(t)}{\delta \eta (t^{\prime} )}\right\rangle +{\int }_{0}^{t}{dt}^{\prime} {K}^{\eta \nu }(t-t^{\prime} )\left\langle \displaystyle \frac{\delta {\rho }_{s}(t)}{\delta \nu (t^{\prime} )}\right\rangle \\ & = & {\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)\right.\\ & & \left.-{U}_{+}(t,0){\rho }_{s}(0){{iU}}_{-}(0,t^{\prime} )f(t^{\prime} ){U}_{-}(t^{\prime} ,t)\right\rangle \\ & & +{\int }_{0}^{t}{dt}^{\prime} {K}^{\eta \nu }(t-t^{\prime} )\left\langle \displaystyle \frac{i}{2}{U}_{+}(t,t^{\prime} )f(t^{\prime} ){U}_{+}(t^{\prime} ,0){\rho }_{s}(0){U}_{-}(0,t)\right.\\ & & \left.+{U}_{+}(t,0){\rho }_{s}(0)\displaystyle \frac{i}{2}{U}_{-}(0,t^{\prime} )f(t^{\prime} ){U}_{-}(t^{\prime} ,t)\right\rangle .\end{array}\end{eqnarray} \tag{ 70 }$$

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

We derive the equations of motion in cylindrical coordinates (r², r_z , ϕ) for general γ. With these coordinates, the x and y components of the coherence vector are given by

$$\begin{eqnarray}&&{r}_{x}=\sqrt{{r}^{2}-{r}_{z}^{2}}\cos \phi \qquad {r}_{y}=\sqrt{{r}^{2}-{r}_{z}^{2}}\sin \phi .\end{eqnarray} \tag{ 71 }$$

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

$$\begin{eqnarray}&&{{dr}}^{2}=\displaystyle \frac{\partial {r}^{2}}{\partial {r}_{z}}{{dr}}_{z}+\displaystyle \frac{1}{2}\displaystyle \frac{{\partial }^{2}{r}^{2}}{\partial {r}_{z}^{2}}{\left({{dr}}_{z}\right)}^{2}+\displaystyle \frac{\partial {r}^{2}}{\partial {r}_{y}}{{dr}}_{y}+\displaystyle \frac{1}{2}\displaystyle \frac{\partial {r}^{2}}{\partial {r}_{y}^{2}}{\left({{dr}}_{y}\right)}^{2}+\displaystyle \frac{\partial {r}^{2}}{\partial {r}_{x}}{{dr}}_{x}+\displaystyle \frac{1}{2}\displaystyle \frac{{\partial }^{2}{r}^{2}}{\partial {r}_{x}^{2}}{\left({{dr}}_{x}\right)}^{2}\end{eqnarray} \tag{ 72 }$$

$$\begin{eqnarray}&&d\phi =\displaystyle \frac{\partial \phi }{\partial {r}_{x}}{{dr}}_{x}+\displaystyle \frac{1}{2}\displaystyle \frac{{\partial }^{2}\phi }{\partial {r}_{x}^{2}}{\left({{dr}}_{x}\right)}^{2}+\displaystyle \frac{\partial \phi }{\partial {r}_{y}}{{dr}}_{y}+\displaystyle \frac{1}{2}\displaystyle \frac{{\partial }^{2}\phi }{\partial {r}_{y}^{2}}{\left({{dr}}_{y}\right)}^{2}+\displaystyle \frac{{\partial }^{2}\phi }{\partial {r}_{x}\partial {r}_{y}}{{dr}}_{x}{{dr}}_{y}.\end{eqnarray} \tag{ 73 }$$

The required derivatives are:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {r}^{2}}{\partial {r}_{x}}=2{r}_{x}\qquad \displaystyle \frac{\partial {r}^{2}}{\partial {r}_{y}}=2{r}_{y}\qquad \displaystyle \frac{\partial {r}^{2}}{\partial {r}_{z}}=2{r}_{z}\end{eqnarray} \tag{ 74 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{r}^{2}}{\partial {r}_{x}^{2}}=\displaystyle \frac{{\partial }^{2}{r}^{2}}{\partial {r}_{y}^{2}}=\displaystyle \frac{{\partial }^{2}{r}^{2}}{\partial {r}_{z}^{2}}=2\end{eqnarray} \tag{ 75 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \phi }{\partial {r}_{x}}=-\displaystyle \frac{{r}_{y}}{{r}_{x}^{2}+{r}_{y}^{2}}\qquad \displaystyle \frac{\partial \phi }{\partial {r}_{y}}=\displaystyle \frac{{r}_{x}}{{r}_{x}^{2}+{r}_{y}^{2}}\end{eqnarray} \tag{ 76 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}\phi }{\partial {r}_{x}^{2}}=\displaystyle \frac{2{r}_{x}{r}_{y}}{{\left({r}_{x}^{2}+{r}_{y}^{2}\right)}^{2}}\qquad \displaystyle \frac{{\partial }^{2}\phi }{\partial {r}_{y}^{2}}=-\displaystyle \frac{2{r}_{x}{r}_{y}}{{\left({r}_{x}^{2}+{r}_{y}^{2}\right)}^{2}}\qquad \displaystyle \frac{{\partial }^{2}\phi }{\partial {r}_{x}{r}_{y}}=\displaystyle \frac{{r}_{y}^{2}-{r}_{x}^{2}}{{\left({r}_{x}^{2}+{r}_{y}^{2}\right)}^{2}},\end{eqnarray} \tag{ 77 }$$

which allows us to simplify equations (72) and (73):

$$\begin{eqnarray}&&{{dr}}^{2}=2{r}_{z}{{dr}}_{z}+{\left({{dr}}_{z}\right)}^{2}+2{r}_{y}{{dr}}_{y}+{\left({{dr}}_{y}\right)}^{2}+2{r}_{x}{{dr}}_{x}+{\left({{dr}}_{x}\right)}^{2}\end{eqnarray} \tag{ 78 }$$

$$\begin{eqnarray}&&d\phi =-\displaystyle \frac{{r}_{y}}{{r}_{x}^{2}+{r}_{y}^{2}}{{dr}}_{x}+\displaystyle \frac{{r}_{x}{r}_{y}}{{\left({r}_{x}^{2}+{r}_{y}^{2}\right)}^{2}}{\left({{dr}}_{x}\right)}^{2}+\displaystyle \frac{{r}_{x}}{{r}_{x}^{2}+{r}_{y}^{2}}{{dr}}_{y}-\displaystyle \frac{{r}_{x}{r}_{y}}{{\left({r}_{x}^{2}+{r}_{y}^{2}\right)}^{2}}{\left({{dr}}_{y}\right)}^{2}+\displaystyle \frac{{r}_{y}^{2}-{r}_{x}^{2}}{{\left({r}_{x}^{2}+{r}_{y}^{2}\right)}^{2}}{{dr}}_{x}{{dr}}_{y}.\end{eqnarray} \tag{ 79 }$$

We can calculate ${\left({{dr}}_{x}\right)}^{2}$, ${\left({{dr}}_{y}\right)}^{2}$, ${\left({{dr}}_{z}\right)}^{2}$ and dr_x dr_y using equations (35)–(37) and the substitution dW_i dW_j = δ_ij dt, keeping linear terms in dt:

$$\begin{eqnarray}&&{\left({{dr}}_{x}\right)}^{2}={\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}\end{eqnarray} \tag{ 80 }$$

$$\begin{eqnarray}&&{\left({{dr}}_{y}\right)}^{2}=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}\end{eqnarray} \tag{ 81 }$$

$$\begin{eqnarray}&&{\left({{dr}}_{z}\right)}^{2}={\lambda }^{2}{r}_{x}^{2}{\left(\beta \epsilon +2{r}_{z}\right)}^{2}{dt}+{\lambda }^{2}{r}_{y}^{2}{\left(\beta \epsilon +2{r}_{z}\right)}^{2}{dt}+4{\lambda }^{2}{\gamma }^{2}{\left(1-{r}_{z}^{2}\right)}^{2}{dt}\end{eqnarray} \tag{ 82 }$$

$$\begin{eqnarray}&&{{dr}}_{x}{{dr}}_{y}=-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}.\end{eqnarray} \tag{ 83 }$$

After some algebra and dropping quadratic terms in β, we obtain

$$\begin{eqnarray}\begin{array}{rcl}{\left({{dr}}_{x}\right)}^{2} & = & 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-2{r}_{x}^{2}\right]{dt}\\ {\left({{dr}}_{y}\right)}^{2} & = & 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-2{r}_{y}^{2}\right]{dt}\\ {\left({{dr}}_{z}\right)}^{2} & = & 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}.\end{array}\end{eqnarray} \tag{ 84 }$$

The sum of equations (84) is given by

$$\begin{eqnarray}\begin{array}{rcl}{\left({{dr}}_{x}\right)}^{2}+{\left({{dr}}_{y}\right)}^{2}+{\left({{dr}}_{z}\right)}^{2} & = & 4{\lambda }^{2}\left[-2{\gamma }^{2}{r}_{z}^{2}+2\beta \epsilon {r}_{z}+{\gamma }^{2}\left(1+{r}_{z}^{4}\right)\right.\\ & & \left.+\left({r}^{2}-{r}_{z}^{2}\right)\left({r}^{2}+{\gamma }^{2}{r}_{z}^{2}-2\right)+2\right]{dt}.\end{array}\end{eqnarray} \tag{ 85 }$$

Moving on to the 2r_i dr_i terms:

$$\begin{eqnarray}\begin{array}{rcl}2{r}_{x}{{dr}}_{x} & = & -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}\\ 2{r}_{y}{{dr}}_{y} & = & 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}\\ 2{r}_{z}{{dr}}_{z} & = & -8{\lambda }^{2}{r}_{z}^{2}{dt}-8{\lambda }^{2}\beta \epsilon {r}_{z}{dt}-2\lambda {r}_{x}{r}_{z}\left(\beta \epsilon +2{r}_{z}\right){{dW}}_{x}-2\lambda {r}_{y}{r}_{z}\left(\beta \epsilon +2{r}_{z}\right){{dW}}_{y}+4\gamma \lambda {r}_{z}\left(1-{r}_{z}^{2}\right){{dW}}_{z}.\end{array}\end{eqnarray} \tag{ 86 }$$

Adding up equations (86) leads to

$$\begin{eqnarray}\begin{array}{rcl}2{r}_{x}{{dr}}_{x}+2{r}_{y}{{dr}}_{y}+2{r}_{z}{{dr}}_{z} & = & -4{\lambda }^{2}\left(2{r}_{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}\\ & & +4\lambda {r}_{y}\left(1-{r}^{2}\right){{dW}}_{y}+4\lambda \gamma {r}_{z}\left(1-{r}^{2}\right){{dW}}_{z}.\end{array}\end{eqnarray} \tag{ 87 }$$

By adding equations (85) and (87) we obtain the final expression for dr²:

$$\begin{eqnarray}\begin{array}{rcl}{{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-\displaystyle \frac{{r}_{z}^{2}}{{r}^{2}}}r\cos \phi \left(1-{r}^{2}\right){{dW}}_{x}\\ & & +4\lambda \sqrt{1-\displaystyle \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{array}\end{eqnarray} \tag{ 88 }$$

For the calculation of d ϕ, the contributions linear in dr_i are:

$$\begin{eqnarray}&&-\displaystyle \frac{{r}_{y}}{{r}_{x}^{2}+{r}_{y}^{2}}{{dr}}_{x}+\displaystyle \frac{{r}_{x}}{{r}_{x}^{2}+{r}_{y}^{2}}{{dr}}_{y}=\displaystyle \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].\end{eqnarray} \tag{ 89 }$$

The terms quadratic in dr_i lead to

$$\begin{eqnarray}&&\displaystyle \frac{{r}_{x}{r}_{y}}{{\left({r}_{x}^{2}+{r}_{y}^{2}\right)}^{2}}{\left({{dr}}_{x}\right)}^{2}-\displaystyle \frac{{r}_{x}{r}_{y}}{{\left({r}_{x}^{2}+{r}_{y}^{2}\right)}^{2}}{\left({{dr}}_{y}\right)}^{2}+\displaystyle \frac{{r}_{y}^{2}-{r}_{x}^{2}}{{\left({r}_{x}^{2}+{r}_{y}^{2}\right)}^{2}}{{dr}}_{x}{{dr}}_{y}=0,\end{eqnarray} \tag{ 90 }$$

allowing us to write the final expression for d ϕ as

$$\begin{eqnarray}&&d\phi =2\epsilon {dt}-\lambda \displaystyle \frac{\left(\beta \epsilon {r}_{z}+2\right)\sin \phi }{\sqrt{{r}^{2}-{r}_{z}^{2}}}{{dW}}_{x}+\lambda \displaystyle \frac{\left(\beta \epsilon {r}_{z}+2\right)\cos \phi }{\sqrt{{r}^{2}-{r}_{z}^{2}}}{{dW}}_{y}.\end{eqnarray} \tag{ 91 }$$

We derive an analytical expression for d〈〈Δs_sys 〉〉, the noise- and coordinate-averaged, incremental stochastic system entropy production, written as

$$\begin{eqnarray}&&d\langle \langle {\rm{\Delta }}{s}_{{sys}}\rangle \rangle ={\int }_{-1}^{1}d\langle {\rm{\Delta }}{s}_{{sys}}\rangle \,p({r}_{z},t){{dr}}_{z},\end{eqnarray} \tag{ 92 }$$

where d〈Δs_sys 〉 is given by the noise averaged form of equation (11). We proceed as follows:

$$\begin{eqnarray}\begin{array}{rcl}d\langle \langle {\rm{\Delta }}{s}_{{sys}}\rangle \rangle & = & {\int }_{-1}^{1}\left(-\displaystyle \frac{\partial \mathrm{ln}p({r}_{z},t)}{\partial t}{dt}-\displaystyle \frac{\partial \mathrm{ln}p({r}_{z},t)}{\partial {r}_{z}}d\langle {r}_{z}\rangle -{D}_{{zz}}\displaystyle \frac{{\partial }^{2}\mathrm{ln}p({r}_{z},t)}{\partial {r}_{z}^{2}}{dt}\right)p({r}_{z},t){{dr}}_{z}\\ & = & {\int }_{-1}^{1}\left(-\displaystyle \frac{\partial \mathrm{ln}p({r}_{z},t)}{\partial t}{dt}-{A}_{z}\displaystyle \frac{\partial \mathrm{ln}p({r}_{z},t)}{\partial {r}_{z}}{dt}-{D}_{{zz}}\displaystyle \frac{{\partial }^{2}\mathrm{ln}p({r}_{z},t)}{\partial {r}_{z}^{2}}{dt}\right)p({r}_{z},t){{dr}}_{z}.\end{array}\end{eqnarray} \tag{ 93 }$$

Integrating the last term by parts, and introducing the r_z component of the probability current,

$$\begin{eqnarray}&&{J}_{z}={A}_{z}p({r}_{z},t)-\displaystyle \frac{\partial \left({D}_{{zz}}p({r}_{z},t)\right)}{\partial {r}_{z}},\end{eqnarray} \tag{ 94 }$$

that appears in the Fokker-Planck equation (8)

$$\begin{eqnarray}&&\displaystyle \frac{\partial p({r}_{z},t)}{\partial t}=-\displaystyle \frac{\partial {J}_{z}}{\partial {r}_{z}},\end{eqnarray} \tag{ 95 }$$

we can write:

$$\begin{eqnarray}\begin{array}{rcl}d\langle \langle {\rm{\Delta }}{s}_{{sys}}\rangle \rangle & = & {\int }_{-1}^{1}\left(\left(-\displaystyle \frac{\partial \mathrm{ln}p({r}_{z},t)}{\partial t}{dt}-{A}_{z}\displaystyle \frac{\partial \mathrm{ln}p({r}_{z},t)}{\partial {r}_{z}}{dt}\right)p({r}_{z},t)+\displaystyle \frac{\partial \left({D}_{{zz}}p({r}_{z},t)\right)}{\partial {r}_{z}}\displaystyle \frac{\partial \mathrm{ln}p({r}_{z},t)}{\partial {r}_{z}}{dt}\right){{dr}}_{z}\\ & & -{\left[{D}_{{zz}}\displaystyle \frac{\partial \mathrm{ln}p({r}_{z},t)}{\partial {r}_{z}}p({r}_{z},t)\right]}_{-1}^{1}{dt}\\ & = & {\int }_{-1}^{1}\left(-\displaystyle \frac{\partial \mathrm{ln}p({r}_{z},t)}{\partial t}p({r}_{z},t){dt}-{J}_{z}\displaystyle \frac{\partial \mathrm{ln}p({r}_{z},t)}{\partial {r}_{z}}{dt}\right){{dr}}_{z}-{\left[{D}_{{zz}}\displaystyle \frac{\partial p({r}_{z},t)}{\partial {r}_{z}}\right]}_{-1}^{1}{dt}.\end{array}\end{eqnarray} \tag{ 96 }$$

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):

$$\begin{eqnarray}\begin{array}{rcl}d\langle \langle {\rm{\Delta }}{s}_{{sys}}\rangle \rangle & = & {\int }_{-1}^{1}\left(-\displaystyle \frac{\partial \mathrm{ln}p({r}_{z},t)}{\partial t}p({r}_{z},t){dt}+\mathrm{ln}p({r}_{z},t)\displaystyle \frac{\partial {J}_{z}}{\partial {r}_{z}}{dt}\right){{dr}}_{z}-{\left[{J}_{z}\mathrm{ln}p({r}_{z},t)\right]}_{-1}^{1}{dt}\\ & & -{\left[{D}_{{zz}}\displaystyle \frac{\partial p({r}_{z},t)}{\partial {r}_{z}}\right]}_{-1}^{1}{dt}\\ & = & -{\int }_{-1}^{1}\left(\displaystyle \frac{\partial \mathrm{ln}p({r}_{z},t)}{\partial t}p({r}_{z},t){dt}+\mathrm{ln}p({r}_{z},t)\displaystyle \frac{\partial p({r}_{z},t)}{\partial t}{dt}\right){{dr}}_{z}-{\left[{J}_{z}\mathrm{ln}p({r}_{z},t)\right]}_{-1}^{1}{dt}\\ & & -{\left[{D}_{{zz}}\displaystyle \frac{\partial p({r}_{z},t)}{\partial {r}_{z}}\right]}_{-1}^{1}{dt}\\ & = & -\left(\displaystyle \frac{d}{{dt}}{\int }_{-1}^{1}p({r}_{z},t)\mathrm{ln}p({r}_{z},t){{dr}}_{z}\right){dt}-{\left[{J}_{z}\mathrm{ln}p({r}_{z},t)\right]}_{-1}^{1}{dt}-{\left[{D}_{{zz}}\displaystyle \frac{\partial p({r}_{z},t)}{\partial {r}_{z}}\right]}_{-1}^{1}{dt}\\ & = & {{dS}}_{G}(t)-{\left[{J}_{z}\mathrm{ln}p({r}_{z},t)\right]}_{-1}^{1}{dt}-{\left[{D}_{{zz}}\displaystyle \frac{\partial p({r}_{z},t)}{\partial {r}_{z}}\right]}_{-1}^{1}{dt},\end{array}\end{eqnarray} \tag{ 97 }$$

where dS_G is the increment of the Gibbs entropy, see equation (13).