A quantum fluctuation theorem for any Lindblad master equation

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.


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 theorems (FTs) extending the second law of thermodynamics to the microscopic realm [1][2][3][4][5][6][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 interacting 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 generalise the FT and the impossibility of monitoring a quantum system without disturbance has generated a long debate (see, for example, the reviews [11][12][13]). Several approaches to the quantum FT employ the formalism of quantum trajectories or quantum Monte Carlo (QMC) [14][15][16][17][18][19][20][21][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. Interestingly, the back-action from the observation of the jump, e.g. through the emission of a photon, gives rise to a genuinely quantum energy contribution, dubbed "quantum heat" [23][24][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][28][29][30]. One can also prove the FT by diagonalizing the instantaneous density matrix, thus gen-erating 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 quantum 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 accounts 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 interaction 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 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 timedependent 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][35][36][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][39][40] ( = 1): for its density matrix ρ(t) with dissipators where α labels the environment and L λ = |j j| and λ = λ(j → j ) denotes a transition between two states |j and |j . 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 }, though arbitrary, is motivated by the physical environments. If the states |j are the eigenstates of H the ME is dubbed "global" or "local" otherwise [41]. We split an operator X = X D + X N D into the sum of its diagonal and nondiagonal parts in the basis {|j } and set X jj = j|X|j . We consider the general case where different baths can drive the same transition λ, if the dissipation rates (potentially time dependent) γ α,λ = 0. In the classical case the FT reads [7] e −∆S B −∆S S = 1, where ∆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 α reads ∆S B = − α β α Q α . 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 α > 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 γ j j because of the interaction with the baths: we thus introduce the elementary current associated with the jump λ ∆s α,j j = − log(γ α,jj /γ α,j j ).
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 ∆S B in Eq. (3): we will elaborate later on this connection. We would thus like to characterize the joint probability distribution Φ j (S, t) of finding the system in the state |j 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 ρ(S, t), such that its diagonal elements describe the desired joint probability ρ jj = Φ j (S, t). It is possible to show that the modified density matrix satisfies the mod-ified ME: has vanishing diagonal terms and only couples nondiagonal terms of the density matrix in the chosen basis. A derivation of Eqs. (5)-(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 → j such a current reads S + ∆s α,j j . Next, we introduce the generating function operator Ψ(ξ, t) akin to the generating function in classical physics, defined as and applying this integral transform to both sides of Eqs. (5)-(6), one obtains an equation for Ψ(ξ, t) where the superoperator L ξ [·] depends parametrically on ξ, and its full expression is given in Appendix A, see Eqs. (A8)-(A9). We will see that the modified density matrix ρ(t), its generating function Ψ(t) and their modified MEs are the essential ingredients to prove the quantum FT. We now discuss the physical initial condition for Eqs. (5)-(6) and Eq. (8). Let ρ (0) be the initial state at t = 0. Since no current S has yet been generated: ρ(S, t = 0) = ρ (0) δ(S), where δ(x) is the Dirac-delta function. Given the definition of Ψ, Eq. (7), its initial condition becomes Ψ(ξ, t = 0) = ρ (0) , ∀ξ and in the following Ψ(ξ, t| ρ (0) ) will indicate the solution to Eqs. (8) with this specific initial condition. The following equalities hold e −ξS πj ,ρ (0) = Tr[π j Ψ(ξ, t|ρ (0) )] = Ψ jj (ξ, t|ρ (0) ), (9) e −ξS where · · · πj ,ρ (0) in Eq. (9) denotes the expectation value constrained by the initial state ρ (0) , the state at time t taken to be π j = |j j|, and in Eq. (10) we summed over all possible final states. In the following we will mostly be interested in the case ξ = 1 to evaluate the expectation value of exp(−S). Letting Ψ (1) (t) = Ψ(ξ = 1, t), one finds where D * α is the dual of the dissipator D α : The details of the derivation of Eq. (11) are given in Appendix A. Thus the time evolution of Ψ (1) , as given by Eq. (11), has the same conservative part as in Eq. (1), but the dissipative part (the dual of D α ) 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 Ψ (1) bears information on both the system state ρ(t) and on the quantity exp(−S). As such, the operator Ψ (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 [16] and later in [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 (5)-(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 Ψ (1) (t) Given that Ψ (1) (t) is a linear combination of solutions of Eq. (8) with ξ = 1, it is a solution itself, with initial condition Ψ (1) (0) = j0 π j0 = I. Inspection of Eq. (11), suggests that Ψ (1) (t) is a stationary solution at any time: The last equality for Ψ (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 Ψ (1) (t | π j0 ) 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 }, and noticing that Eq. (8) is linear, we can write the solution at time t, with the specific initial condition π b0 as Furthermore Eqs. (10), (13), and (14) imply 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: where we have used Eqs. (13)- (15) and introduced an arbitrary but normalised final state ρ (f) . This feature of the FT was already noticed in Ref. [7] for the classical case: Eq. (16) holds for any normalised quantum final state ρ (f) , not necessarily the solution of Eq. (1) at time t.
Eq. (16) has the form of the expectation value of the Hermitian operator Ψ (1) over the final state of the system. We then recall the physical interpretation of the operator Ψ (1) (t| ρ (0) ): its jth diagonal element represents the expectation value of exp(−S) constrained by the initial ρ (0) and final state π j , Eq. (9). We make no specific assumption on the coherence in the initial and the final states ρ (0) and ρ (f) : they can both exhibit coherence in the jump-operators' basis {|j } as well as in the energy eigenbasis. Let {|k 0 } and {|r f } be the eigenbases of ρ (0) and ρ (f) , respectively. Changing basis, from {|j } to {|r f } and recalling definition (7), we can rewrite Eq. (16) as 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 Φ r f (S). Moreover, Eqs. (16)-(17) exhibit the structure of a double trace, one over the initial state ρ (0) and one over the final one ρ (f) , and as such their values do not depend on the chosen basis. Notice that performing the initial projection of ρ (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: We can now establish the relation between the dynamics of the system, as expressed by Eqs. (1)-(2) and the thermodynamics significance of our results in Eqs. (14), (16), (17). Such a relation follows immediately if the dissipation rates γ α,λ are taken to obey the local detailed balance condition (LDBC): and β α is the inverse temperature of the bath α. If the LDBC holds then ∆s α,j j in Eq. (4) can be immediately interpreted as the entropy change in the bath, given that −ω j j = β −1 α log(γ α,j j /γ α,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 = ∆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 −ω 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 γ α,j j and their relation to the ω j j need be estimated a priori.
Assuming the LDBC, since the heat flowing into a bath as a consequence of a jump is −ω j j , involving only the diagonal part of the system Hamiltonian, we have S = ∆S B = − α β α Q D,α and we can rewrite Eq. (18) as 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 asQ D,α = Tr {ρD * α [H D ]}, enters the differential version of the second law: d t S S ≥ α β αQD,α 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 [41]. However the nondiagonal part of the Hamiltonian, specifically the quan-tityQ N D,α = Tr {ρD * α [H N D ]}, 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 definitioṅ Q D,α =Q α = Tr {D α [ρ]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 ρ (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 at t = 0 and ending in the state |j . From Eq. (16) or (17) one then recovers immediately the classical FT Our results also hold when the system's Hamiltonian and/or the dissipation rates γ α,γ , entering Eq. (4), are time-dependent. In the special case of a single bath at inverse temperature β and choosing the final and initial states to be the Gibbs states ρ (f) = exp(−βH(t))/Z t , and where ∆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 where we have set ∆S B = −βQ D , and have defined work as W D ≡ ∆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 adjust instantaneously to the value of the energy level, which in turn holds only for slow driving.

A. 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, characterised by Pauli . 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 non thermal, with the dissipation rates not satisfying the LDBC. The jump operators "flip" the individual spins: L λ = σ −,a ⊗ I b or L λ = I a ⊗ σ −,b . When J = 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 |ψ(0) to |ψ(t) 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: 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 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 ρ(t) to the ME (1) as ρ(t) = |ψ(t) ψ(t)| tr , where · · · tr is the average over the QMC trajectories. Analogously, from the definition of Ψ(ξ, t) in Eq. (7), and from the sampling of the generalized entropy (22) along the trajectories, we can evaluate Ψ (1) (t | π k0 ), for any chosen initial basis {|k 0 }, as discussed in Appendix C.
The system's entropy change ∆S S can be also evaluated along the QMC trajectories as detailed in Appendix C, and we can thus evaluate exp (−∆S tot ) tr , with ∆S tot = S + ∆S s .
We consider both the case where the baths are thermal, with rates γ α,λ 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 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, exp(−∆S B ) ≥ 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 (16)-(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 Ψ (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)- (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α : the total Hamiltonian H thus reads where V α,λ is the interaction Hamiltonian with the α-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 In panels (a)-(d) the LDBC is used, and thus S = ∆SB. (e) Non-thermal jump rates γ α,λ that do not satisfy the LDBC, with diagonal ρ (0) , and J = 0.1, h = 0.2. See Appendix C for further details on the numerical simulations.
Furthermore, we introduce the operators akin to an entropy operator for the single bath and for all the baths, respectively. The case where some of eigenvalues of the ρ Bα (0) are zero does not pose an issue on the discussion below, as we will be interested in the exponential of S Bα . We do not assume, a priori, that the initial states of the baths ρ Bα (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 θ, such that [θ, H] = 0 and −iHθ• = θiH•, the latter equality expressing the antilinearity of the operator θ [50].
We assume that the initial state for the total system is a pure state |J 0 , (eigenstate of some observable). Here and in the following |J x indicates states of the total system and Π x = |J x J x | the corresponding projectors, while |j x and π x = |j x j x | are states and projectors of the system of interest only as in the previous section.
The probability to observe the state |J f after the total system evolves until time t is where the unitary operator U t reads and T is the time-ordering operator. We also introduce the unitary operator when the global system evolves under the time reversed Hamiltonian We now assume we can prepare the total system in the state |J f = θ|J f and evaluate the probability to find the system in the state |J 0 = θ|J 0 after evolving under the time reversed dynamics for a time t where we have used the antilinearity of θ. 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, ρ B (0). We are then interested in the probability 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˜ (0) = θπ f ⊗ ρ B (0)θ −1 , and ask what is the probability that the system of interest is in the state |j 0 after a time t. We remind the reader that θ is the time reversal operator for the whole system. A straightforward calculation gives = P (j f , t|j 0 ).
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 A j f ,j0 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 Using the definition of S B in Eq. (29) we consider the two-time correlation 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 β α , Eq. (37) becomes and ∆S B = α β α ∆H Bα is the change of entropy in the baths, along any system trajectory connecting j 0 to j f . From Eq. (37) we obtain immediately We then continue along the lines of the previous section and consider the arbitrary basis {|b 0 } and arbitrary initial and final states ρ (0) and ρ (f) , and from Eq. Tr Eqs.(39)-(40) represent the FT for the joint unitary evolution of the system and the baths: they generalise Eqs. (16)-(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 measurements on the system is required, represented by the projection onto the arbitrary basis state |b 0 . As discussed in the previous section, while the final state ρ (f) is arbitrary, one can of course take the specific choice ρ (f) = ρ(t) = Tr B (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., [52].

A. 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 α-th bath, introduced in Eq. (25), we choose the general form where L λ are system operators, and A α,λ 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α , ρ Bα (0)] = 0. Furthermore, we choose the bath operators to be eigenoperators of the bath entropy operator [54] [S Bα , A α,λ ] = Ω α,λ A α,λ .
It is always possible to express the interaction Hamiltonian (41) in terms of eigenoperators of S Bα . Indeed if A α,λ are not eigenoperators of S Bα , by expressing them in terms of the eigenbasis of S Bα one can write V α in terms of new (rotated) operators A α,λ which are eigenoperators of S Bα . 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 with and where is a modified, non-Hermitian interaction Hamiltonian, see Appendix D.
Let U t be the nonunitary operator obtained from Eq. (31) with the substitution V → V . Then comparing (43) with (37), we finally find resembling Eq. (33) that expresses the conservation of probability under time reversal for the total system. However, in (46) the forward time evolution occurs with U t .
Let us now consider the evolution of the system operator given an arbitrary (possibly mixed) initial state of the system Ψ(t). The projection of such an operator on the state j f , is the two-time correlation function after a time t + τ , as introduced in Eqs. (37) and (46), given the system state at time t Ψ(t) 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 = Ψ(t) ⊗ ρ B , and that the baths' state is time independent ρ B = ρ B (0) [53,55]. Thus in Eq. (47) we take the expansion up to the second order in τ : some extra care is needed as V (see Eq. (45)) and thus H, is not Hermitian. We thus obtain A lengthy but straightforward calculation gives We now introduce the rates We see that, by virtue of Eq. (42), A α,λ A † α,λ and A † α,λ A α,λ are projectors onto the eigenstates of S Bα (28). Furthermore, by taking A α,λ = |N α N α |, we have that the quantity Ω α,λ introduced in Eq.(42) becomes an entropy gap Ω α,λ = −(s N α − s Nα ) = log r N α − log r Nα , where r Nα are the eigenvalues of the bath state ρ Bα . 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 We remind the reader that, given our choice of the in-teraction Hamiltonian (41) and the jump operators L λ and A α,λ , a jump λ in the system j → j corresponds to a jump in the bath N α → N α . Eq. (53) determines the asymmetric part of the system rates, and involves the entropy change in the bath for a transition (j → j , N α → N α ), or equivalently the ratio of the population between the two states N α and N α in the initial bath state. Neither the entropy nor the energy of the system appear in this relation.
Let j = H jj be the diagonal elements of the system Hamiltonian in the basis {|j }. Thus the quantities j are eigenvalues of H only if the chosen basis is an eigenbasis of H. Eq. (53) becomes the standard local detailed balance condition involving the system energy gaps introduced in Sec. II γ α (j → j )/γ α (j → j) = exp(−β α ω λ ) only under two additional conditions: (i) the baths are initially in the thermal state, thus Ω α,λ = −β α (E N α − E Nα ), where E Nα are the eigenvalues of H Bα , and (ii) the energy gaps in the baths and in the system are resonant, i.e., E N α −E Nα = ω λ = 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 α ) 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 } 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 (51)-(52), the detailed balance condition (53), and taking the limit τ → 0, we obtain , i.e. the generalized ME Eq. Thus comparing this last result with Eqs. (47), (48), we conclude that the operator Ψ (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 → j , N α → N α ), the change in the bath entropy is given by Ω α,λ . We recall the definition of Ψ (1) (t) = Ψ(ξ = 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 Eq. (53). If the baths are initially thermal, and conserve their state for the entire duration of the system dynamics, Ω α,λ = −β α (E N α − E Nα ), and thus Eq. (53) implies that the quantity S is the change in the baths' energy multiplied by the corresponding inverse temperature: S = α β α ∆H α . 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]. 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. All data created during this research is openly available from QUB-Pure at https://doi.org/ 10.17034/b3e856e6-87dc-4849-8c0b-5785856cb618 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 ρ(t) is described by the GKSL ME (1). We notice that the diagonal terms of the density matrix ρ(t) express the time evolution of the populations, i.e., the probability p j (t) = ρ jj (t) of observing the system in the state |j at time t. We are interested in the joint probability distribution Φ j (S, t) of finding the system in the state j with a total generalised 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 ρ(S, t), such that its diagonal elements describe the desired joint probability ρ jj = Φ 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 nondiagonal elements. The coupling between the diagonal and nondiagonal elements of ρ occurs through the coherent part of the dynamics, expressed by the commutator in Eq. (1). Thus we can write where we remind the reader that λ = λ(j → j ) denotes a transition between two states |j and |j . Isolating the diagonal part of (A2), one obtains {γ α,j j ρ jj − γ α,jj ρ j j } .
(A3) The jump between the states |j and |j occurs with rate γ α (j → j ) = γ α,j j , and the corresponding entropy change in the bath α is given by ∆s α,j j = − log(γ α,jj /γ α,j j ).
(A4) which corresponds to Eq. (4) in the main text. Such jumps in the ME (1) are described by D D [ρ D ] implicitly introduced in Eq. (A3), thus only the coupling between the diagonal elements in the ME contributes to the change in the bath generalised entropy S. After the jump j → j the generalised entropy changes as S → S + ∆s α,j j . Let us now suppose that we know the modified density matrix ρ(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 + τ , we have For the nondiagonal part of the dynamics, that does not contribute to S, we have Taking the limit τ → 0 in both (A5) and (A6), and expanding the right-hand side of Eq. (A5) in Taylor series of the microscopic entropy change ∆s α,j j , we obtain the modified quantum ME Eqs. (5)-(6) in the main text. We introduce the operator Ψ(ξ, t) defined as and applying this integral transform to both sides of Eqs. (5) and (6) in the main text, one obtains the equations for Ψ(ξ, t) We now take ξ = 1 in Eqs. (A8)-(A9), yielding The dual of the dissipator D α [·] applied on an operator X reads and considering the diagonal element one finds that it corresponds to the term on the righthand side of Eq. (A10). For the nondiagonal part, one All in all, one finds that Ψ (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α , interacting through a Hamiltonian V , see Eqs. (23)- (25). Given an arbitrary orthonormal basis {|j } for the system, we split its Hamiltonian into its diagonal and nondiagonal parts The time evolution of the total density matrix is thus given by for any t and τ . 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 are eigenstates of the system Hamiltonian.
We recall the definition of the interacting Hamiltonians introduced in Eq. (25) V α,λ = g α,λ (L † λ A α,λ + L λ A † α,λ ), with the system operators L λ and the bath operators A α,λ satisfying with Eq. (B3) follows immediately from the choice of the jump operators L λ = |j j|, while in this Appendix we choose the operators A α,λ to be eigenoperators of the α-th Hamiltonian, with A α,λ = |N α N α |, and |N α 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 → j , as long as g α,λ = 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 τ : To calculate the second order term in ∆E B we need to calculate, where the individual commutators read where λ and µ are different transition indices. For future reference, we notice that only the last commutator contains quadratic terms of the types A α,λ A † α,λ or A † α,λ A α,λ . We now turn our attention to the system energy balance. Let us introduce the operators From Eqs. (B8)-(B10), we see that to the first order in τ we need to evaluate To calculate the second order term in (B10) we need the We have We notice that the only commutators that are quadratic in the bath operators A λ and A † λ 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, ρ tot (t) = ρ S (t) ⊗ ρ B (t), that the bath density matrix is time independent: ρ B (t) = N b α=1 ρ Bα , and that [ρ Bα , H Bα ] = 0 even though ρ Bα 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 A similar approach for the interaction Hamiltonian leads to 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 ω α,λ , i.e. the difference between the diagonal elements of the system Hamiltonian, see Eq. (B5). The terms multiplying the energy differences ω α,λ in the first line of Eq. (B21) or (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)-(B23).
As in the main text we define the transition rate for a transition j → j and its reverse Taking the limit τ → 0 in (B22), with τ g 2 α,λ = const., we can distinguish the two heat currents on the RHS of Eq. (B22) where π j = L † λ L λ = |j j| as in the main text. Thus, from Eq. (B22) finally we obtaiṅ 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/into the system and the baths? If we assume the energy gaps in the baths and in the system to be resonant, i.e.
a part of the energy currents precisely matches the energy current of the baths: such an energy current is the diagonal heat current (B26) andQ D,α = −Ė Bα holds. Yet, even if (B29) holds, a part of the energy flows toward the interaction mechanism, represented by the Hamiltonian V , and we findQ N D,α = −Ė Vα . Only when the chosen basis {|j } is the eigenbasis of H and (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Ė s = αQ D,α +Q N D,α . The identical result is obtained by consideringĖ S = Ḣ and the dynamics in Eq. (1). Now, in the steady statė 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 equalitẏ E Bα = −Ė Vα 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 → j jump (L λ ) and its inverse j → j (L † λ ). Thus the sum over λ in Eqs.(B26) entails both the transitions, while in the main text the sum over λ in, e.g., Eq. (2), involves only the transition j → j .
It is also worth noticing that, starting from Eq. (B2) and using the Born-Markov approximation discussed above, after calculating the difference ρ(t + τ ) − ρ(t) up to the second order in τ , one recovers the local master equation (1) 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 } 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): where now · · · tr is the average over the QMC trajectories.
Given that the generalised 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 ∆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 ∆S S,y = − log ( ψ(t)|ρ(t)|ψ(t) ) + log( ψ(0)|ρ (0) |ψ(0) )). (C4) where ρ(t) is the solution of the ME (1). In the classical case, these two quantities are equivalent, as ρ(t) is diagonal, and the state |ψ(t) 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 ρ(t) to the ME as Strictly speaking, the equality is exact when one considers the whole ensemble of possible trajectories. In a QMC trajectory, from |ψ(0) = |k 0 to |ψ(t) , 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 Ψ (1) (t | π k0 ) with an average over the MC trajectories. Inspection of Eq. (C5), and of the definition of Ψ(ξ, t), Eq. (7), suggests Comparison of the last equality with the FT in the form of Eqs. (16) and (C2), suggests that ∆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 twospin 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 ρ (0) , and the two possible definitions of system entropy change (C3)-(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 panels (a)-(d) in Fig. 3 correspond to the case where the baths are thermal, with dissipation rates obeying the LDBC, in panel (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 } that diagonalises the initial state ρ (0) , and (ii) the basis {|j } defining the jump operators introduced in Eq. (2). We see that the FT is numerically satisfied for both choices.

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σ x,a ⊗ σ x,b − h(σ z,a ⊗ I 2,b + ⊗I 2,a ⊗ σ z,b ) (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 λ = σ −,a ⊗ I b or L λ = I a ⊗ σ −,b .
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 τ → 0. By using the definition for the rates γ α (j → j ) = γ α,j j introduced in Eqs. (51)-(52), Eq. (50) becomes 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 (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).