# Quantum fluctuation theorem for dissipative processes

Gabriele de Chiara, Alberto Imparato

► To cite this version:

Gabriele de Chiara, Alberto Imparato. Quantum fluctuation theorem for dissipative processes. Phys.Rev.Res., 2022, 4 (2), pp.023230. 10.1103/PhysRevResearch.4.023230. hal-03708478

HAL Id: hal-03708478

https://hal.science/hal-03708478

Submitted on 13 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Quantum fluctuation theorem for dissipative processes

Gabriele De Chiara<sup>1,2</sup> and Alberto Imparato<sup>1,3</sup>

<sup>1</sup>Centre for Theoretical Atomic, Molecular and Optical Physics, Queen's University Belfast, Belfast BT7 1NN, United Kingdom

<sup>2</sup>Laboratoire Charles Coulomb (L2C), UMR 5221 CNRS-Université de Montpellier, F-34095 Montpellier, France

<sup>3</sup>Department of Physics and Astronomy, Aarhus University, 8000 Aarhus, Denmark

<figure><img src="image_2.png" /></figure>

(Received 18 August 2021; revised 13 January 2022; accepted 25 May 2022; published 21 June 2022)

We present a general quantum fluctuation theorem for the entropy production of an open quantum system coupled to multiple environments, not necessarily at equilibrium. Such a general theorem, when restricted to the weak-coupling and Markovian regime, holds for both local and global master equations, corroborating the thermodynamic consistency of local quantum master equations. The theorem is genuinely quantum, as it can be expressed in terms of conservation of a Hermitian operator, describing the dynamics of the system state operator and of the entropy change in the baths. The integral fluctuation theorem follows from the properties of such an operator. Furthermore, it is also valid when the system is described by a time-dependent Hamiltonian. As such, the quantum Jarzynski equality is a particular case of the general result presented here. Moreover, our result can be extended to nonthermal baths, as long as microreversibility is preserved. We present some numerical examples to showcase the exact results previously obtained. We finally generalize the fluctuation theorem to the case where the interaction between the system and the bath is explicitly taken into account. We show that the fluctuation theorem amounts to a relation between time-reversed dynamics of the global density matrix and a two-time correlation function along the forward dynamics involving the baths' entropy alone.

DOI: 10.1103/PhysRevResearch.4.023230

# I. INTRODUCTION

In classical stochastic thermodynamics the physics of work and heat fluctuations in small, out-of-equilibrium systems is now well understood in the framework of fluctuation the-orems (FTs) extending the second law of thermodynamics to the microscopic realm [1–7]. Their importance originates from their generality, relying on very few assumptions, and from providing connections between equilibrium quantities and fluctuations of entropy, heat, and work. While the proof of the FT in the classical regime often relies on the stochastic trajectories that a system performs in its phase space while in-teracting with the external environment [7], other approaches based on the symmetries of the classical master equation [8] or of the Fokker-Planck equation have been proposed [9,10].

In the quantum realm, several approaches have been put forward to generalize the FT and the impossibility of moni-toring a quantum system without disturbance has generated a long debate (see, for example, the reviews [11–13]). Several approaches to the quantum FT employ the formalism of quan-tum trajectories or quantum Monte Carlo (QMC) [14–22]. In this context, thermodynamic quantities, e.g., entropy and energy, can be sampled along quantum trajectories generated

by a continuous evolution and random quantum jumps. In-terestingly, the back-action from the observation of the jump, e.g., through the emission of a photon, gives rise to a gen-üinely quantum energy contribution, dubbed “quantum heat” [23–25]. Very often, the so-called two-point measurement scheme, where the system is observed at the start and end of the protocol, is used [26], potentially destroying useful quantum coherences. Alternative approaches exist [27–30]. One can also prove the FT by diagonalizing the instantaneous density matrix, thus generating a quantum counterpart of the classical trajectories [31]. A fully coherent quantum FT in the framework of quantum resource theory has also been proposed [32] and recently experimentally verified [33].

In this paper, we first present a general and unified approach that naturally extends the classical FT to open quan-tum systems and is only based on the quantum Lindblad master equation. We prove the FT by introducing an auxiliary quantum master equation and a Hermitian operator that ac-counts for the time evolution of both the system's state and the bath entropy. Our formalism, based on the change of entropy in the baths and in the system, is general: it is valid at all times and not just at steady state; it is valid for both local and global quantum master equations (for which a heated debate has arisen in recent years), for an arbitrary number of baths, and for time-dependent system Hamiltonians.

We then consider the fate of the FT when the interac-tion between the system of interest and its environment is explicitly taken into account, and the total dynamics for the combined system is unitary. We show that the FT in this case can be expressed in terms of a two-time correlation function

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
<figure><img src="image_3.png" /><figcaption>FIG. 1. Schematics of the setup: a system is interacting with $N_B$ baths which are continuously monitored by their respective measurement apparatus.</figcaption></figure>

involving the baths' entropy alone. Such a correlation function turns out to be equal to the transition probability of the system along the time-reversed dynamics. This result confirms the intimate connection between the FT and the lack of symmetry between the forward and backward dynamics.

Moreover, we validate the results obtained in the first part of the paper: starting from the unitary dynamics we rederive the auxiliary master equation for the system in the same limit where the Lindblad master equation holds. We discuss explicitly the requirements for the local detailed balance to hold (or not), in connection with the interaction mechanisms between the system and the baths, and with the baths' properties. This analysis provides an operational approach to evaluate the entropy change in the bath and thus to test experimentally or numerically the proposed FT.

Our quantum FT goes beyond the two-point measurement scheme: it requires an initial projection of the system in an arbitrary basis to estimate the initial system's entropy along a specific trajectory. The initial projection may also occur in the eigenbasis of the initial density matrix, thus preventing any measurement back-action. It also requires the continuous monitoring of the baths but not the final system's projection or the continuous measurement of the work done on the system for time-dependent Hamiltonians (see Fig. 1). As such, its experimental implementation may be easier than schemes based on two-point measurements.

Moreover, a dissipative quantum Jarzynski equality (JE) [34–37] can be derived in operator form, and the result of the classical stochastic thermodynamics can be immediately recovered for diagonal density matrices. Our approach is unifying as it provides a proof of the quantum FT for relevant physical situations, namely, manipulated systems in contact with multiple heat baths, and allows a full description of the system's thermodynamics

## II. THE FLUCTUATION THEOREM

We assume the evolution of a d-level system weakly coupled to $N_b$ environments to follow the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) master equation (ME) [38–40]

(ħ = 1):

$$ \frac{d\rho(t)}{dt} = -i[H(t), \rho(t)] + \sum_{\alpha=1}^{N_b} D_{\alpha}[\rho(t)], \quad (1) $$

for its density matrix $\rho(t)$ with dissipators

$$ D_{\alpha}[\rho] = \sum_{\lambda} \gamma_{\alpha, \lambda} \left( L_{\lambda} \rho L_{\lambda}^{\dagger} - \frac{1}{2} \{ L_{\lambda}^{\dagger} L_{\lambda}, \rho \} \right), \quad (2) $$

where $\alpha$ labels the environment and $L_\lambda = |j'\rangle\langle j|$ and $\lambda = \lambda(j \to j')$ denotes a transition between two states $|j\rangle$ and $|j'\rangle$. In the following, in order to lighten the notation, we will omit the initial and final states of such a transition. The orthonormal basis $\{|j\rangle\}$, though arbitrary, is motivated by the physical environments. If the states $|j\rangle$ are the eigenstates of $H$ the ME is dubbed “global” or “local” otherwise [41]. We split an operator $X = X_D + X_{ND}$ into the sum of its diagonal and nondiagonal parts in the basis $\{|j\rangle\}$ and set $X_{jj} = \langle j|X|j\rangle$. We consider the general case where different baths can drive the same transition $\lambda$, if the dissipation rates (potentially time dependent) $\gamma_{\alpha, \lambda} \neq 0$.

In the classical case the FT reads [7]

$$ \langle e^{-\Delta S_B - \Delta S_S} \rangle = 1, \quad (3) $$

where $\Delta S_S$ is the entropy change of the system along a given stochastic trajectory, and, if the detailed balance is fulfilled, the total entropy change in the baths due to the heat $Q_\alpha$ reads $\Delta S_B = -\sum_\alpha \beta_\alpha Q_\alpha$. In Eq. (3) the average runs over all the possible stochastic trajectories in the phase space, given the system dynamics and the time protocol $H(t)$. We adopt the convention $Q_\alpha > 0$ when the heat flows from the bath into the system. We will first derive a quantum FT without assuming the local detailed balance condition for the jump rates: such an assumption will be later introduced in the paper, and its consequences discussed.

The classical FT (3) requires the characterization of the bath entropy statistics. To this end, we will follow closely the approach discussed in Ref. [8] for classical stochastic systems. The “jumps” between two states in the system occur at the rates $\gamma_{jj}$ because of the interaction with the baths: we thus introduce the elementary current associated with the jump $\lambda$:

$$ \Delta s_{\alpha, j'j} = - \log(\gamma_{\alpha, j'j} / \gamma_{\alpha, jj}). \quad (4) $$

We introduce a total “jump” current $S$ given by the sum of the contributions (4) for all the baths and for all the jumps up to time $t$. Such a quantity is akin to the entropy change $\Delta S_B$ in Eq. (3): we will elaborate later on this connection. We would thus like to characterize the joint probability distribution $\Phi_j(S, t)$ of finding the system in the state $|j\rangle$ with a total jump current $S$ up to the time $t$. This can be done by introducing an extended quantum ME that takes into account the dynamic evolution of both the density operator and of the quantity $S$. In practice we introduce a modified density matrix $\bar{\rho}(S, t)$, such that its diagonal elements describe the desired joint probability $\bar{\rho}_{jj} = \Phi_j(S, t)$. It is possible to show that
the modified density matrix satisfies the modified ME:

$$ \partial_t \bar{\rho}_{j'j}(S, t) = -i[H, \bar{\rho}]_{j'j} + \sum_{\alpha, j} \{\gamma_{\alpha, j'} \times \left[ \sum_{n=0}^{\infty} \frac{(-\Delta s_{\alpha, j'}^n)^n}{n!} \frac{\partial^n \bar{\rho}_{jj}}{\partial \Delta S_B^n} \right] - \gamma_{\alpha, j'} \bar{\rho}_{j'j} \}, \quad (5) $$

$$ \partial_t \bar{\rho}_{lk}(S, t) = -i[H, \bar{\rho}]_{lk} + \sum_{\alpha} D_{ND,\alpha}[\bar{\rho}_{ND}]_{lk}, \quad l \neq k, \quad (6) $$

where $D_{ND,\alpha}[\cdot]$ is a superoperator such that $D_{ND}[\bar{\rho}_{ND}]$ has vanishing diagonal terms and only couples nondiagonal terms of the density matrix in the chosen basis. A derivation of Eqs. (5) and (6) based on the GKSL (1) alone is discussed in Appendix A. Such a derivation starts from the fact that, if $S$ is the total jump current at time $t$, after the transition $j \to j'$ such a current reads $S + \Delta s_{\alpha, j'j}$. Next, we introduce the generating function operator $\Psi(\xi, t)$ akin to the generating function in classical physics, defined as

$$ \Psi(\xi, t) = \int dS \bar{\rho}(S, t) e^{-\xi S}, \quad (7) $$

and applying this integral transform to both sides of Eqs. (5) and (6), one obtains an equation for $\Psi(\xi, t)$:

$$ \partial_t \Psi(\xi, t) = \mathcal{L}_{\xi}[\Psi(\xi, t)], \quad (8) $$

where the superoperator $\mathcal{L}_{\xi}[\cdot]$ depends parametrically on $\xi$, and its full expression is given in Appendix A [see Eqs. (A8) and (A9)]. We will see that the modified density matrix $\bar{\rho}(t)$, its generating function $\Psi(t)$, and their modified MEs are the essential ingredients to prove the quantum FT.

We now discuss the physical initial condition for Eqs. (5) and (6) and Eq. (8). Let $\rho^{(0)}$ be the initial state at $t=0$. Since no current $S$ has yet been generated, $\bar{\rho}(S, t=0) = \rho^{(0)}\delta(S)$, where $\delta(x)$ is the Dirac delta function. Given the definition of $\Psi$, Eq. (7), its initial condition becomes $\Psi(\xi, t=0) = \rho^{(0)}, \forall \xi$, and in the following $\Psi(\xi, t|\rho^{(0)})$ will indicate the solution to Eq. (8) with this specific initial condition. The following equalities hold:

$$ \langle e^{-\xi S} \rangle_{\pi_j, \rho^{(0)}} = \mathrm{Tr}[\pi_j \Psi(\xi, t|\rho^{(0)})] = \Psi_{jj}(\xi, t|\rho^{(0)}), \quad (9) $$

$$ \langle e^{-\xi S} \rangle_{\rho^{(0)}} = \mathrm{Tr}[\Psi(\xi, t|\rho^{(0)})], \quad (10) $$

where $\langle \cdots \rangle_{\pi_j, \rho^{(0)}}$ in Eq. (9) denotes the expectation value constrained by the initial state $\rho^{(0)}$, the state at time $t$ taken to be $\pi_j = |j\rangle\langle j|$, and in Eq. (10) we summed over all possible final states.

In the following we will mostly be interested in the case $\xi = 1$ to evaluate the expectation value of $\exp(-S)$. Letting $\Psi^{(1)}(t) = \Psi(\xi = 1, t)$, one finds

$$ \partial_t \Psi^{(1)}(t) = -i[H, \Psi^{(1)}(t)] + \sum_{\alpha} D_{\alpha}^{*}[\Psi^{(1)}(t)], \quad (11) $$

where $D_{\alpha}^{*}$ is the dual of the dissipator $D_{\alpha}$:

$$ D_{\alpha}^{*}[\cdot] = \sum_{\lambda} \gamma_{\alpha, \lambda} \left( L_{\lambda}^{\dagger} \cdot L_{\lambda} - \frac{1}{2} \{ L_{\lambda}^{\dagger} L_{\lambda}, \cdot \} \right). \quad (12) $$

The details of the derivation of Eq. (11) are given in Appendix A. Thus the time evolution of $\Psi^{(1)}$, as given by Eq. (11), has the same conservative part as in Eq. (1), but the dissipative part (the dual of $D_{\alpha}$) is the same as the one found in the time evolution of an operator in the Heisenberg picture: this observation reflects the fact that the operator $\Psi^{(1)}$ bears information on both the system state $\rho(t)$ and on the quantity $\exp(-S)$. As such, the operator $\Psi^{(1)}$ is Hermitian at any time if it is Hermitian at $t=0$. Modified MEs of the type (8) emerge quite naturally through the large deviation approach used to study the long-time limit of thermodynamic currents. Such an approach was first introduced in Ref. [42], and became later quite consolidated in the field of stochastic thermodynamics (see, e.g., Refs. [8,43,44]). A modified quantum ME was first introduced in Ref. [16] and later in Ref. [45] to study the long-time counting statistics in dissipative quantum systems. While the modified ME (11) above has been obtained by applying the integral transformation (7) on the modified ME, Eqs. (5) and (6), in Sec. III we will provide an alternative derivation of Eq. (11), based on the unitary evolution of a system explicitly interacting with a set of baths.

Let us now introduce the operator $\bar{\Psi}^{(1)}(t)$:

$$ \bar{\Psi}^{(1)}(t) = \sum_{j_0} \Psi^{(1)}(t|\pi_{j_0}). \quad (13) $$

Given that $\bar{\Psi}^{(1)}(t)$ is a linear combination of solutions of Eq. (8) with $\xi = 1$, it is a solution itself, with initial condition $\bar{\Psi}^{(1)}(0) = \sum_{j_0} \pi_{j_0} = \mathbb{I}$. Inspection of Eq. (11) suggests that $\bar{\Psi}^{(1)}(t)$ is a stationary solution at any time:

$$ \bar{\Psi}^{(1)}(t) = \bar{\Psi}^{(1)}(t|\mathbb{I}) = \mathbb{I}, \quad \forall t \ge 0. \quad (14) $$

The last equality for $\bar{\Psi}^{(1)}(t)$ is the first important result of the present paper. It is the operatorial counterpart of the integral theorem (3), and it involves the Hermitian operators $\Psi^{(1)}(t|\pi_{j_0})$ expressing the joint dynamics of the system density operator and of the total jump current $S$. It does not depend on the choice of the initial basis. Remarkably by introducing an arbitrary basis $\{|b_0\rangle\}$, and noticing that Eq. (8) is linear, we can write the solution at time $t$, with the specific initial condition $\pi_{b_0}$ as

$$ \Psi^{(1)}(t|\pi_{b_0}) = \sum_{j_0, j'_0} \langle j_0|b_0\rangle \langle b_0|j'_0\rangle \Psi^{(1)}(t||j_0\rangle\langle j'_0|). \quad (15) $$

From Eq. (15) one obtains $\sum_{b_0} \Psi^{(1)}(t|\pi_{b_0}) = \sum_{j_0} \Psi^{(1)}(t|\pi_{j_0}) = \bar{\Psi}^{(1)}(t)$, using $\sum_{b_0} |b_0\rangle\langle b_0| = \mathbb{I}$. Furthermore Eqs. (10), (13), and (14) imply $\sum_{j_0} (\exp(-\Delta S_B))_{\pi_{j_0}} = d$.

We now present the mathematical statement of the integral FT which is the second main result of this paper, and later on we discuss its physical significance. Such a FT reads

$$ \mathrm{Tr} \left[ \sum_{b_0} \rho_{b_0 b_0}^{(0)} \Psi^{(1)}(t | \pi_{b_0}) e^{-\log \rho_{b_0 b_0}^{(0)}} e^{\log \rho^{(f)}} \right] = \mathrm{Tr}[\bar{\Psi}^{(1)}(t) \rho^{(f)}] = 1, \quad (16) $$

where we have used Eqs. (13)–(15) and introduced an arbitrary but normalized final state $\rho^{(f)}$. This feature of the FT was already noticed in Ref. [7] for the classical case: Eq. (16) holds
for any normalized quantum final state $\rho^{(f)}$, not necessarily the solution of Eq. (1) at time $t$.

Equation (16) has the form of the expectation value of the Hermitian operator $\bar{\Psi}^{(1)}$ over the final state of the system. We then recall the physical interpretation of the operator $\Psi^{(1)}(t|\rho^{(0)})$: its jth diagonal element represents the expectation value of $\exp(-S)$ constrained by the initial $\rho^{(0)}$ and final state $\pi_j$, Eq. (9). We make no specific assumption on the coherence in the initial and the final states $\rho^{(0)}$ and $\rho^{(f)}$: they can both exhibit coherence in the jump operators' basis $\{|j\rangle\}$ as well as in the energy eigenbasis. Let $\{|k_0\rangle\}$ and $\{|r_f\rangle\}$ be the eigenbases of $\rho^{(0)}$ and $\rho^{(f)}$, respectively. Changing basis from $\{|j\rangle\}$ to $\{|r_f\rangle\}$ and recalling definition (7), we can rewrite Eq. (16) as

$$1 = \sum_{r_f, j_0} \langle e^{-S} \rangle_{\pi_{r_f}, \pi_{j_0}} p_{j_0}^{(0)} e^{-\log \rho_{j_0 j_0}^{(0)}} e^{\log \rho_{r_f r_f}^{(f)}} \\ = \sum_{r_f, k_0} \int dS e^{-S + \log \rho_{r_f r_f}^{(f)} - \log \rho_{k_0 k_0}^{(0)}} \Phi_{r_f}(S, t|\pi_{k_0}) p_{k_0}^{(0)}. \quad (17)$$

We notice that the last expression for the FT has the form of the classical counterpart Eq. (3), with an average over the joint probability distribution $\Phi_{r_f}(S)$. Moreover, Eqs. (16) ad (17) exhibit the structure of a double trace, one over the initial state $\rho^{(0)}$ and one over the final one $\rho^{(f)}$, and as such their values do not depend on the chosen basis. Notice that performing the initial projection of $\rho^{(0)}$ in its eigenbasis preserves its quantum coherence in other bases, e.g., in the energy eigenbasis. Using Jensen's inequality, one recovers the second law:

$$\langle S \rangle - (\mathrm{Tr}[\rho^{(f)} \log \rho^{(f)}] - \mathrm{Tr}[\rho^{(0)} \log \rho^{(0)}]) \ge 0. \quad (18)$$

We can now establish the relation between the dynamics of the system, as expressed by Eqs. (1) and (2) and the thermodynamics significance of our results in Eqs. (14), (16), and (17). Such a relation follows immediately if the dissipation rates $\gamma_{\alpha, \lambda}$ are taken to obey the local detailed balance condition (LDBC): $\gamma_{\alpha, \lambda}(\omega_{\lambda})/\gamma_{\alpha, \lambda}(-\omega_{\lambda}) = \exp(-\beta_{\alpha}\omega_{\lambda})$, where $\omega_{\lambda} = \omega_{j'j} = H_{j'j} - H_{jj} = \Delta H_{D,\lambda}$, and $\beta_{\alpha}$ is the inverse temperature of the bath $\alpha$. If the LDBC holds then $\Delta s_{\alpha, j'j}$ in Eq. (4) can be immediately interpreted as the entropy change in the bath, given that $-\omega_{j'j} = \beta_{\alpha}^{-1}\log(\gamma_{\alpha, j'j}/\gamma_{\alpha, jj'})$ is heat flowing into a bath as a consequence of a jump. Thus the jump current $S$ entering Eqs. (5)–(7) can be identified, in this case, with the total change of entropy in the baths, $S = \Delta S_B$. While the LDBC is a standard assumption that allows to relate the dynamics to the thermodynamics in systems in contact with thermal baths, both in the classical and in the quantum regime [5,22,42], our proof of the FT does not rely on it. When the LDBC holds, evaluating the energy change $-\omega_{j'j}$ in the bath is sufficient to evaluate the entropy change in it according to Eq. (4). In case the LDBC does not hold, the rates $\gamma_{\alpha, j'j}$ and their relation to the $\omega_{j'j}$ need be estimated *a priori*.

Assuming the LDBC, since the heat flowing into a bath as a consequence of a jump is $-\omega_{j'j}$, involving only the diagonal part of the system Hamiltonian, we have $\langle S \rangle = \langle \Delta S_B \rangle = -\sum_{\alpha} \beta_{\alpha} Q_{D,\alpha}$ and we can rewrite Eq. (18) as

$$\Delta S_S = -(\mathrm{Tr}[\rho^{(f)} \log \rho^{(f)}] - \mathrm{Tr}[\rho^{(0)} \log \rho^{(0)}]) \ge \sum_{\alpha} \beta_{\alpha} Q_{D,\alpha}. \quad (19)$$

This is in accordance with the findings of Ref. [41], where it was found that, for local MEs, only the diagonal part of the heat currents, defined as $\dot{Q}_{D,\alpha} = \mathrm{Tr}\{\rho D_{\alpha}^*[H_D]\}$, enters the differential version of the second law, $d_t S_S \ge \sum_{\alpha} \beta_{\alpha} \dot{Q}_{D,\alpha}$, and flows to the environments as shown in Appendix B in the Born-Markov approximation. Notice that both the diagonal and nondiagonal heat contributions enter the first law already at operatorial level as $d_t H = \sum_{\alpha} D_{\alpha}^*[H_D + H_{ND}]$ [41]. However, the nondiagonal part of the Hamiltonian, specifically the quantity $\dot{Q}_{ND,\alpha} = \mathrm{Tr}\{\rho D_{\alpha}^*[H_{ND}]\}$, enters the energy balance of the interaction Hamiltonian, once a microscopic model made of the system, the baths, and their interaction mechanism is considered as detailed in Appendix B. While, within the collisional model framework, this corresponds to the work done when switching on and off the interaction with the environment, how this current can be interpreted in autonomous systems is left for future investigations.

Comparing Eqs. (17) and (19) we reach the conclusion (and the third important result in this paper) that the change in entropy in the bath, within the local ME framework, is only determined by the diagonal part of the system's Hamiltonian, and only $H_D$ enters the FT (17) and the second law (19). When the ME is global, $H_D = H$, one recovers the standard definition $\dot{Q}_{D,\alpha} = \dot{Q}_{\alpha} = \mathrm{Tr}\{D_{\alpha}[\rho]H\}$ (see Ref. [46]).

In this respect, it is interesting to consider the case of absence of coherence at $t=0$ and $t$. This situation occurs when the ME (1) is global, and the initial state $\rho^{(0)}$ has no coherence (i.e., the system is classical at any time), or when one performs measurements of the system state at $t=0$ and $t$, as in the two-point measurement scheme. In this case, it makes sense to introduce the concept of classical trajectories starting from the state $|j_0\rangle$ at $t=0$ and ending in the state $|j\rangle$. From Eq. (16) or Eq. (17) one then recovers immediately the classical FT (3) with $\Delta S_B + \Delta S_S = -\sum_{\alpha} \beta_{\alpha} Q_{\alpha} + \log \rho_{jj}^{(f)} - \log \rho_{j_0j_0}^{(0)}$.

Our results also hold when the system's Hamiltonian and/or the dissipation rates $\gamma_{\alpha, \gamma}$, entering Eq. (4), are time dependent. In the special case of a single bath at inverse temperature $\beta$ and choosing the final and initial states to be the Gibbs states $\rho^{(f)} = \exp(-\beta H(t))/Z_t$ and $\rho^{(0)} = \exp(-\beta H(0))/Z_0$, Eq. (16) becomes

$$\mathrm{Tr} \left[ \sum_{j_0} \rho_{j_0 j_0}^{(0)} \Psi^{(1)}(t | \pi_{j_0}) e^{\beta H_{j_0 j_0}^{(0)}} e^{-\beta H(t)} \right] = Z_t/Z_0 = e^{-\beta \Delta F}, \quad (20)$$

where $\Delta F$ is the difference in equilibrium free energy between the final and the initial thermal states. Furthermore, by assuming the LDBC the last equation reduces to the JE. Indeed, following the same procedure that leads to Eq. (17) we find

$$\sum_{j, j_0} \int dQ e^{\beta(Q_D - H_{jj}(t) + H_{j_0 j_0}^{(0)})} \Phi_j(Q_D, t | \pi_{j_0}) p_{j_0}^{(0)} \\ = \int dW_D P(W_D) e^{-\beta W_D} = e^{-\beta \Delta F}, \quad (21)$$

where we have set $\Delta S_B = -\beta Q_D$, and have defined work as $W_D \equiv \Delta H_D - Q_D$. Yet, it is worth noting that requiring the local detailed balance for a time-dependent Hamiltonian as we do in Eq. (21) is equivalent to requiring that the jump rates
<figure><img src="image_4.png" /><figcaption>FIG. 2. QMC simulations of the two-spin system: $\langle \exp(-\Delta S_{\text{tot}}) \rangle_{\text{tr}}$ and $\langle \exp(-\Delta S_B) \rangle_{\text{tr}}$ as functions of $t$, with $\Delta S_{\text{tot}} = S + \Delta S_S$. Averages over $10^6$ trajectories, $g = 0.1$, $T_a = 1$. (a) Diagonal $\rho^{(0)}$, $J = 0$, $h = 0.2$, $T_b = 1$ (classical case at equilibrium). (b) Diagonal $\rho^{(0)}$, $J = 0.1$, $h = 0.2$, $T_b = 1.2$. (c) Nondiagonal $\rho^{(0)}$, $J = 0.1$, $h = 0.2$, $T_b = 1.2$. (d) Nondiagonal $\rho^{(0)}$, $J = 0.2$, $T_b = 1.2$, and time-dependent field $h(t) = 0.4t/t_f$, with $t_f = 15$. In panels (a)–(d) the LDBC is used, and thus $S = \Delta S_B$. (e) Nonthermal jump rates $\gamma_{\alpha,\lambda}$ that do not satisfy the LDBC, with diagonal $\rho^{(0)}$, and $J = 0.1$, $h = 0.2$. See Appendix C for further details on the numerical simulations.</figcaption></figure>

adjust instantaneously to the value of the energy level, which in turn holds only for slow driving.

### Numerical examples

While results (16) and (17) are exact, we resort to numerical simulations to exemplify them and gain physical insight into their implications. Specifically, we consider a system of two spin-1/2 particles, characterized by Pauli operators $\sigma_{x,y,z}$, with Hamiltonian $H = -J\sigma_{x,a}\sigma_{x,b} - h(\sigma_{z,a} + \sigma_{z,b})$. We consider both the cases where each spin is connected to an equilibrium reservoir at temperatures $T_a$ and $T_b$, respectively, and where the two baths are nonthermal, with the dissipation rates not satisfying the LDBC. The jump operators “flip” the individual spins: $L_\lambda = \sigma_{-,a} \otimes \mathbb{I}_b$ or $L_\lambda = \mathbb{I}_a \otimes \sigma_{-,b}$. When $J \neq 0$ the quantum ME (1) is local (see Appendix C for further details on the numerical simulations).

We numerically evolve the system's state using the QMC algorithm for the unraveling of the ME [47] provided by QUTIP [48,49]. The QMC algorithm evolves the system's state from $|\psi(0)\rangle$ to $|\psi(t)\rangle$ with an alternation of continuous dynamics and stochastic jumps [47]. Along each trajectory, we monitor the individual jumps, driven by one of the two baths, and collect the statistics for the generalized bath's entropy:

$$S = \sum_{l=1}^{n_a} \Delta s_{a,j_{l+1},j_l}(t_l) + \sum_{m=1}^{n_b} \Delta s_{b,j_{m+1},j_m}(t_m), \quad (22)$$

where $t_l$ is the time at which the $l$th jump (out of $n_a$ total jumps) induced by the bath $a$ occurs, and analogously for bath $b$. When the LDBC holds, the previous equation has the simple interpretation $S = \Delta S_B = -(\beta_a Q_{a,D} + \beta_b Q_{b,D})$.

We evaluate the quantum FT in the form of Eq. (16) where now the average is taken over the trajectories generated by QMC. Within QMC, one can evaluate the solution $\rho(t)$ to the ME (1) as $\rho(t) = \langle |\psi(t)\rangle \langle \psi(t)| \rangle_{\text{tr}}$, where $\langle \cdot \cdot \rangle_{\text{tr}}$ is the average over the QMC trajectories. Analogously, from the definition of $\Psi(\xi, t)$ in Eq. (7), and from the sampling of the generalized entropy (22) along the trajectories, we can evaluate $\Psi^{(1)}(t|\pi_{k_0})$, for any chosen initial basis $\{|k_0\rangle\}$, as discussed in Appendix C.

The system's entropy change $\Delta S_S$ can be also evaluated along the QMC trajectories as detailed in Appendix C, and we can thus evaluate $\langle \exp(-\Delta S_{\text{tot}}) \rangle_{\text{tr}}$, with $\Delta S_{\text{tot}} = S + \Delta S_S$.

We consider both the case where the baths are thermal, with rates $\gamma_{\alpha,\lambda}$ obeying the LDBC, and the case of rates not satisfying the LDBC. Our numerical results in Fig. 2 confirm

the quantum FT (16) for initially diagonal and nondiagonal states in contact with thermal baths at the same temperatures ($T_a = T_b$) or different ones. Moreover, in Fig. 2(d) we consider the case of a time-dependent Hamiltonian, where the external field is changed according to $h(t) = h_0 + (h_1 - h_0)t/t_f$. In Fig. 2(e) we show the results for nonthermal baths violating the LDBC. In the same figure, we also show that, for all the considered cases, $\langle \exp(-\Delta S_B) \rangle \ge 1$, thus demonstrating the relevance of the system entropy variation for the FT to hold.

### III. FT AND LOCAL DETAILED BALANCE FOR UNITARY DYNAMICS

The aim of this section is twofold. First we want to generalize the results or Eqs. (16) and (17) of the previous section to the case of unitary dynamics, in a setup made of the system of interest interacting with a set of baths. We will then rederive Eq. (11) governing the time evolution of the operator $\Psi^{(1)}$ under the same assumptions that lead the standard ME (1). This approach will provide the reader with a more physical intuition of the current $S$, that we introduced in the previous section, and that enters the FT in the form of Eqs. (16) and (17).

We consider a total system made of the system of interest interacting through a Hamiltonian $V$ with a set of baths with Hamiltonian $H_{B_a}$: the total Hamiltonian $\mathcal{H}$ thus reads

$$\mathcal{H} = H_B + H + V, \quad (23)$$

$$H_B = \sum_{\alpha} H_{B_a}, \quad (24)$$

$$V = \sum_{\alpha, \lambda} V_{\alpha, \lambda}, \quad (25)$$

where $V_{\alpha,\lambda}$ is the interaction Hamiltonian with the $\alpha$th bath. The strength of the bath system is taken to be arbitrary in the following discussion, the standard weak-coupling limit being considered only later in this section.

We assume that initially the system and the baths are in a product state

$$\rho(0) = \rho(0) \otimes \rho_B(0), \quad (26)$$

$$\rho_B(0) = \rho_{B_1}(0) \otimes \rho_{B_2}(0) \otimes \cdots \otimes \rho_{B_{N_b}}(0). \quad (27)$$
Furthermore, we introduce the operators

$$S_{B_{\alpha}} = -\log \rho_{B_{\alpha}}(0), \qquad (28)$$

$$S_B = \sum_{\alpha=1}^{N_b} S_{B_{\alpha}}, \qquad (29)$$

akin to an entropy operator for the single bath and for all the baths, respectively. The case where some of eigenvalues of the $\rho_{B_{\alpha}}(0)$ are zero does not pose an issue on the discussion below, as we will be interested in the exponential of $S_{B_{\alpha}}$. We do not assume, a priori, that the initial states of the baths $\rho_{B_{\alpha}}(0)$ are thermal, although we will consider this specific case at a later stage. We also assume that the baths’ Hamiltonians, and the interaction Hamiltonian $V$, are time independent while the system Hamiltonian can be time dependent as discussed later.

Having introduced the total system, we start the discussion by deriving a rather obvious result, which will, however, set the stage for later derivations. We introduce the time-reversal operator for the total system $\theta$, such that $[\theta, \mathcal{H}] = 0$ and $-i\mathcal{H}\theta = \theta i\mathcal{H}\theta$, the latter equality expressing the antilinearity of the operator $\theta$ [50].

We assume that the initial state for the total system is a pure state $|\mathcal{J}_0\rangle$ (eigenstate of some observable). Here and in the following $|\mathcal{J}_x\rangle$ indicates states of the total system and $\Pi_x = |\mathcal{J}_x\rangle\langle\mathcal{J}_x|$ the corresponding projectors, while $|j_x\rangle$ and $\pi_x = |j_x\rangle\langle j_x|$ are states and projectors of the system of interest only as in the previous section.

The probability to observe the state $|\mathcal{J}_f\rangle$ after the total system evolves until time $t$ is

$$P(\mathcal{J}_f, t | \mathcal{J}_0, t = 0) = \mathrm{Tr}\{\Pi_f U_f \Pi_0 U_t^\dagger \Pi_f\}, \qquad (30)$$

where the unitary operator $U_t$ reads

$$U_t = \mathcal{T} \exp \left[ -i \int_0^t dt' (H_B + H(t') + V) \right], \qquad (31)$$

and $\mathcal{T}$ is the time-ordering operator. We also introduce the unitary operator when the global system evolves under the time-reversed Hamiltonian $H_B + H(t - t') + V$:

$$U_{\tilde{t}} = \mathcal{T} \exp \left[ -i \int_{0}^{t} dt' (H_{B} + H(t-t') + V) \right]. \qquad (32)$$

We now assume we can prepare the total system in the state $|\tilde{\mathcal{J}}_f\rangle = \theta|\mathcal{J}_f\rangle$ and evaluate the probability to find the system in the state $|\tilde{\mathcal{J}}_0\rangle = \theta|\mathcal{J}_0\rangle$ after evolving under the time-reversed dynamics for a time $t$,

$$\begin{aligned} P(\tilde{\mathcal{J}}_0, t | \tilde{\mathcal{J}}_f, t = 0) &= \mathrm{Tr}\{\tilde{\Pi}_0 U_{\tilde{t}} \tilde{\Pi}_f U_t^\dagger \tilde{\Pi}_0\} \\ &= \mathrm{Tr}\{\Pi_0 U_t^\dagger \Pi_f U_t \Pi_0\} \\ &= P(\tilde{\mathcal{J}}_f, t | \tilde{\mathcal{J}}_0, t = 0), \end{aligned} \qquad (33)$$

where we have used the antilinearity of $\theta$. This result is not surprising and expresses the conservation of the probability under time reversal and for unitary dynamics.

We now imagine that only the system is accessible to our measures, and assume the system and the baths are initially (only initially) in a product state. Furthermore, we want to do thermodynamics, and introduce a temperature (or a set of temperatures); thus we assume the bath(s) to be in the thermal

state, $\rho_B(0)$. We are then interested in the probability

$$\begin{aligned} P(j_f, t | j_0) &= \mathrm{Tr}\{\pi_f \otimes \mathbb{I}_B U_f \pi_0 \otimes \rho_B(0) U_f^\dagger \pi_f \otimes \mathbb{I}_B\} \\ &= \mathrm{Tr}\{\pi_f \rho_s(t) \pi_f\}, \end{aligned} \qquad (34)$$

and where, here and in the following, we have omitted to indicate the initial time $t = 0$. We then prepare the total system in the time-reversed state $\tilde{\rho}(0) = \theta\pi_f \otimes \rho_B(0)\theta^{-1}$, and ask what is the probability that the system of interest is in the state $|\tilde{j}_0\rangle$ after a time $t$. We remind the reader that $\theta$ is the time-reversal operator for the whole system. A straightforward calculation gives

$$\begin{aligned} P(\tilde{j}_0, t | \tilde{j}_f) &= \mathrm{Tr}\{\tilde{\pi}_0 \otimes \mathbb{I}_B U_{\tilde{t}} \tilde{\pi}_f \otimes \tilde{\rho}_B(0) U_{\tilde{t}}^\dagger \tilde{\pi}_0 \otimes \mathbb{I}_B\} \\ &= \mathrm{Tr}\{\pi_0 \otimes \mathbb{I}_B U_t^\dagger \pi_f \otimes \rho_B(0) U_t \pi_0 \otimes \mathbb{I}_B\} \\ &\neq P(j_f, t | j_0). \end{aligned} \qquad (35)$$

Compared to the case discussed above in this section, we see that the lack of symmetry in the probability under time reversal is due to the fact that we are now considering states of the system alone. In other words, we are tracing over the bath, before evaluating the transition probability of the system: this is the origin of the irreversibility in the transition probability. The goal is now to express the transition probability (35) in terms of the forward dynamics, and find a relation between the probability for the forward and backward dynamics that can lead us to a FT, as is the case for classical dynamics [51].

In the following we will use the notation $\langle A \rangle_{j_f, j_0}$ to indicate the time evolution of an observable with the forward dynamics, with constraints on the initial and final states. We also notice that the two-time correlation of the observables A and B reads

$$\langle A(t)B(0) \rangle = \mathrm{Tr}[U_t^\dagger A U_t B \mathcal{O}(0)]. \qquad (36)$$

Using the definition of $S_B$ in Eq. (29) we consider the two-time correlation

$$\begin{aligned} \langle e^{-S_B(t)} e^{S_B(0)} \rangle_{j_f, j_0} &= \mathrm{Tr}\{\pi_f \otimes \rho_B(0) U_f \pi_0 \otimes \mathbb{I}_B U_t^\dagger \pi_f \otimes \mathbb{I}_B\} \\ &= P(\tilde{j}_0, t | \tilde{j}_f). \end{aligned} \qquad (37)$$

This is a generalization of the Crooks’ quantum fluctuation relation [52] to the case of arbitrary interaction strength between the system and the baths, and for an arbitrary number of baths. Such a generalization relates the probability of the time-reversed transition to the two-time correlation function of the exponential of the baths’ entropy. For the classical case, with each bath initially in a thermal state at temperature $\beta_{\alpha}$, Eq. (37) becomes

$$P(j_f, t | j_0) e^{-\Delta S_B} = P(\tilde{j}_0, t | \tilde{j}_f) \qquad (38)$$

and $\Delta S_B = \sum_{\alpha} \beta_{\alpha} \Delta H_{B_{\alpha}}$ is the change of entropy in the baths, along any system trajectory connecting $j_0$ to $j_f$.

From Eq. (37) we obtain immediately

$$\begin{aligned} &\sum_{j_0} \langle e^{-S_B(t)} e^{S_B(0)} \rangle_{j_f, j_0} \\ &= \sum_{j_0} \mathrm{Tr}\{\pi_f \otimes \rho_B(0) U_f \pi_0 \otimes \mathbb{I}_B U_t^\dagger \pi_f \otimes \mathbb{I}_B\} = 1. \end{aligned} \qquad (39)$$

We then continue along the lines of the previous section and consider the arbitrary basis $\{|b_0\rangle\}$ and arbitrary initial and
final states $\rho^{(0)}$ and $\rho^{(f)}$, and from Eq. (39) we obtain

$$ 
\begin{aligned}
& \sum_{b_0} \langle e^{\log \rho^{(f)}(t)} e^{-S_B(t)} e^{S_B(0)} e^{-\log \rho^{(0)}_{b_0, b_0}} \rangle_{b_0} \rho^{(0)}_{b_0, b_0} \\
&= \sum_{b_0} \mathrm{Tr} \{ \rho^{(f)} \otimes \rho_B(0) U_t \pi_{b_0} \otimes \mathbb{I}_B U_t^\dagger \pi_f \otimes \mathbb{I}_B \} \\
&= \sum_{b_0} P(\tilde{b}_0, t | \tilde{\rho}^{(f)}) = 1.
\end{aligned} 
\quad (40) $$

Equations (39) and (40) represent the FT for the joint unitary evolution of the system and the baths: they generalize Eqs. (16) and (17) obtained starting from the master equation of the system alone. This result for the global unitary dynamics involves the two-time correlation of the entropy in the bath, which is related to the backward dynamics as expressed by Eq. (37). Furthermore, Eq. (40) involves the two-time correlation function of the system entropy too. Yet an initial measurement on the system is required, represented by the projection onto the arbitrary basis state $|b_0\rangle$. As discussed in the previous section, while the final state $\rho^{(f)}$ is arbitrary, one can of course take the specific choice $\rho^{(f)} = \rho(t) = \mathrm{Tr}_{B\rho}(t)$, i.e., the final state of the system along the forward dynamics. In this case, the FT involves the change of entropy in the system.

In deriving the above results we explicitly assumed that the system Hamiltonian is time dependent, $H(t)$, while the bath and the interaction Hamiltonian are time independent. Thus, even when work is done on or by the system on the external agent changing the system Hamiltonian, the relevant quantities to consider in the FT are the correlation between the exponential of the baths' entropy, as appearing on the right-hand side (RHS) of Eqs. (39) and (40). As such, these equations do not require the monitoring of the system energy change, or of the work done on or by the system, that might disturb the quantum dynamics of the system. Such an approach, requiring the monitoring of the environment only, was already put forward in, e.g., Ref. [52].

### Connection to the modified ME

We now connect the fluctuation relations for the unitary dynamics, Eqs. (37)–(40), with the modified master equation (11) introduced and discussed in the previous section [see also Eq. (A9) in Appendix A].

For the interaction Hamiltonian between the system and the $\alpha$th bath, introduced in Eq. (25), we choose the general form

$$ V_{\alpha, \lambda} = g_{\alpha, \lambda} (L_{\lambda}^{\dagger} A_{\alpha, \lambda} + L_{\lambda} A_{\alpha, \lambda}^{\dagger}), \quad (41) $$

where $L_\lambda$ are system operators, and $A_{\alpha, \lambda}$ are bath operators.

In the following we make a number of assumptions on the baths' properties and on the system-baths interaction, that are usually introduced in the derivation of the standard quantum ME (1) [53]. In particular we make the requirement that the baths' initial states commute with their corresponding Hamiltonians $[H_{B_\alpha}, \rho_{B_\alpha}(0)] = 0$. Furthermore, we choose the bath operators to be eigenoperators of the bath entropy operator [54],

$$ [S_{B_\alpha}, A_{\alpha, \lambda}] = \Omega_{\alpha, \lambda} A_{\alpha, \lambda}. \quad (42) $$

It is always possible to express the interaction Hamiltonian (41) in terms of eigenoperators of $S_{B_\alpha}$. Indeed if $A_{\alpha, \lambda}$ are not eigenoperators of $S_{B_\alpha}$, by expressing them in terms of the eigenbasis of $S_{B_\alpha}$ one can write $V_\alpha$ in terms of new (rotated) operators $A'_{\alpha, \lambda}$, which are eigenoperators of $S_{B_\alpha}$. In order to ease the notation, in the following we consider the case where the system Hamiltonian is time independent; the derivation of the corresponding results for the time-dependent case follows the same lines discussed below.

By using the assumption in Eq. (42), one can show the following equality to hold:

$$ e^{-\frac{S_B}{2}} e^{-i t \mathcal{H}} = e^{-i t \bar{\mathcal{H}}} e^{-\frac{S_B}{2}} \quad (43) $$

with

$$ \bar{\mathcal{H}} = H_B + H + \bar{V} \quad (44) $$

and where

$$ 
\begin{aligned}
\bar{V} &= \sum_{\alpha, \lambda} \bar{V}_{\alpha, \lambda} \\
&= \sum_{\alpha, \lambda} g_{\alpha, \lambda} \left( e^{-\beta_{\alpha} \Omega_{\alpha, \lambda} / 2} L_{\lambda}^{\dagger} A_{\alpha, \lambda} + e^{\beta_{\alpha} \Omega_{\alpha, \lambda} / 2} L_{\lambda} A_{\alpha, \lambda}^{\dagger} \right)
\end{aligned} 
\quad (45) $$

is a modified, non-Hermitian interaction Hamiltonian (see Appendix D).

Let $\bar{U}_f$ be the nonunitary operator obtained from Eq. (31) with the substitution $V \to \bar{V}$. Then comparing Eq. (43) with Eq. (37), we finally find

$$ 
\begin{aligned}
P(\tilde{j}_0, t | \tilde{j}_f, t = 0) &= \langle e^{-S_B(t)} e^{S_B(0)} \rangle_{\tilde{j}_f, \tilde{j}_0} \\
&= \mathrm{Tr} \{ \pi_f \otimes \mathbb{I}_B \bar{U}_f \pi_0 \otimes \rho_B(0) \bar{U}_f^\dagger \pi_f \otimes \mathbb{I}_B \},
\end{aligned} 
\quad (46) $$

resembling Eq. (33) that expresses the conservation of probability under time reversal for the total system. However, in Eq. (46) the forward time evolution occurs with $\bar{U}_f$.

Let us now consider the evolution of the system operator

$$ \Psi(t + \tau) = \mathrm{Tr}_B \{ \bar{U}_\tau \Psi(t) \otimes \rho_B(0) \bar{U}_\tau^\dagger \}, \quad (47) $$

given an arbitrary (possibly mixed) initial state of the system $\Psi(t)$. The projection of such an operator on the state $j_f$ is the two-time correlation function after a time $t + \tau$, as introduced in Eqs. (37) and (46), given the system state at time $t$, $\Psi(t)$:

$$ \langle e^{-S_B(t+\tau)} e^{S_B(t)} \rangle_{j_f, \psi(t)} = \langle j_f | \Psi(t+\tau) | j_f \rangle. \quad (48) $$

We make the Born-Markov approximation valid in the weak-coupling limit, and assume that the global system is in a product state at any time $\rho = \Psi(t) \otimes \rho_B$, and that the baths' state is time independent, $\rho_B = \rho_B(0)$ [53,55]. Thus in Eq. (47) we take the expansion up to the second order in $\tau$: some extra care is needed as $\bar{V}$ [see Eq. (45)], and thus $\bar{\mathcal{H}}$, is not Hermitian. We thus obtain

$$ 
\begin{aligned}
\Psi(t + \tau) &= \mathrm{Tr}_B \{ \bar{U}_\tau \Psi(t) \otimes \rho_B \bar{U}_\tau^\dagger \\
&\approx \mathrm{Tr}_B \left\{ \left( \mathbb{I} - i \tau \bar{\mathcal{H}} - \frac{\tau^2}{2} \bar{V}^2 \right) \right. \\
&\qquad \left. \times \Psi(t) \otimes \rho_B \left( \mathbb{I} + i \tau \bar{\mathcal{H}}^\dagger - \frac{\tau^2}{2} (\bar{V}^\dagger)^2 \right) \right\}.
\end{aligned} 
\quad (49) $$
A lengthy but straightforward calculation gives

$$ \Psi(t + \tau) \approx \Psi(t) - i\tau[H, \Psi(t)] + \tau^2 \sum_{\alpha, \lambda} g_{\alpha, \lambda}^2 (e^{\beta_a \Omega_{\alpha, \lambda}} L_{\lambda} \Psi(t) L_{\lambda}^{\dagger} \langle A_{\alpha, \lambda} A_{\alpha, \lambda}^{\dagger} \rangle_B + e^{-\beta_a \Omega_{\alpha, \lambda}} L_{\lambda}^{\dagger} \Psi(t) L_{\lambda} \langle A_{\alpha, \lambda}^{\dagger} A_{\alpha, \lambda} \rangle_B) \\ - \frac{\tau^2}{2} \sum_{\alpha, \lambda} g_{\alpha, \lambda}^2 (\langle L_{\lambda}^{\dagger} L_{\lambda}, \Psi(t) \rangle \langle A_{\alpha, \lambda} A_{\alpha, \lambda}^{\dagger} \rangle_B + \langle L_{\lambda} L_{\lambda}^{\dagger}, \Psi(t) \rangle \langle A_{\alpha, \lambda}^{\dagger} A_{\alpha, \lambda} \rangle_B). \qquad (50) $$

We now introduce the rates

$$ \gamma_{\alpha}(j \to j') = \tau g_{\alpha, \lambda}^2 \langle A_{\alpha, \lambda} A_{\alpha, \lambda}^{\dagger} \rangle_B, \qquad (51) $$

$$ \gamma_{\alpha}(j' \to j) = \tau g_{\alpha, \lambda}^2 \langle A_{\alpha, \lambda}^{\dagger} A_{\alpha, \lambda} \rangle_B. \qquad (52) $$

We see that, by virtue of Eq. (42), $A_{\alpha, \lambda}A_{\alpha, \lambda}^{\dagger}$ and $A_{\alpha, \lambda}^{\dagger}A_{\alpha, \lambda}$ are projectors onto the eigenstates of $S_{B_{\alpha}}$, Eq. (28). Furthermore, by taking $A_{\alpha, \lambda} = |N'_{\alpha}\rangle\langle N_{\alpha}|$, we have that the quantity $\Omega_{\alpha, \lambda}$ introduced in Eq. (42) becomes an entropy gap $\Omega_{\alpha, \lambda} = -(s_{N'_{\alpha}} - s_{N_{\alpha}}) = \log r_{N'_{\alpha}} - \log r_{N_{\alpha}}$, where $r_{N_{\alpha}}$ are the eigenvalues of the bath state $\rho_{B_{\alpha}}$. Thus, we find that the jump rates for the system dynamics appearing both in the system ME (1) and in the modified ME (11) obey the following relation:

$$ \frac{\gamma_{\alpha}(j \to j')}{\gamma_{\alpha}(j' \to j)} = \frac{\langle A_{\alpha, \lambda} A_{\alpha, \lambda}^{\dagger} \rangle_B}{\langle A_{\alpha, \lambda}^{\dagger} A_{\alpha, \lambda} \rangle_B} = e^{\Omega_{\alpha, \lambda}}. \qquad (53) $$

We remind the reader that, given our choice of the interaction Hamiltonian (41) and the jump operators $L_{\lambda}$ and $A_{\alpha, \lambda}$, a jump $\lambda$ in the system $j \to j'$ corresponds to a jump in the bath $N'_{\alpha} \to N_{\alpha}$. Equation (53) determines the asymmetric part of the system rates, and involves the entropy change in the bath for a transition ($j \to j'$, $N'_{\alpha} \to N_{\alpha}$), or equivalently the ratio of the population between the two states $N_{\alpha}$ and $N'_{\alpha}$ in the initial bath state. Neither the entropy nor the energy of the system appears in this relation.

Let $\epsilon_j = H_{jj}$ be the diagonal elements of the system Hamiltonian in the basis $\{|j\rangle\}$. Thus the quantities $\epsilon_j$ are eigenvalues of $H$ only if the chosen basis is an eigenbasis of $H$. Equation (53) becomes the standard local detailed balance condition involving the system energy gaps introduced in Sec. II, $\gamma_{\alpha}(j \to j')/\gamma_{\alpha}(j' \to j) = \exp(-\beta_{\alpha}\omega_{\lambda})$, only under two additional conditions: (i) the baths are initially in the thermal state, thus $\Omega_{\alpha, \lambda} = -\beta_{\alpha}(E_{N'_{\alpha}} - E_{N_{\alpha}})$, where $E_{N_{\alpha}}$ are the eigenvalues of $H_{B_{\alpha}}$, and (ii) the energy gaps in the baths and in the system are resonant, i.e., $E_{N'_{\alpha}} - E_{N_{\alpha}} = \omega_{\lambda} = H_{j'j} - H_{jj}$. We call such a condition local, because we envisage the typical physical setup where a bath is locally connected to a subpart of the system, possibly a single particle, through some mechanism as embodied by the interaction Hamiltonian (41). Condition (i) reminds us that the temperature (or the energy scale $k_B T_{\alpha}$) is a physical quantity associated with a bath at thermal equilibrium. Condition (ii) leads to the global detailed balance condition if, within the weak-coupling assumptions, no energy is stored in the coupling mechanism. We see immediately that this is the case when $\{|j\rangle\}$ is an eigenbasis for $H$: condition (ii) is then equivalent to assume that, given a change of energy in the system, there is a corresponding change of energy in the bath, with no energy released by or stored in the interaction mechanism represented by the Hamiltonian $V$ (see also Appendix B).

Having settled the connection between the system dynamic evolution and the underlying thermodynamic processes, as represented by Eq. (53), we continue our analysis of Eq. (50). Using the definition of the rates in Eqs. (51) and (52), the detailed balance condition (53), and taking the limit $\tau \to 0$, we obtain

$$ \partial_t \Psi(t) = -i[H, \Psi(t)] + \sum_{\alpha} D_{\alpha}^{*} [\Psi(t)], \qquad (54) $$

i.e., the generalized ME Eq. (11) [see Eq. (A9) and Appendix D for the details]. Had we taken the standard unitary operator $U_t$ in Eq. (49) [i.e., $\bar{V}$ replaced by $V$ on the RHS of Eq. (48)], we would have obtained the standard master equation (1) in the same limit.

Thus comparing this last result with Eqs. (47) and (48), we conclude that the operator $\Psi^{(1)}(t)$ appearing in Eq. (11) describes the joint evolution of the system state, and of the two-time correlation of the baths' entropy. For a single transition ($j \to j'$, $N'_{\alpha} \to N_{\alpha}$), the change in the bath entropy is given by $\Omega_{\alpha, \lambda}$. We recall the definition of $\Psi^{(1)}(t) = \Psi(\xi = 1, t)$, Eq. (7), and thus the quantity $S$ introduced in Sec. II is the total entropy change of the baths alone along a trajectory of the total system, the contribution of a single jump being set by the condition in Eq. (53). If the baths are initially thermal, and conserve their state for the entire duration of the system dynamics, $\Omega_{\alpha, \lambda} = -\beta_{\alpha}(E_{N'_{\alpha}} - E_{N_{\alpha}})$, thus Eq. (53) implies that the quantity $S$ is the change in the baths' energy multiplied by the corresponding inverse temperature: $S = \sum_{\alpha} \beta_{\alpha} \Delta H_{\alpha}$.

Even if the system Hamiltonian is time dependent, i.e., work is done on or by the system, making the standard assumption that the baths are at equilibrium implies that the FT (16) entails only the change of energy in the baths, which is of course affected by the work done on the system.

## IV. CONCLUSIONS

We have proved a quantum fluctuation theorem valid for a quantum system evolving with a dissipative process and with a time-dependent Hamiltonian. Our theorem goes beyond the two-point measurement scheme by preserving quantum correlations created by the nonequilibrium dynamics. Specifically, our results are also valid for local master equations, and their physical consistency is restored when one discerns between the entropic and energetic balance between the baths, the system, and the interaction mechanism. Our theorem provides a convenient tool for the numerical and experimental exploration of irreversible quantum coherent thermodynamics.

We have also shown that the derivation of the FT based on the ME can be extended to the case of unitary dynamics, when the explicit interaction between the systems and the bath is taken into account. This approach has the advantage
of highlighting the quantity that needs to be measured in a possible experimental validation of the theorem: the two-time correlation function of the exponential of the baths' entropy. In the case of weak coupling between the system and the baths, under the standard assumptions that lead to the GKSL ME, such a correlation function becomes a function of the jump rates appearing in the ME.

Furthermore, starting from the joint unitary dynamics of the system and the environment has the advantage of immediately setting the connection between the system dynamics and the underlying thermodynamic processes. The ratio between forward and backward jump rates originally depends on the properties of the baths alone: whether this relation leads to the global, local, or even to no detailed balance condition depends on further assumptions on the baths' and system's properties. Yet the FT does not require any detailed balance condition for the jump rates appearing in the ME.

The experimental validation of our theorem requires the monitoring of the environments and the observation of the corresponding quantum jumps. This can be achieved, for instance, in circuit and cavity QED [56], trapped-ion setups [57], or noisy intermediate-scale quantum computers [58].

All data created during this research are openly available from QUB-Pure [59].

## ACKNOWLEDGMENTS

G.D.C. acknowledges support by the UK EPSRC EP/S02994X/1. The numerical results presented in this work were obtained at the Centre for Scientific Computing, Aarhus. A.I. gratefully acknowledges the financial support of the Faculty of Science and Technology at Aarhus University through a Sabbatical scholarship and the hospitality of the Quantum Technology group, the Centre for Theoretical Atomic, Molecular and Optical Physics and the School of Mathematics and Physics, during his stay at Queen's University Belfast. This research was supported in part by the International Centre for Theoretical Sciences (ICTS) through the online program "Classical and Quantum Transport Processes: Current State and Future Directions" (code: ICTS/ctqp2022/1).

## APPENDIX A: THE MODIFIED QUANTUM MASTER EQUATION

Here we discuss the derivation of Eq. (5) in the main text. Such a derivation follows the corresponding derivation of the modified ME for classical stochastic systems discussed in Ref. [8]. The time evolution of $\rho(t)$ is described by the GKSL ME (1). We notice that the diagonal terms of the density matrix $\rho(t)$ express the time evolution of the populations, i.e., the probability $p_j(t) = \rho_{jj}(t)$ of observing the system in the state $|j\rangle$ at time $t$. We are interested in the joint probability distribution $\Phi_j(S, t)$ of finding the system in the state $j$ with a total generalized entropy $S$ which has flowed into the reservoirs at time $t$ as a consequence of the jumps between states in the system. In the main text we have already introduced the modified density matrix $\bar{\rho}(S, t)$, such that its diagonal elements describe the desired joint probability $\bar{\rho}_{jj} = \Phi_j(S, t)$. We first observe that the dissipators in the ME, as defined in

Eq. (2), only couple diagonal elements of the density matrix with diagonal elements and nondiagonal elements with non-diagonal elements. The coupling between the diagonal and nondiagonal elements of $\rho$ occurs through the coherent part of the dynamics, expressed by the commutator in Eq. (1). Thus we can write

$$D_{\alpha}[\rho] = D_{D,\alpha}[\rho_D] + D_{ND,\alpha}[\rho_{ND}], \quad (A1)$$

where $D_D[\rho_D]$ and $D_{ND}[\rho_{ND}]$ contain only diagonal and nondiagonal entries in their matrix representation, respectively, and $D_D[\rho_{ND}] = D_{ND}[\rho_D] = 0$. This becomes evident when one rewrites Eq. (1) as

$$\frac{d\rho(t)}{dt} = -i[H(t), \rho(t)] + \sum_{\alpha, \lambda} \gamma_{\lambda, \alpha} (|j'\rangle\langle j|\rho|j\rangle\langle j'| \\ -\frac{1}{2}(\pi_j\rho + \rho\pi_j)), \quad (A2)$$

where we remind the reader that $\lambda = \lambda(j \to j')$ denotes a transition between two states $|j\rangle$ and $|j'\rangle$. Isolating the diagonal part of Eq. (A2), one obtains

$$\partial_t \rho_{j'j} = -i[H, \rho]_{j'j} + \sum_{\alpha,j} \{\gamma_{\alpha,j'j} \rho_{jj} - \gamma_{\alpha,j'j'} \rho_{j'j}\}. \quad (A3)$$

The jump between the states $|j\rangle$ and $|j'\rangle$ occurs with rate $\gamma_{\alpha}(j \to j') = \gamma_{\alpha,j'j}$, and the corresponding entropy change in the bath $\alpha$ is given by

$$\Delta s_{\alpha,j'j} = -\log(\gamma_{\alpha,j'j} / \gamma_{\alpha,j'j}), \quad (A4)$$

which corresponds to Eq. (4) in the main text. Such jumps in the ME (1) are described by $D_D[\rho_D]$ implicitly introduced in Eq. (A3); thus only the coupling between the diagonal elements in the ME contributes to the change in the bath generalized entropy $S$. After the jump $j \to j'$ the generalized entropy changes as $S \to S + \Delta s_{\alpha,j'j}$. Let us now suppose that we know the modified density matrix $\bar{\rho}(S, t)$ at time $t$, and let us consider the time evolution of its diagonal and nondiagonal parts. For the diagonal part, at time $t + \tau$, we have

$$\begin{aligned} \bar{\rho}_{j'j}(S, t + \tau) \approx & \bar{\rho}_{j'j}(S, t) - i\tau[H, \bar{\rho}(S)]_{j'j} \\ & + \tau \sum_{\alpha,j} \gamma_{\alpha,j'j} \bar{\rho}_{jj}(S - \Delta s_{\alpha,j'j}, t) \\ & - \gamma_{\alpha,j'j} \bar{\rho}_{j'j}(S, t). \end{aligned} \quad (A5)$$

For the nondiagonal part of the dynamics, that does not contribute to $S$, we have

$$\begin{aligned} \bar{\rho}_{lk}(S, t + \tau) \approx & \bar{\rho}_{lk}(S, t) - i\tau[H, \bar{\rho}(S)]_{lk} \\ & + \tau \sum_{\alpha} D_{ND}[\bar{\rho}_{ND}(S, t)]_{lk}, \quad l \neq k. \end{aligned} \quad (A6)$$

Taking the limit $\tau \to 0$ in both Eqs. (A5) and (A6), and expanding the right-hand side of Eq. (A5) in Taylor series of the microscopic entropy change $\Delta s_{\alpha,j'j}$, we obtain the modified quantum ME, Eqs. (5) and (6) in the main text.

We introduce the operator $\Psi(\xi, t)$ defined as

$$\Psi(\xi, t) = \int dS \bar{\rho}(S, t)e^{-\xi S}, \quad (A7)$$
and, applying this integral transform to both sides of Eqs. (5) and (6) in the main text, one obtains the equations for $\Psi(\xi, t)$

$$ \partial_t \Psi_{j'j}(\xi, t) = -i[H, \Psi]_{j'j} + \sum_{\alpha, j} \left\{ \gamma_{\alpha, j'j} \left( \frac{\gamma_{\alpha, j'j}}{\gamma_{\alpha, j'j}} \right)^{\xi} \Psi_{jj}(\xi, t) - \gamma_{\alpha, j'j} \Psi_{j'j}(\xi, t) \right\}, \quad (A8) $$

$$ \partial_t \Psi_{lk}(\xi, t) = -i[H, \Psi]_{lk} + D_{ND}[\Psi_{ND}]_{lk}, \quad l \neq k. \quad (A9) $$

We now take $\xi = 1$ in Eqs. (A8) and (A9), yielding

$$ \partial_t \Psi_{j'j}^{(1)}(t) = -i[H, \Psi^{(1)}]_{j'j} + \sum_{\alpha, j} \gamma_{\alpha, j'j} \{\Psi_{jj}^{(1)}(t) - \Psi_{j'j}^{(1)}(t)\}, \quad (A10) $$

$$ \partial_t \Psi_{lk}^{(1)}(t) = -i[H, \Psi^{(1)}]_{lk} + D_{ND}[\Psi_{ND}^{(1)}]_{lk}, \quad l \neq k. \quad (A11) $$

The dual of the dissipator $D_\alpha[\cdot]$ applied on an operator X reads

$$ D_{\alpha}^{*}[X] = \sum_{j \to j''} \gamma_{\alpha, j''j} \left( X_{j''j''} \pi_j - \frac{1}{2} \{\pi_j, X\} \right), \quad (A12) $$

and considering the diagonal element

$$ D_{\alpha}^{*}[X]_{j'j} = \sum_{j''} \gamma_{\alpha, j''j} (X_{j''j''} - X_{j'j}), \quad (A13) $$

one finds that it corresponds to the term on the right-hand side of Eq. (A10). For the nondiagonal part, one finds $D_{\alpha,ND}^*[X_{ND}] = D_{\alpha,ND}[X_{ND}]$. All in all, one finds that $\Psi^{(1)}(t)$ obeys Eq. (11) in the main text.

## APPENDIX B: ON THE HEAT EXCHANGED BETWEEN THE SYSTEM AND THE BATH(S)

As in the main text, we consider the total system made of the system of interest with Hamiltonian $H$ and a set of baths with Hamiltonians $H_{B_\alpha}$, interacting through a Hamiltonian $V$ [see Eqs. (23)–(25)]. Given an arbitrary orthonormal basis $\{|j\rangle\}$ for the system, we split its Hamiltonian into its diagonal and nondiagonal parts,

$$ H = H_D + H_{ND}. \quad (B1) $$

The time evolution of the total density matrix is thus given by

$$ \rho(t + \tau) = e^{-i\tau\mathcal{H}}\rho(t)e^{i\tau\mathcal{H}}, \quad (B2) $$

for any $t$ and $\tau$. To lighten the notation, here and in the following we assume that all the above Hamiltonians, in particular $V$, are time independent. Here we do not consider the case where $H$ depends explicitly on the time, since in this section we are only interested in the energy exchange between the system and the baths, and not on the work done on the system by an external agent. However, the case of a time-dependent $H(t)$ only requires a straightforward modification. As discussed in the main text, we do not assume a priori that the states $|j\rangle$ are eigenstates of the system Hamiltonian.

We recall the definition of the interacting Hamiltonians introduced in Eq. (25), $V_{\alpha,\lambda} = g_{\alpha,\lambda}(L_\lambda^\dagger A_{\alpha,\lambda} + L_\lambda A_{\alpha,\lambda}^\dagger)$, with the system operators $L_\lambda$ and the bath operators $A_{\alpha,\lambda}$

satisfying

$$ [H_D, L_\lambda] = \omega_\lambda L_\lambda, \quad (B3) $$

$$ [H_{B_\alpha}, A_{\alpha,\lambda}] = \Omega_{\alpha,\lambda} A_{\alpha,\lambda}, \quad (B4) $$

with

$$ \omega_\lambda = \omega_{j'j} = \langle j'|H|j'\rangle - \langle j|H|j\rangle = H_{D,j'j} - H_{D,jj}, \quad (B5) $$

$$ \Omega_{\alpha,\lambda} = \langle N'_{\alpha} | H_{B_{\alpha}} | N'_{\alpha} \rangle - \langle N_{\alpha} | H_{B_{\alpha}} | N_{\alpha} \rangle = E_{N'_{\alpha}} - E_{N_{\alpha}}. \quad (B6) $$

Equation (B3) follows immediately from the choice of the jump operators $L_\lambda = |j'\rangle\langle j|$, while in this Appendix we choose the operators $A_{\alpha,\lambda}$ to be eigenoperators of the $\alpha$th Hamiltonian, with $A_{\alpha,\lambda} = |N'_{\alpha}\rangle\langle N_{\alpha}|$, and $|N_{\alpha}\rangle$ the corresponding energy eigenstates. We also notice that the choice of the interaction Hamiltonian, Eq. (25), implies that different baths can induce the same transition $j \to j'$, as long as $g_{\alpha,\lambda} \neq 0$. Physically, this corresponds to the case where a bath can be connected to more than a single subpart of the system, and is therefore more general than the case where one has the same number of baths and subparts, each bath only inducing transitions in the corresponding subpart of the system. Yet, we deem the approach as “local,” as the eigenkets defining the jump operators are not taken a priori to be eigenkets of $H_S$.

We can now study the energy balance for the baths and the system. First we notice that, since the total Hamiltonian is time independent, the total energy is always conserved for any $t$ and $\tau$:

$$ \Delta E_{\text{tot}}(t, t + \tau) = \text{Tr}[\mathcal{H}(\rho(t + \tau) - \rho(t))] = 0. \quad (B7) $$

If we define

$$ \Delta E_x = \langle H_x \rangle_{t+\tau} - \langle H_x \rangle_t = \text{Tr}[(e^{i\tau\mathcal{H}}H_x e^{-i\tau\mathcal{H}} - H_x)\rho(t)], \quad (B8) $$

with $x = B, S, V$, energy conservation imposes

$$ \Delta E_B + \Delta E_S = -\Delta E_V. \quad (B9) $$

We have

$$ e^{i\tau\mathcal{H}}H_x e^{-i\tau\mathcal{H}} = H_x + i\tau[\mathcal{H}, H_x] + \frac{(i\tau)^2}{2}[\mathcal{H}, [\mathcal{H}, H_x]] + \dots \quad (B10) $$

Let us first consider $E_B(t + \tau) - E_B(t)$ up to the second order in $\tau$. From Eq. (B10) we see that, in order to evaluate the first-order contribution to the baths’ energy flow, we need to evaluate $[\mathcal{H}, H_B]$. To this end, we introduce the operator

$$ \begin{aligned} V'_B &= [\mathcal{H}, H_B] = [V, H_B] \\ &= -\sum_{\alpha, \lambda} g_{\alpha, \lambda} \Omega_{\alpha, \lambda} \lambda (L_{\lambda}^{\dagger} A_{\alpha, \lambda} - L_{\lambda} A_{\alpha, \lambda}^{\dagger}). \end{aligned} \quad (B11) $$

To calculate the second-order term in $\Delta E_B$ we need to calculate

$$ [\mathcal{H}, V'_B] = \left[ \sum_{\alpha} H_{B_{\alpha}} + H_S + V, V'_B \right], \quad (B12) $$
where the individual commutators read

$$[H, V'_B] = \sum_{\alpha, \lambda} g_{\alpha, \lambda} \Omega_{\alpha, \lambda} (\omega_{\lambda} (L_{\lambda}^{\dagger} A_{\alpha, \lambda} + L_{\lambda} A_{\alpha, \lambda}^{\dagger}) - ([H_{ND}, L_{\lambda}^{\dagger}] A_{\alpha, \lambda} - [H_{ND}, L_{\lambda}] A_{\alpha, \lambda}^{\dagger})), \quad (B13)$$

$$\left[\sum_{\alpha} H_{B_{\alpha}}, \tilde{V}_{B}\right] = -\sum_{\alpha, \lambda} g_{\alpha, \lambda} \Omega_{\alpha, \lambda}^{2} (L_{\lambda}^{\dagger} A_{\alpha, \lambda} + L_{\lambda} A_{\alpha, \lambda}^{\dagger}), \quad (B14)$$

$$[V, \tilde{V}_B] = - \sum_{\alpha} \sum_{\lambda, \mu} [g_{\alpha, \lambda} (L_{\lambda}^{\dagger} A_{\alpha, \lambda} + L_{\lambda} A_{\alpha, \lambda}^{\dagger}), g_{\alpha, \mu} \Omega_{\alpha, \mu} \times (L_{\mu}^{\dagger} A_{\alpha, \mu} - L_{\mu} A_{\alpha, \mu}^{\dagger})], \quad (B15)$$

where $\lambda$ and $\mu$ are different transition indices. For future reference, we notice that only the last commutator contains quadratic terms of the types $A_{\alpha, \lambda}A_{\alpha, \lambda}^{\dagger}$ or $A_{\alpha, \lambda}^{\dagger}A_{\alpha, \lambda}$.

We now turn our attention to the system energy balance. Let us introduce the operators

$$V'_{D} = [V, H_{D}] = \sum_{\alpha, \lambda} g_{\alpha, \lambda} \omega_{\lambda} (L_{\lambda}^{\dagger} A_{\alpha, \lambda} - L_{\lambda} A_{\alpha, \lambda}^{\dagger}), \quad (B16)$$

$$V'_{ND} = [V, H_{ND}] = \sum_{\alpha, \lambda} g_{\alpha, \lambda} ([L_{\lambda}^{\dagger}, H_{ND}] A_{\alpha, \lambda} + [L_{\lambda}, H_{ND}] A_{\alpha, \lambda}^{\dagger}). \quad (B17)$$

From Eqs. (B8)–(B10), we see that to the first order in $\tau$ we need to evaluate

$$[\mathcal{H}, H] = [V, H_S] = V'_D + V'_{ND}. \quad (B18)$$

To calculate the second-order term in Eq. (B10) we need the commutator

$$[\mathcal{H}, V'_D + V'_{ND}] = [H_B + H_D + H_{ND} + V, V'_D + V'_{ND}].$$

We have

$$[H_B, V'_D] = \sum_{\lambda, \lambda} g_{\alpha, \lambda} \omega_{\lambda} \Omega_{\alpha, \lambda} (L_{\lambda}^{\dagger} A_{\alpha, \lambda} + L_{\lambda} A_{\alpha, \lambda}^{\dagger}),$$

$$[H_B, V'_{ND}] = -\sum_{\alpha, \lambda} g_{\alpha, \lambda} \Omega_{\alpha, \lambda} ([H_{ND}, L_{\lambda}^{\dagger}] A_{\alpha, \lambda} - [H_{ND}, L_{\lambda}] A_{\alpha, \lambda}^{\dagger}),$$

$$[H_S, V'_D] = \sum_{\alpha, \lambda} g_{\alpha, \lambda} \omega_{\lambda} (-\omega_{\lambda} (L_{\lambda}^{\dagger} A_{\alpha, \lambda} + L_{\lambda} A_{\alpha, \lambda}^{\dagger}) + [H_{ND}, L_{\lambda}^{\dagger}] A_{\alpha, \lambda} - [H_{ND}, L_{\lambda}] A_{\alpha, \lambda}^{\dagger}),$$

$$[V, V'_D] = \sum_{\alpha} \sum_{\lambda, \mu} [g_{\alpha, \lambda} (L_{\lambda}^{\dagger} A_{\alpha, \lambda} + L_{\lambda} A_{\alpha, \lambda}^{\dagger}), g_{\alpha, \mu} \omega_{\mu} (L_{\mu}^{\dagger} A_{\alpha, \mu} - L_{\mu} A_{\alpha, \mu}^{\dagger})], \quad (B19)$$

$$[V, V'_{ND}] = \sum_{\alpha} \sum_{\lambda, \mu} [g_{\alpha, \lambda} (L_{\lambda}^{\dagger} A_{\alpha, \lambda} + L_{\lambda} A_{\alpha, \lambda}^{\dagger}), \times g_{\alpha, \mu} ([L_{\mu}^{\dagger}, H_{ND}] A_{\alpha, \mu} + [L_{\mu}, H_{ND}] A_{\alpha, \mu}^{\dagger})]. \quad (B20)$$

We notice that the only commutators that are quadratic in the bath operators $A_\lambda$ and $A_\lambda^\dagger$ are the last two. In accordance with the Born-Markov approximation used to derive the Markovian master equation [55], we assume that the total density matrix is factorized at any time, $\rho_{\text{tot}}(t) = \rho_S(t) \otimes \rho_B(t)$, that the bath density matrix is time independent, $\rho_B(t) = \bigotimes_{\alpha=1}^{N_b} \rho_{B_\alpha}$, and that $[\rho_{B_\alpha}, H_{B_\alpha}] = 0$ even though $\rho_{B_\alpha}$ needs not be an equilibrium state. Without these assumptions the master equation would not be valid or would be modified.

Armed with these assumptions on the baths' state and operators, we can finally calculate the energy flow for the bath and for the system. From Eqs. (B8)–(B15), we have

$$\begin{aligned} \Delta E_B &= i\tau \langle V'_B \rangle - \frac{\tau^2}{2} \langle [\mathcal{H}, V'_B] \rangle \\ &= -\tau^2 \sum_{\alpha, \lambda} g_{\alpha, \lambda}^2 \Omega_{\alpha, \lambda} (\langle L_{\lambda}^{\dagger} L_{\lambda} \rangle_{S,t} \langle A_{\alpha, \lambda} A_{\alpha, \lambda}^{\dagger} \rangle_B - \langle L_{\lambda} L_{\lambda}^{\dagger} \rangle_{S,t} \langle A_{\alpha, \lambda}^{\dagger} A_{\alpha, \lambda} \rangle_B) + O(\tau^4), \end{aligned} \quad (B21)$$

where $\langle \cdot \rangle_{S,t} = \text{Tr}[\cdot \rho(t)]$ and $\langle \cdot \rangle_B = \text{Tr}[\cdot \rho_B]$. Similarly from Eqs. (B16)–(B20) we obtain

$$\begin{aligned} \Delta E_S = \tau^2 &\sum_{\alpha, \lambda} g_{\alpha, \lambda}^2 \omega_{\lambda} (\langle L_{\lambda}^{\dagger} L_{\lambda} \rangle_{S,t} \langle A_{\alpha, \lambda} A_{\alpha, \lambda}^{\dagger} \rangle_B - \langle L_{\lambda} L_{\lambda}^{\dagger} \rangle_{S,t} \langle A_{\alpha, \lambda}^{\dagger} A_{\alpha, \lambda} \rangle_B) \\ &+ \frac{g_{\alpha, \lambda}^2}{2} (\langle L_{\lambda}^{\dagger} [H_{ND}, L_{\lambda}] + [L_{\lambda}^{\dagger}, H_{ND}] L_{\lambda} \rangle_{S,t}) \langle A_{\alpha, \lambda} A_{\alpha, \lambda}^{\dagger} \rangle_B \\ &+ \frac{g_{\alpha, \lambda}^2}{2} (\langle L_{\lambda} [H_{ND}, L_{\lambda}^{\dagger}] + [L_{\lambda}, H_{ND}] L_{\lambda}^{\dagger} \rangle_{S,t}) \langle A_{\alpha, \lambda}^{\dagger} A_{\alpha, \lambda} \rangle_B + O(\tau^4). \end{aligned} \quad (B22)$$

A similar approach for the interaction Hamiltonian leads to

$$\begin{aligned} \Delta E_V = -\tau^2 &\sum_{\alpha, \lambda} \frac{g_{\alpha, \lambda}^2}{2} [2(\omega_{\alpha, \lambda} - \Omega_{\alpha, \lambda})(\langle L_{\alpha, \lambda}^{\dagger} L_{\alpha, \lambda} \rangle_{S,t} \langle A_{\alpha, \lambda} A_{\alpha, \lambda}^{\dagger} \rangle_{B,t} - \langle L_{\alpha, \lambda} L_{\alpha, \lambda}^{\dagger} \rangle_{S,t} \langle A_{\alpha, \lambda}^{\dagger} A_{\alpha, \lambda} \rangle_{B,t}) \\ &+ (\langle L_{\lambda}^{\dagger} [H_{ND}, L_{\lambda}] + [L_{\lambda}^{\dagger}, H_{ND}] L_{\lambda} \rangle_{S,t}) \langle A_{\alpha, \lambda} A_{\alpha, \lambda}^{\dagger} \rangle_B + (\langle L_{\lambda} [H_{ND}, L_{\lambda}^{\dagger}] + [L_{\lambda}, H_{ND}] L_{\lambda}^{\dagger} \rangle_{S,t}) \langle A_{\alpha, \lambda}^{\dagger} A_{\alpha, \lambda} \rangle_B] + O(\tau^4). \end{aligned} \quad (B23)$$

The same result can be obtained by invoking energy conservation, Eqs. (B7)–(B9).
Comparing Eqs. (B21) and (B22) we see that the amount of energy exchanged by the baths and the system is a linear combination of the $\omega_{\alpha,\lambda}$, i.e., the difference between the diagonal elements of the system Hamiltonian [see Eq. (B5)]. The terms multiplying the energy differences $\omega_{\alpha,\lambda}$ in the first line of Eq. (B21) or Eq. (B22) play the role of kinetic terms expressing the rate of such an energy exchange. The nondiagonal part of the system Hamiltonian does not contribute to heat currents to the baths but is associated with an energy flow between the system and the mechanism connecting the system to the baths, as represented by the Hamiltonian V [see Eqs. (B22) and (B23)].

As in the main text we define the transition rate for a transition $j \to j'$ and its reverse:

$$ \gamma_{\alpha}(j \to j') = \tau g_{\alpha,\lambda}^2 \langle A_{\alpha,\lambda} A_{\alpha,\lambda}^{\dagger} \rangle_B, \quad (B24) $$

$$ \gamma_{\alpha}(j' \to j) = \tau g_{\alpha,\lambda}^2 \langle A_{\alpha,\lambda}^{\dagger} A_{\alpha,\lambda} \rangle_B. \quad (B25) $$

Taking the limit $\tau \to 0$ in Eq. (B22), with $\tau g_{\alpha,\lambda}^2 = \text{const}$, we can distinguish the two heat currents on the RHS of Eq. (B22):

$$ \dot{Q}_{D,\alpha} = \sum_{\lambda} \omega_{\lambda}(\gamma_{\alpha}(j \to j')\langle\pi_{j'}\rangle_{S,t} - \gamma_{\alpha}(j' \to j)\langle\pi_{j'}\rangle_{S,t}), \quad (B26) $$

$$ \dot{Q}_{ND,\alpha} = \frac{1}{2} \sum_{\lambda} [\gamma_{\alpha}(j \to j') (\langle L_{\lambda}^{\dagger}[H_{ND}, L_{\lambda}] + [L_{\lambda}^{\dagger}, H_{ND}]L_{\lambda} \rangle_{S,t}) \\ + \gamma_{\alpha}(j' \to j) (\langle L_{\lambda}[H_{ND}, L_{\lambda}^{\dagger}] + [L_{\lambda}, H_{ND}]L_{\lambda}^{\dagger} \rangle_{S,t})], \quad (B27) $$

where $\pi_j = L_{\lambda}^{\dagger}L_{\lambda} = |j\rangle\langle j|$ as in the main text. Thus, from Eq. (B22) finally we obtain

$$ \begin{aligned} \dot{E}_S &= \sum_{\alpha} \dot{Q}_{D,\alpha} + \dot{Q}_{ND,\alpha} \\ &= \sum_{\alpha} \text{Tr}(\rho D_{\alpha}^*[H_D] + \rho D_{\alpha}^*[H_{ND}]), \end{aligned} \quad (B28) $$

which corresponds to the findings of Ref. [41], and it is the result mentioned in the main text.

Inspection of Eqs. (B21)–(B23) can finally help us to address the following question: which are the energy currents that flow from and/or into the system and the baths? If we assume the energy gaps in the baths and in the system to be resonant, i.e.,

$$ \Omega_{\alpha,\lambda} = \omega_{\lambda}, \quad (B29) $$

a part of the energy currents precisely matches the energy current of the baths: such an energy current is the diagonal heat current (B26) and $\dot{Q}_{D,\alpha} = -\dot{E}_{B,\alpha}$ holds. Yet, even if Eq. (B29) holds, a part of the energy flows toward the interaction mechanism, represented by the Hamiltonian V, and we find $\dot{Q}_{ND,\alpha} = -\dot{E}_{V,\alpha}$. Only when the chosen basis $\{|j\rangle\}$ is the eigenbasis of H and Eq. (B29) holds, the latter contribution vanishes, and all the energy current flows from the system into the baths. Furthermore, from Eq. (B22) we also obtain $\dot{E}_S = \sum_{\alpha} \dot{Q}_{D,\alpha} + \dot{Q}_{ND,\alpha}$. The identical result is obtained by considering $\dot{E}_S = \langle H \rangle$ and the dynamics in Eq. (1). Now, in the steady state $\dot{E}_S = 0$, and thus the RHS of Eq. (B22) vanishes.

Thus by comparing Eqs. (B21) and (B23) we conclude that when the system reaches the steady state the equality $\dot{E}_{B,\alpha} = -\dot{E}_{V,\alpha}$ holds too.

We notice a difference with the notation of the main text: the interaction Hamiltonian V as given by Eq. (25) entails both the forward $j \to j'$ jump ($L_{\lambda}$) and its inverse $j' \to j$ ($L_{\lambda}^{\dagger}$). Thus the sum over $\lambda$ in Eqs. (B26) entails both the transitions, while in the main text the sum over $\lambda$ in, e.g., Eq. (2) involves only the transition $j \to j'$.

It is also worth noticing that, starting from Eq. (B2) and using the Born-Markov approximation discussed above, after calculating the difference $\rho(t+\tau) - \rho(t)$ up to the second order in $\tau$, one recovers the local master equations (1) and (2) of the main text, with dissipation rates given by Eqs. (52) and (51).

## APPENDIX C: ADDITIONAL INFORMATION ON THE NUMERICAL RESULTS

For each trajectory, the QMC algorithm evolves the system's state from an initial pure state $|\psi(0)\rangle$ to a final state $|\psi(t)\rangle$. To initiate a QMC trajectory one thus needs to choose an initial state $|\psi(0)\rangle$ compatible with the initial density matrix $\rho^{(0)}$. One special choice is to take $|\psi(0)\rangle$ to be one of the $\{|k_0\rangle\}$, i.e., the eigenstates of $\rho^{(0)}$, with probability $p_{k_0}$, such that

$$ \rho^{(0)} = \sum_{k_0} p_{k_0} |k_0\rangle\langle k_0|. \quad (C1) $$

This choice corresponds to a noninvasive measurement that does not induce any back-action on the system. This is the approach that we followed for the simulations whose results are shown in Fig. 2 in the main text, and in Fig. 3 in this Appendix. Given that the FT, Eq. (16), is independent of the initial basis, we also consider below the case where the system is initially projected on the basis $\{|j\rangle\}$ defining the jump operators introduced in Eq. (2). In contrast to the previous choice, in this case the measurement is indeed invasive and induces measurement back-action. Nonetheless, our fluctuation theorem is still valid.

We now discuss the evaluation of the system's entropy change with the QMC. We write the quantum FT, Eq. (16), in the form of its classical counterpart, Eq. (3):

$$ \langle e^{-S - \Delta S_S} \rangle_{\text{tr}} = 1, \quad (C2) $$

where now $\langle \cdots \rangle_{\text{tr}}$ is the average over the QMC trajectories.

Given that the generalized entropy $S$ along a trajectory is well defined [see Eq. (22) in the main text], we are left with the question of how to sample the system's entropy change $\Delta S_S$ along the quantum trajectories. Inspired by what one does in the classical case one can possibly sample two different quantities along a single trajectory:

$$ \begin{aligned} \Delta S_{S,x} = & -\langle \psi(t) | \log \rho(t) | \psi(t) \rangle \\ & +\langle \psi(0) | \log \rho^{(0)} | \psi(0) \rangle, \end{aligned} \quad (C3) $$

$$ \begin{aligned} \Delta S_{S,y} = & -\log(\langle \psi(t) | \rho(t) | \psi(t) \rangle) \\ & + \log(\langle \psi(0) | \rho^{(0)} | \psi(0) \rangle), \end{aligned} \quad (C4) $$
<figure><img src="image_5.png" /><figcaption>FIG. 3. QMC simulations of the two-spin system (C7): $\langle \exp(-\Delta S_{tot}) \rangle_{tr}$ and $\langle \exp(-\Delta S_B) \rangle_{tr}$ as functions of $t$, with $\Delta S_{tot} = S + \Delta S_S$. For the system's entropy variation we use the two possible definitions in Eq. (C3) and (C4). Averages over $10^6$ trajectories, $g = 0.1$, $T_a = 1$. (a) Diagonal $\rho^{(0)}$, $J = 0$, $h = 0.2$, $T_b = 1$ (classical case at equilibrium). (b) Diagonal $\rho^{(0)}$, $J = 0.1$, $h = 0.2$, $T_b = 1.2$. (c) Nondiagonal $\rho^{(0)}$, $J = 0.1$, $h = 0.2$, $T_b = 1.2$. (d) Nondiagonal $\rho^{(0)}$, $J = 0.2$, $T_b = 1.2$, and time-dependent field $h(t) = 0.4t/t_f$, with $t_f = 15$. See Appendix C for the detailed description of $\rho^{(0)}$. In panels (a)–(d) the jump rates (C10) are used, which satisfy the LDBC, and thus $S = \Delta S_B$. (e) Nonthermal jump rates $\gamma_{\alpha,\lambda}$, as reported in Table I, that do not satisfy the LDBC, with diagonal $\rho^{(0)}$, and $J = 0.1$, $h = 0.2$.</figcaption></figure>

where $\rho(t)$ is the solution of ME (1). In the classical case, these two quantities are equivalent, as $\rho(t)$ is diagonal, and the state $|\psi(t)\rangle$ is one of the states of the chosen basis. However, these two possible definitions for the system entropy provide a completely different result for a quantum nondiagonal state. This is also evident in our numerical demonstration of the FT, and we anticipate that the second definition is the correct one. We can also provide a theoretical argument for this result. With the QMC, one can evaluate the solution $\rho(t)$ to the ME as

$$ \rho(t) = \langle |\psi(t)\rangle \langle \psi(t)| \rangle_{tr}. \qquad (C5) $$

Strictly speaking, the equality is exact when one considers the whole ensemble of possible trajectories. In a QMC trajectory, from $|\psi(0)\rangle = |k_0\rangle$ to $|\psi(t)\rangle$, one can sample $S$, Eq. (22). Let us now consider the FT in the form of Eq. (16) and address the question of how we can evaluate $\Psi^{(1)}(t|\pi_{k_0})$ with an average over the MC trajectories. Inspection of Eq. (C5), and of the definition of $\Psi(\xi, t)$, Eq. (7), suggests

$$ \Psi^{(1)}(t|\pi_{k_0}) = \langle e^{-S} |_{|\psi(t)|, k_0} | \psi(t) \rangle \langle \psi(t) | \rangle_{tr}. \qquad (C6) $$

Comparison of the last equality with the FT in the form of Eqs. (16) and (C2) suggests that $\Delta S(t)_{S,y}$ as given by Eq. (C4) must be used in the numerical evaluation of the FT.

This conclusion is confirmed by the results for the two-spin system shown in Fig. 3 where the trajectory average in Eq. (C2) is shown for different values of the system parameters and of the initial state $\rho^{(0)}$, and the two possible definitions of system entropy change, Eqs. (C3) and (C4), are used. In particular, in Fig. 3(d) we consider the case of a time-dependent Hamiltonian, where the external field is changed according to $h(t) = h_0 + (h_1 - h_0)t/t_f$. While Figs. 3(a)–3(d) correspond to the case where the baths are thermal, with dissipation rates obeying the LDBC, in Fig. 3(e) we consider the case where the dissipation rates violate the detailed balance (see below for further details).

We can also check that the FT holds for an arbitrary initial basis, as predicted by Eq. (16). In Fig. 4 we start the simulations from (i) the basis $\{|k_0\rangle\}$ that diagonalizes the initial state $\rho^{(0)}$, and (ii) the basis $\{|j\rangle\}$ defining the jump operators introduced in Eq. (2). We see that the FT is numerically satisfied for both choices.

<figure><img src="image_6.png" /><figcaption>FIG. 4. QMC simulations of the two-spin system (C7): $\langle \exp(-\Delta S_{tot}) \rangle_{tr}$ and $\langle \exp(-\Delta S_B) \rangle_{tr}$ as functions of $t$, with $\Delta S_{tot} = \Delta S_B + \Delta S_S$. Nondiagonal $\rho^{(0)}$, $J = 0.1$, $h = 0.2$, $T_b = 1.2$. See Appendix C for the detailed description of $\rho^{(0)}$. We start the simulations with (i) the basis $\{|k_0\rangle\}$ that diagonalizes the initial state $\rho^{(0)}$, and (ii) the basis $\{|j\rangle\}$ defining the jump operators introduced in Eq. (2). The dissipation rates (C10) have been used.</figcaption></figure>
<caption>TABLE I. Dissipation rates for the different transitions without the LDBC. We enumerate the standard basis as follows: $\{|0\rangle, |1\rangle, |2\rangle, |3\rangle\} = \{|\uparrow\rangle, |\uparrow\downarrow\rangle, |\downarrow\rangle, |\downarrow\downarrow\rangle\}$.</caption>

$$
\begin{align*}
\gamma_a (0 \to 2) &= 0.0813 \\
\gamma_a (2 \to 0) &= 0.1335 \\
\gamma_a (1 \to 3) &= 0.0976 \\
\gamma_a (3 \to 1) &= 0.1407 \\
\gamma_a (j \to j) &= 0.1 \\
\gamma_b (0 \to 1) &= 0.1011 \\
\gamma_b (1 \to 0) &= 0.1270 \\
\gamma_b (2 \to 3) &= 0.1112 \\
\gamma_b (3 \to 2) &= 0.1411 \\
\gamma_b (j \to j) &= 0.12
\end{align*}
$$

Other information on the numeric results

We remind the reader that in the main text we consider a
system of two spin-1/2 particles with Hamiltonian

$$
H = -J \sigma_{x,a} \otimes \sigma_{x,b} - h(\sigma_{z,a} \otimes \mathbb{I}_{2,b} \otimes \mathbb{I}_{2,a} \otimes \sigma_{z,b}), \quad (C7)
$$

where each spin is connected to an equilibrium reservoir
at temperatures $T_a$ and $T_b$, respectively. The jump operators
“flip” the individual spins: $L_\lambda = \sigma_{-,a} \otimes \mathbb{I}_b$ or $L_\lambda = \mathbb{I}_a \otimes \sigma_{-,b}$.

The diagonal initial state for Figs. 2(a), 2(b), and 2(d) in the main text and Figs. 3(a), 3(b), and 3(d) in this Appendix reads $\rho_D^{(0)} = (0.4, 0.275, 0.175, 0.15)$ in the basis $\{|\uparrow\uparrow\rangle, |\uparrow\downarrow\rangle, |\downarrow\uparrow\rangle, |\downarrow\downarrow\rangle\}$.

The nondiagonal initial state in Figs. 2(c) and 3(c) was
obtained as follows. We introduced the operator

$$
H_0 = -\frac{1}{2}\pi_{|\uparrow\uparrow\rangle} - \frac{1}{4}\pi_{|\downarrow\downarrow\rangle} - \frac{1}{3}\pi_{|\uparrow\downarrow\rangle} \\ -h\sigma_{x,a} \otimes \mathbb{I}_2 + h\mathbb{I}_2 \otimes \sigma_{x,b}, \qquad (\text{C8})
$$

then

$$
\rho^{(0)} = e^{-H_0/2} / \text{Tr}[e^{-H_0/2}]. \quad (C9)
$$

For the case where the LDBC holds, the bosonic bath
dissipation rates used to solve numerically the quantum ME
(1) and in the QMC algorithm read [55]

$$
\gamma_{\alpha, \lambda} = \frac{g|\omega_{\lambda}|}{1 - e^{\beta_{\alpha}|\omega_{\lambda}|}} \begin{cases} e^{-\beta_{\alpha} \omega_{\lambda}}, & \omega_{\lambda} \ge 0 \\ 1, & \omega_{\lambda} \le 0, \end{cases} \quad (\text{C10})
$$

where g is a microscopic frequency and $\omega_\lambda$ is defined in Eq. (B5). When the LDBC holds, the generalized entropy $S$, Eq. (22), takes the simple form

$$
\begin{align}
S = \Delta S_B &= -(\beta_a Q_{a,D} + \beta_b Q_{b,D}) \nonumber \\
&= -\sum_{l=1}^{n_a} \beta_a \omega_{j_{l+1}, j_l}(t_l) - \sum_{m=1}^{n_b} \beta_b \omega_{j_{m+1}, j_m}(t_m), \tag{C11}
\end{align}
$$

where $\Delta S_B$ is the entropy change in the bath.

For the case where the LDBC does not hold [Fig. 3(e) in
this Appendix and Fig. 2(e) in the main text], we choose the
dissipation rates as in Table I.

APPENDIX D: COMPLEMENT TO SECTION III

This Appendix complements the derivation of the results
in Sec. III.

We first prove Eq. (43). We start by introducing the non-
Hermitian operator

$$
V'_{\alpha, \lambda} = g_{\alpha, \lambda} (L_{\lambda}^{\dagger} A_{\alpha, \lambda} - L_{\lambda} A_{\alpha, \lambda}^{\dagger}), \quad (\text{D1})
$$

and notice that

$$
[H_B, \mathcal{H}] = [H_B, V] = \sum_{\alpha, \lambda} \Omega_{\alpha, \lambda} V'_{\alpha, \lambda}, \quad (\text{D2})
$$

$$
[H_B, [H_B, \mathcal{H}]] = \left[ H_B, \sum_{\alpha, \lambda} \Omega_{\alpha, \lambda} V'_{\alpha, \lambda} \right] = \sum_{\alpha, \lambda} \Omega_{\alpha, \lambda}^2 V_{\alpha, \lambda}. \quad (\text{D3})
$$

We proceed by noticing that for two noncommuting opera-
tors X and Y, the following equality holds:

$$
e^X e^Y = e^{Y+[X,Y]} + \frac{1}{2!}[X,[X,Y]] + \frac{1}{3!}[X,[X,[X,Y]]] + \dots e^X. \quad (\text{D4})
$$

We now turn our attention back to Eq. (35), and taking $X = -S_B/2$ and $Y = \mathcal{H}$ in the previous equation, we find

$$
\begin{align}
e^{-\frac{S_B}{2}} e^{-i\mathcal{H}} &= e^{-i\mathcal{H}} [H+H_B+\sum_{\alpha,\lambda} \sum_{n=0}^{\infty} g_{\alpha,\lambda} (y_{\alpha,\lambda}^{2n} \frac{V_{\alpha,\lambda}}{2n!} - y_{\alpha,\lambda}^{2n+1} \frac{V'_{\alpha,\lambda}}{(2n+1)!})] e^{-\beta_{\alpha} H_B/2} \nonumber \\
&= e^{-i\mathcal{H}} [H+H_B+\sum_{\alpha,\lambda} \bar{V}_{\alpha,\lambda}] e^{-\frac{S_B}{2}}, \tag{D5}
\end{align}
$$

where $y_{\alpha, \lambda} = \beta_{\alpha} \Omega_{\alpha, \lambda} / 2$ and

$$
\bar{V}_{\alpha, \lambda} = g_{\alpha, \lambda} (e^{-\beta_{\alpha} \Omega_{\alpha, \lambda} / 2} L_{\lambda}^{\dagger} A_{\alpha, \lambda} + e^{\beta_{\alpha} \Omega_{\alpha, \lambda} / 2} L_{\lambda} A_{\alpha, \lambda}^{\dagger}) \quad (\text{D6})
$$

is a modified, non-Hermitian interaction Hamiltonian that en-
ters in Eq. (45).

We now turn our attention to Eq. (50), and show how it is
equivalent to Eq. (54) and thus to Eq. (11) in the limit $\tau \to 0$. By using the definition for the rates $\gamma_a (j \to j') = \gamma_{a, j'j}$ introduced in Eqs. (51) and (52), Eq. (50) becomes

$$
\begin{align}
& \frac{\Psi(t + \tau) - \Psi(t)}{\tau} \nonumber \\
& \approx -i[H, \Psi(t)] + \sum_{\alpha, \lambda} \gamma_{\alpha, j'j} e^{\beta_{\alpha} \Omega_{\alpha, \lambda} L_{\lambda}} \Psi(t) L_{\lambda}^{\dagger} \nonumber \\
& \qquad + \sum_{\alpha, \lambda} \gamma_{\alpha, j'j} e^{-\beta_{\alpha} \Omega_{\alpha, \lambda} L_{\lambda}^{\dagger}} \Psi(t) L_{\lambda} \nonumber \\
& \qquad - \frac{1}{2} \sum_{\alpha, \lambda} \gamma_{\alpha, j'j} \{L_{\lambda}^{\dagger} L_{\lambda}, \Psi(t)\} + \gamma_{\alpha, j'j} \{L_{\lambda} L_{\lambda}^{\dagger}, \Psi(t)\} \nonumber \\
& = -i[H, \Psi(t)] \nonumber \\
& \qquad + \sum_{\alpha, \lambda} \gamma_{\alpha, j'j} L_{\lambda} \Psi(t) L_{\lambda}^{\dagger} + \gamma_{\alpha, j'j} L_{\lambda}^{\dagger} \Psi(t) L_{\lambda} \nonumber \\
& \qquad - \frac{1}{2} \sum_{\alpha, \lambda} \gamma_{\alpha, j'j} \{L_{\lambda}^{\dagger} L_{\lambda}, \Psi(t)\} + \gamma_{\alpha, j'j} \{L_{\lambda} L_{\lambda}^{\dagger}, \Psi(t)\}, \tag{D7}
\end{align}
$$

where in the last equality we have used Eq. (53). A careful inspection shows that the diagonal part of the RHS of Eq. (53) corresponds to Eq. (A10), while the nondiagonal part of Eq. (53) corresponds to Eq. (A11). This proves that the nonunitary dynamics introduced in Eq. (47) leads to the modified ME (11) under the same approximations that lead to the GKSL ME (1).
[1] C. Jarzynski, Nonequilibrium Equality for Free Energy Differences, Phys. Rev. Lett. 78, 2690 (1997).

[2] J. Kurchan, Fluctuation theorem for stochastic dynamics, J. Phys. A: Math. Gen. 31, 3719 (1998).

[3] G. E. Crooks, Entropy production fluctuation theorem and the nonequilibrium work relation for free energy differences, Phys. Rev. E 60, 2721 (1999).

[4] G. E. Crooks, Path-ensemble averages in systems driven far from equilibrium, Phys. Rev. E 61, 2361 (2000).

[5] U. Seifert, Entropy Production Along a Stochastic Trajectory and an Integral Fluctuation Theorem, Phys. Rev. Lett. 95, 040602 (2005).

[6] M. Esposito and C. Van den Broeck, Three Detailed Fluctuation Theorems, Phys. Rev. Lett. 104, 090601 (2010).

[7] U. Seifert, Stochastic thermodynamics, fluctuation theorems and molecular machines, Rep. Prog. Phys. 75, 126001 (2012).

[8] A. Imparato and L. Peliti, The distribution function of entropy flow in stochastic systems, J. Stat. Mech.: Theory Exp. (2007) L02001.

[9] H. C. Fogedby and A. Imparato, Heat flow in chains driven by thermal noise, J. Stat. Mech.: Theory Exp. (2012) P04005.

[10] H. C. Fogedby and A. Imparato, Heat fluctuations and fluctuation theorems in the case of multiple reservoirs, J. Stat. Mech.: Theory Exp. (2014) P11011.

[11] M. Esposito, U. Harbola, and S. Mukamel, Nonequilibrium fluctuations, fluctuation theorems, and counting statistics in quantum systems, Rev. Mod. Phys. 81, 1665 (2009).

[12] M. Campisi, P. Hänggi, and P. Talkner, Colloquium: Quantum fluctuation relations: Foundations and applications, Rev. Mod. Phys. 83, 771 (2011).

[13] P. Hänggi and P. Talkner, The other QFT, Nat. Phys. 11, 108 (2015).

[14] J. Dereziński, W. De Roeck, and C. Maes, Fluctuations of quantum currents and unravelings of master equations, J. Stat. Phys. 131, 341 (2008).

[15] F. Liu and J. Xi, Characteristic functions based on a quantum jump trajectory, Phys. Rev. E 94, 062133 (2016).

[16] J. P. Garrahan and I. Lesanovsky, Thermodynamics of Quantum Jump Trajectories, Phys. Rev. Lett. 104, 160601 (2010).

[17] J. M. Horowitz, Quantum-trajectory approach to the stochastic thermodynamics of a forced harmonic oscillator, Phys. Rev. E 85, 031110 (2012).

[18] Y. Subaşı and B. L. Hu, Quantum and classical fluctuation theorems from a decoherent histories, open-system analysis, Phys. Rev. E 85, 011112 (2012).

[19] J. M. Horowitz and J. M. R. Parrondo, Entropy production along nonequilibrium quantum jump trajectories, New J. Phys. 15, 085028 (2013).

[20] B. Leggio, A. Napoli, A. Messina, and H.-P. Breuer, Entropy production and information fluctuations along quantum trajectories, Phys. Rev. A 88, 042111 (2013).

[21] G. Manzano, J. M. Horowitz, and J. M. R. Parrondo, Nonequilibrium potential and fluctuation theorems for quantum maps, Phys. Rev. E 92, 032129 (2015).

[22] G. Manzano, J. M. Horowitz, and J. M. R. Parrondo, Quantum Fluctuation Theorems for Arbitrary Environments: Adiabatic and Nonadiabatic Entropy Production, Phys. Rev. X 8, 031037 (2018).

[23] C. Elouard, D. A. Herrera-Martí, M. Clusel, and A. Auffèves, The role of quantum measurement in stochastic thermodynamics, npj Quantum Inf. 3, 9 (2017).

[24] C. Elouard, N. K. Bernardes, A. R. R. Carvalho, M. F. Santos, and A. Auffèves, Probing quantum fluctuation theorems in engineered reservoirs, New J. Phys. 19, 103011 (2017).

[25] A. Auffèves, A short story of quantum and information thermodynamics, SciPost Phys. Lect. Notes 27 (2021).

[26] H. Kwon and M. S. Kim, Fluctuation Theorems for a Quantum Channel, Phys. Rev. X 9, 031029 (2019).

[27] R. Chetrite and K. Mallick, Quantum fluctuation relations for the Lindblad master equation, J. Stat. Phys. 148, 480 (2012).

[28] G. Manzano, R. Fazio, and E. Roldán, Quantum Martingale Theory and Entropy Production, Phys. Rev. Lett. 122, 220602 (2019).

[29] K. Micadei, G. T. Landi, and E. Lutz, Quantum Fluctuation Theorems Beyond Two-Point Measurements, Phys. Rev. Lett. 124, 090602 (2020).

[30] G. Manzano, D. Subero, O. Maillet, R. Fazio, J. P. Pekola, and E. Roldán, Thermodynamics of Gambling Demons, Phys. Rev. Lett. 126, 080603 (2021).

[31] M. Esposito and S. Mukamel, Fluctuation theorems for quantum master equations, Phys. Rev. E 73, 046129 (2006).

[32] J. Åberg, Fully Quantum Fluctuation Theorems, Phys. Rev. X 8, 011019 (2018).

[33] K. Micadei, J. P. Peterson, A. M. Souza, R. S. Sarthour, I. S. Oliveira, G. T. Landi, R. M. Serra, and E. Lutz, Experimental validation of fully quantum fluctuation theorems, Phys. Rev. Lett. 127, 180603 (2021).

[34] G. E. Crooks, On the Jarzynski relation for dissipative quantum dynamics, J. Stat. Mech.: Theory Exp. (2008) P10023.

[35] M. Campisi, P. Talkner, and P. Hänggi, Fluctuation Theorem for Arbitrary Open Quantum Systems, Phys. Rev. Lett. 102, 210401 (2009).

[36] T. Albash, D. A. Lidar, M. Marvian, and P. Zanardi, Fluctuation theorems for quantum processes, Phys. Rev. E 88, 032146 (2013).

[37] M. Silaev, T. T. Heikkilä, and P. Virtanen, Lindblad-equation approach for the full counting statistics of work and heat in driven quantum systems, Phys. Rev. E 90, 022103 (2014).

[38] E. Davies, Markovian master equation, Commun. Math. Phys. 39, 91 (1974).

[39] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Completely positive dynamical semigroups of N-level systems, J. Math. Phys. 17, 821 (1976).

[40] G. Lindblad, On the generators of quantum dynamical semigroups, Commun. Math. Phys. 48, 119 (1976).

[41] A. Hewgill, G. De Chiara, and A. Imparato, Quantum thermodynamically consistent local master equations, Phys. Rev. Research 3, 013165 (2021).

[42] J. L. Lebowitz and H. Spohn, A Gallavotti–Cohen-type symmetry in the large deviation functional for stochastic dynamics, J. Stat. Phys. 95, 333 (1999).

[43] M. Esposito, U. Harbola, and S. Mukamel, Entropy fluctuation theorems in driven open systems: Application to electron counting statistics, Phys. Rev. E 76, 031132 (2007).

[44] R. Chetrite and H. Touchette, Nonequilibrium Markov processes conditioned on large deviations, Ann. Henri Poincaré 16, 2005 (2015).
[45] F. Carollo, J. P. Garrahan, I. Lesanovsky, and C. Pérez-Espigares, Making rare events typical in Markovian open quantum systems, Phys. Rev. A 98, 010103(R) (2018).

[46] R. Alicki, The quantum open system as a model of the heat engine, J. Phys. A: Math. Gen. 12, L103 (1979).

[47] K. Mølmer, Y. Castin, and J. Dalibard, Monte Carlo wave-function method in quantum optics, J. Opt. Soc. Am. B 10, 524 (1993).

[48] J. Johansson, P. Nation, and F. Nori, QUTIP: An open-source PYTHON framework for the dynamics of open quantum systems, Comput. Phys. Commun. 183, 1760 (2012).

[49] J. Johansson, P. Nation, and F. Nori, QUTIP 2: A PYTHON framework for the dynamics of open quantum systems, Comput. Phys. Commun. 184, 1234 (2013).

[50] J. J. Sakurai and J. Napolitano, Modern Quantum Mechanics, 2nd ed. (Cambridge University Press, Cambridge, U.K., 2017).

[51] C. Maes, Local detailed balance, SciPost Phys. Lect. Notes 32 (2021).

[52] G. E. Crooks, Quantum operation time reversal, Phys. Rev. A 77, 034101 (2008).

[53] C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg, Atom-Photon Interactions: Basic Processes and Applications (Wiley, New York, 1992).

[54] Given the assumption that the bath Hamiltonian commutes with the corresponding initial state this is equivalent to assume that $A_{\alpha,\lambda}$ are eigenoperators of $H_{B\alpha}$.

[55] H. P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University Press, Oxford, U.K., 2002).

[56] S. Haroche, M. Brune, and J. Raimond, From cavity to circuit quantum electrodynamics, Nat. Phys. 16, 243 (2020).

[57] P. Schindler, M. Müller, D. Nigg, J. T. Barreiro, E. A. Martinez, M. Hennrich, T. Monz, S. Diehl, P. Zoller, and R. Blatt, Quantum simulation of dynamical maps with trapped ions, Nat. Phys. 9, 361 (2013).

[58] F. V. Melo, N. Sá, I. Roditi, G. T. Landi, A. M. Souza, I. S. Oliveira, and R. S. Sarthour, Experimental implementation of a two-stroke quantum heat engine, arXiv:2203.13773.

[59] G. De Chiara and A. Imparato, Dataset for “A quantum fluctuation theorem for dissipative processes”, Queen’s University Belfast, 2022, doi: 10.17034/b3e856e6-87dc-4849-8c0b-5785856cb618.