Quantum thermodynamically consistent local master equations

Local master equations are a widespread tool to model open quantum systems, especially in the context of many-body systems. These equations, however, are believed to lead to thermodynamic anomalies and violation of the laws of thermodynamics. In contrast, here we rigorously prove that local master equations are consistent with thermodynamics and its laws without resorting to a microscopic model, as done in previous works. In particular, we consider a quantum system in contact with multiple baths and identify the relevant contributions to the total energy, heat currents and entropy production rate. We show that the second law of thermodynamics holds when one considers the proper expression we derive for the heat currents. We confirm the results for the quantum heat currents by using a heuristic argument that connects the quantum probability currents with the energy currents, using an analogous approach as in classical stochastic thermodynamics. We finally use our results to investigate the thermodynamic properties of a set of quantum rotors operating as thermal devices and show that a suitable design of three rotors can work as an absorption refrigerator or a thermal rectifier. For the machines considered here, we also perform an optimisation of the system parameters using an algorithm of reinforcement learning.

When studying the thermodynamic properties of a system, one has to appropriately model the interaction with the environment [37,38]. A popular approach employs the wellstudied Gorini-Kossakowski-Lindblad-Sudarshan (GKLS) master equation (ME), which describes the interaction of a system with Markovian environment [39][40][41]. In the so-called local approximation to the ME, the corresponding jump operators act only on local subsystems rather than on the eigenstates of the whole Hamiltonian as in the global approach. Many works have analysed the accuracy of local versus global ME [42][43][44][45][46][47][48][49][50][51][52]. Besides, there has been much debate about the thermodynamic consistency of the local ME [53,54], but these issues can be addressed by a careful microscopic modelling of the master equation in question [55][56][57].
In the present paper, we demonstrate the compatibility of local ME with thermodynamics without resorting to microscopic models as done in Refs. [55][56][57]. We achieve this result by identifying, in the time evolution of the energy, one of two terms which alone enters the second law of thermodynamics and can thus be identified with the heat exchanged with the environments. The other term, arising because of the non-compatibility of jump operators and energy eigenstates, gives rise to an additional energy current, that can be identified as the work rate within the framework of a microscopic collisional model. In fact, our results agree with the energy splitting that was previously proposed in Refs. [55][56][57] using a collisional model approach. Finally, our approach allows us to transparently recover the standard definition of heat, and thus of entropy production when considering the corresponding global ME.
We also present a semiclassical argument that relates the energy currents in an open quantum system, whose dynamics is described by a ME, to the probability currents, analogously to what one customarily does in classical stochastic thermodynamics.
We use our findings on the local ME to study the thermodynamic properties of two types of thermal machines, namely absorption refrigerators and thermal rectifiers. We show that one can build efficient thermal devices consisting of quantum rotors, that interact through a clock model Hamiltonian, as working media. The clock model is a generalisation of the spin-1/2 model [58], and has attracted considerable attention in condensed matter physics, with works focusing on the rich phase behaviour of such systems [59][60][61][62][63] and more recently in the context of time crystals [64,65]. Recently it has been shown that the chiral version of the clock model in contact with multiple baths at different temperatures, can convert heat currents into rotational motion. This conversion is the result of the lack of rotational symmetry in the Hamiltonian and of thermal equilibrium, with the device working as an autonomous thermal motor both in the classical [66,67] and in the quantum regime [68].
After introducing the working principles of the proposed thermal machines we will discuss a procedure to optimise their output or efficiency. The performance of the devices is analysed and optimum setups for specific interesting cases are found. In particular we utilise the differential evolution approach [69,70] to find optimum parameter choices. Such a scheme, also called reinforcement learning, has been used, e.g., to find the optimal network topology in interacting elec-tronic systems working as thermoelectric nanoscale engines [71]. This paper is organised as follows: in Sec. II we will introduce the master equation, and discuss the different contributions to the system energy evolution. We derive the second law of thermodynamics building on this analysis. In Sec. III, we present an alternative derivation of the heat currents inspired by classical stochastic thermodynamics. The clock model is reviewed in Sec. IV. We then set the stage of the applicative part by reviewing the properties of a rotor dimer in Sec. V. In Sec. VI we consider a trimer system that works as a refrigerator while in Sec. VII we study a rectifier as a thermal control device. We then conclude this work in Sec. VIII.

II. THE QUANTUM MASTER EQUATION AND ITS ENERGETICS
In this section, we consider the general case of a system with Hamiltonian H, in contact with N b baths, each of which will be denoted with the letter α at the respective inverse temperature β α . The system dynamics is described by a standard GKLS ME for the density matrix ρ [37] ( = 1): with dissipators We choose an arbitrary set of kets | j , and assume that they form an orthogonal basis for the system. We also choose the jump operators L λ to be expressed in terms of such kets, namely, they are of the form L λ = | j j|. Here and in the following λ( j → j ) denotes a transition between two states | j and | j . To lighten the notation we omit the initial and final state of such a transition: so in the following | j always indicates the initial state of an arbitrary transition λ, and | j its final state. In Eq. (2) all the dependencies on the specific bath α are contained in the dissipation rates γ λ,α , which must obey the (local) detailed balance condition: where In this way we consider the general case where different baths can drive the same transition λ, as long as γ λ,α 0. In the following the dependency of γ λ,α on ω λ will be implicitly understood. The Hamiltonian H may or may not be diagonal in the basis | j . The resulting ME is sometimes dubbed global, in the former case, or local, in latter case.
The evolution of an operator A in the Heisenberg picture is given by where D * α [·] is the dual of D α : In particular, if one is interested in studying the system thermodynamics, it is relevant to consider the evolution of the energy operator H. It is convenient to split the Hamiltonian in its diagonal and non-diagonal parts H = H D + H ND , for reasons that will become apparent below. The time evolution of H is thus given by In this work we have implicitly assumed that the system Hamiltonian does not depend explicitly on time. Upon relaxing this condition, one would have to add a term ∂ t H to the right hand side (RHS) of Eq. (7) which would be responsible for an additional work source [11]. After a straightforward manipulation one finds It is worth noting that, within this framework, from Eq. (7) one obtains the standard definition of heat current, which readsQ where we have introduced the expectation valueṡ We see that while the time evolution of the Hamiltonian in the Heisenberg picture, as given by Eq. (7), is general, the subdivision on the RHS of such an equation in two terms which depend on H D and H ND is completely arbitrary and depends on the chosen basis. It is however relevant to note that when such a basis has been chosen, it is the heat current associated with the diagonal part of the Hamiltonian,Q D,α that enters the second law of thermodynamics in terms of the irreversible entropy produc-tionΣ:Σ where we have defined the system entropy S = − tr ρ ln ρ. Eq. (13) is the first relevant result of this paper, and we proceed now with the proof of such an inequality.
The time derivative of the entropy reads where is the total Liouvillian superoperator appearing in Eq. (1). Let us also introduce the partial , where we have used our assumption that there are N b thermal baths at inverse temperature β α . The thermal equilibrium local state ρ α = exp(−β α H D )/Z α is a steady state for the partial superoperator L α . We can now use Spohn's inequality [72], that states that for any superoperator L x of Lindblad form, with steady state ρ ss , (i.e. L x [ρ ss ] = 0), the following inequality holds Let us inspect the second term in this inequality. We have where L * x is the dual of L x . Considering now the specific case of L α , Eq. (15) and (16) give that proves Eq. (13), and thus the second law for the diagonal part of the heat aloneQ D,α . We remark again that, similarly to the RHS of Eq. (7), the detailed expression of the inequality Eq. (13) depends on the arbitrarily chosen basis. We cannot use Eqs. (15) and (16) for the total superoperator L that contains both H D and H ND as the expression of its steady state is in general not an equilibrium state, and furthermore is not known in most of the cases. On the other hand, for the partial superoperator L α we can use in Spohn's inequality the local reference state ρ α , carrying information about the temperature of the bath α.
Spohn inequality was used in ref. [1] to derive an expression of the second law similar to Eq. (13), but withQ α (Eq. (10)) instead ofQ D,α . However, here we have shown that such an expression of the second law is correct only if one deals with a global ME, for which the identityQ α =Q D,α holds. When the Hamiltonian is not diagonal in the basis | j defining the jump operators, Eq. (13) is the correct form of the second law. A similar approach to prove the second law in presence of a single bath was used in [56], but there the ME was explicitly taken to be global.
A posteriori, one should also conclude that, within the framework of the local ME, it makes sense that of the two components appearing on the RHS of Eq. (10), onlyQ D,α enters the second law: the energy exchange with environment is encoded in the ME by the detailed balance condition (3), where only the diagonal components of H in the chosen basis appear.
Let us make a few considerations. We mentioned in the introduction that Ref. [53,54] shows that using Eq. (10) as a definition for the heat currents can lead to apparent nonphysical results, for instance spontaneous heat flow from the cold to an hot bath in a system of two coupled harmonic oscillators. This inconsistency has been resolved in [55][56][57], where a collisional microscopic model was used to realise the interaction of a system with multiple baths. We can now make contact between the quantities defined in this work and the findings of Ref. [57] clarifying even further the origin of each term.
The quantityQ D,α corresponds, in a collisional model, to the heat current exchanged by the system with the colliding environmental particles. Its mathematical expression, as given by Eq. (11), matches Eq. (41) in Ref. [57]. Similarly, the quantity αQND,α corresponds, in a collisional model, to the work done or produced when switching on and off the interaction of the system with the colliding particles. Its mathematical expression, as given by Eq. (12) matches Eq. (43) in Ref. [57].
Finally, no ambiguity appears when considering the global master equation, since H ND = 0, for which jumps occur between the energy eigenstates.

III. A HEURISTIC APPROACH TO OBTAIN THE ENERGY CURRENTS
In this section, we retrieve the results contained in the previous section by using a semiclassical heuristic approach. In classical stochastic thermodynamics it is quite straightforward to associate energy currents to a stochastic process. Let us consider a system with discrete state space, and let us assume its dynamics to be described by a continuous-time Markov process. The corresponding (classical) master equation readṡ where p j represents the probability for the system to be in state j and W j j are the transition rates from state j to state j . The probability current between any two states reads One can consider the states { j} as lying on a graph where the vertexes are labelled by the state index j. A typical example of such a stochastic process is a particle hopping on a discrete lattice, where p j (t) gives the probability of finding the particle on the site j at time t. The probability current (19) coincides in this case with the particle current between any two sites j and j . Let E j indicate the energy of the state j. A jump j → j occurs thus at the expenses of an energy E j − E j absorbed or injected from/into the surrounding environment. After inspecting Eq. (19), one can then write the heat current flowing from the vertex j to the vertex j along the link j → j . Such a heat current readsQ( j → j ) = (E j − E j )(W j j p j − W j j p j ). Thus the total heat current flowing from/into the state j, while the system interacts with its environment, reads [73] where the sum runs over the set Ω j of nodes connected to the node j.
The quantum analogue of Eq. (18) is Eq. (1) while the quantum analogue of the probability current (19) was introduced in [68]. In the quantum case such a current is composed of two where x j = | j j| is the projector onto the state j. The first component is the current resulting from state-to-state jumps due to the interaction with the environment as embodied by the dissipators (2) and we remind the reader that the definition of ω j j , is given by Eq. (4). Given our choice for the jump operators L λ = | j j|, the following equality holds x j = L † λ L λ . By using this result, the first anticommutator on the RHS of Eq. (23) reads Thus we conclude that Eq. (23) and (8) are equivalent, and we can write this is the first main result of this section, and expresses the fact that part of the energy current operator d t H, as given by Eq. (7), can be expressed in terms of the probability current (21), analogously to what one standardly does in classical stochastic thermodynamics to derive the heat current (20). The corresponding result for J (tun) Q requires a somewhat more elaborate analysis, given that the operator D * α [H ND ], Eq. (9), cannot be directly identified with the current operator J tun Q,α (λ), Eq. (24). As a matter of fact there is no equation, analogue to Eq. (25), relating D * α [H ND ] and J tun Q,α (λ). However we show in the following that such an equality can be found for the expectation values of the two operators. We first notice that a straightforward manipulation gives Furthermore, according to Eq. (5), we can write We thus see that the current (24) can be associated with the coherent part of the dynamics of H ND . Finally comparing Eqs. (24), (26) and (27), we conclude that in the steady state the equality d t H ND where we have used the definition (12) in the last equality. This is the second main result of this section: the tunnelling probability current (22) is associated with an energy current (24), which in turn is equal to the second contribution to the energy current d t H (7) in the steady state. Based on the results of the previous section (specifically on Eq. (9 and the subsequent discussion), we also conclude that λ J (tun) Q (λ) ss corresponds to the steady state work rate that can be obtained within a collisional model framework.

IV. THE CLOCK MODEL
We now apply the results developed in the previous sections to a physical system consisting of a chain of N rotors with N s levels each denoted as |n ∈ {|0 , |1 , |2 ... |N s − 1 }, each level n corresponding to a specific position of the clock's hand. The Hamiltonian of the system, that generalises the spin-1/2 case, can be written as [60]: where the matrices σ and µ are N s dimensioned and have the following form: where ν = exp(iψ), ψ = 2π/N s , the operator σ (σ † ) is a tunnelling term that rotates the spin anti-clockwise (clockwise) while µ is a measurement of the clock's position. The phase φ i j kπ makes the interaction between the rotors i and j chiral.
Of particular interest for our discussion here will be the case τ i = 0, in which there is no local tunnelling term and the Hamiltonian becomes diagonal in the position eigenbasis.
In the specific case of the rotors, the probability currents defined in Eqs. (21) and (22) can be interpreted as rotational currents that express the rotational rates of the rotors.
We will now show that the chiral interaction part of the Hamiltonian (29) proportional to K i, j is not invariant under a specific rotation. In order to fix the ideas let us consider the case of a dimer (N = 2) in a state |n 1 , n 2 , the interaction energy reads U(n 1 , n 2 ) = K 1,2 /2 cos[ψ(n 1 − n 2 ) + φ 12 ]. Let us analyse the energy of the two following states: first, the state |−n 1 , −n 2 , obtained from |n 1 , n 2 after a reflection about the origin of both rotors, and, second, the state obtained with a rotation of the first rotor alone: |n 1 + m, n 2 , where m ∈ Z is an integer, with the prescription that the indexes are taken to be cyclic in the sense that |n k + N s ≡ |n k . A straightforward manipulation, using standard trigonometric equalities, shows that it is possible to find an integer translation m, such that U(−n 1 , −n 2 ) = U(n 1 + m, n 2 ) for any n 1 , n 2 = 0, . . . N s − 1, if and only if φ 12 = ψ/2, with ∈ Z. In this case one finds m = N s − . Thus, the condition φ 12 ψ/2 results in an interaction energy that is not rotationally invariant, in the sense that U(−n 1 , −n 2 ) U(n 1 + k, n 2 ).
The above argument can be generalised to N rotors, by taking a reflection of all the spins about the origin {−n i } and a single spin translation. One then finds the same condition for the phases φ i j . Thus for φ i j ψ/2 the system interaction energy is not rotationally invariant. This broken symmetry, together with the thermal disequilibrium, obtained by putting different rotors in contact with different baths at different temperatures, results in a non-vanishing steady state rotational current J th ss , both for classical [66,67] and quantum [68] systems of rotators. The model (29) will thus behave as an autonomous thermal motor, converting heat currents into mechanical currents (rotations in this specific case).
To simplify the analysis of the dynamics, in the following we only allow jumps between states where only one spin is flipped, for example | j = |n 1 . . . n k . . . n N → | j = |n 1 . . . n k ± 1 . . . n N . We find the system steady state by solving the ME (1) with jumps operators L λ = | j j|, and bosonic bath dissipation rates [37] γ α (|n 1 . . . n k . . . n N → |n 1 . . . n k ± 1 . . . n N ) with ω λ given by Eq. (4). The choice of the dissipation rates (31) entails the transition of the k-th rotor to be driven by the α = k bath alone.

V. DIMER SYSTEM
We begin our discussion by looking at the simplest nontrivial system for which our results can be illustrated, namely the dimer (N = 2, N s = 3, N b = 2), which was extensively studied in [68]. We assume τ i = τ, K 1,2 = K and φ 1,2 = φ.
This system is connected to two baths at temperatures T α = 1/β α (α = 1, 2) . Without loss of generality we assume T 2 > T 1 . We first analyze the thermal and the tunnelling probability currents J th and J tun given by Eqs. (21) and (22), respectively. In Fig. 1, we plot the two currents against the phase φ for two different values of the transverse field τ. In Fig. 2 we plot the currents as a function of τ at fixed φ. We see that the currents are 2π/N s periodic in φ, and the tunnelling current (22) vanishes for τ = 0, as expected.
We now turn our attention to the two energy currentsQ D,α andQ ND,α given in Eq. (11) and Eq. (12), respectively. As discussed in Sec. II, these are the contributions to the total energy current associated with the bath α, Eq. (10), arising from the diagonal and non-diagonal part of the Hamiltonian, respectively. We use the conventionQ x > 0 when a heat current flows from the bath into the system. The heat currents show a similar behaviour as the probability currents depicted in Figs. 1 and 2. In particular in Fig. 3 we plot the heat currents as a function of φ and again we see that these quantities are 2π/N s periodic. By taking the steady state expectation value of both sides of Eq. (7), and keeping in mind the definitions (11)-(12), one finds αQD,α +Q ND,α = 0. Inspection of Fig. 3 shows that when the Hamiltonian (29) is diagonal (τ i = 0), the two heat currentsQ D,1 andQ D,2 sum up exactly to zero, given thatQ ND,α = 0. Conversely, when τ 0, αQD,α = − αQND,α 0 holds, see Fig. 3-(b).
In Fig. 4 we plot, as functions of the transverse field τ, the energy currentsQ D,α andQ ND,α at fixed φ. The steady state results αQD,α +Q ND,α = 0 is also confirmed by this diagram. Incidentally, we remind the reader that the quantity αQND,α can be identified as an input work rate within the collisional model, see discussion in Sec. II.
To see what possible operating modes are available for the dimer system we generate a list of 10,000 parameters choices and calculate the resulting heat currents which are shown in Fig. 5. A special operation regime is achieved in the classical case when τ i = τ = 0, for which αQND,α = 0 and so the current from each bathQ D,α is equal and opposite. Interestingly, we do not find a regime where the dimer works as a refrigerator (Q D,1 > 0) therefore in the next section we move to a trimer system.

VI. TRIMER SYSTEM
We consider now a system of three rotors as depicted in Fig. 6. Again without loss of generality we take T 3 > T 2 > T 1 . The addition of the third bath in the trimer system (N = 3, N s = 3, N b = 3) opens up a lot of new possibilities compared to the dimer case in terms of thermal machine construction, one of the most notable being that of absorption devices [35,43,. These perform some task without the requirement of external work input, thus operating as autonomous devices. For example an absorption refrigerator performs refrigeration without external work. Its performance is measured using the coefficient of performance (COP) which measures the refrigeration power from the cold bath (Q 1 here) with respect to the heat input required from the hot bath (Q 3 ), COP =Q 1 /Q 3 . Absorption refrigerators have been found in a selection of different quantum systems (see for example [5]). The previously conjecture that three-body interactions were necessary for quantum absorption refrigerators has been proven wrong [35]. To help us understand whether a trimer of rotors can be designed so as to work as an absorption refrigerator, we start by considering a specific set up in which two of the rotors are non-interacting, namely in Fig. 6 we take K 2,3 = 0. We also take τ = 0, implyingQ α =Q D,α , the case of non-vanishing transverse field being considered later in this section.
If we inspect the behaviour of a dimer at τ = 0, as exemplified by Figs.1-(a) and 3-(a), we see that if we isolate the dimer 1 − 3 (K 1,2 = K 2,3 = 0), with 0 < φ 1,3 < π/3, with all the other parameters being fixed, the heat flows from 3 to 1, Fig. 3-(a), as expected, while the rotational current of the spin 1 is negative, Fig. 1-(a). Similarly if one considers only the dimer 1 − 2 (K 1,3 = K 2,3 = 0), with −π/3 < φ 1,2 < 0 the rotational current for the rotor 1 is now positive, but the heat current still flows into the cold bath at temperature T 1 . Our first attempt has thus been to choose φ 1,2 > 0 and φ 1,3 < 0, Solid lines simply join the points and are a guide to the eye. so as to have two conflicting effects on the rotational current of rotor 1, that might invert the sign ofQ 1 , thus resulting in a refrigerator.
Based on these arguments, in the following we take φ 1,3 = −φ 1,2 . Fig. 6 shows that absorption refrigeration is indeed possible (Q D,1 > 0) albeit in a rather narrow range of interaction strengths, given the choice of the other parameters.
As one would expect there is a trade-off between the maximum refrigeration power that one can obtain and its COP. This is confirmed by inspection of Fig. 7 where we compare the two quantities.
We then consider the case of non-vanishing field τ 0, and find that adding a coherent term to the Hamiltonian reduces the refrigerator's effectiveness, for any value of τ, see Fig. 8.
We now optimize the functioning of this system as an absorption refrigerator, using a method termed differential evolution [69,70]. The basic essence of this method is to vary some of the parameters of the system with the aim of maximizing a given quantity of interest. In our framework, we focus on the heat currentQ D,1 extracted from the coldest reservoir, and use a differential evolution algorithm to improve its value by changing a (sub-)set of system parameters. A description of the algorithm used to implement this is found in Appendix A. In this cases we fixed the temperature of the baths to T 1 = 1, T 2 = 1.5, T 3 = 2.5 and the rates g i = g = 1 and allowed the other parameters to vary. The resulting optimal parameters are: τ = 0, Notice that the four curves sum up exactly to zero for any τ.
These optimization results yield a much larger refrigeration power thanks to the condition K 2,3 0 compared to what we found previously in Fig. 7 with K 2,3 = 0. These increased power is however at the expense of the COP which is now lower. Interestingly, the optimization outputs the parameter τ = 0 corresponding to the classical regime.

VII. THE TRIMER AS A THERMAL CONTROL DEVICE
The trimer system is interesting from a thermal control aspect as well, with possibilities of rectification and switching [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36]. In the following we consider the linear chain geometry with K 1,3 = 0. We start with the case of a switching device. In the case with τ i = 0, if the middle bath at temperature T 2 is not attached to the corresponding rotor 2, then the latter will be unable to rotate and will stop any current flowing between the baths at temperatures T 1 and T 3 , see Fig. 9. One can therefore use the coupling g 2 between rotor 2 and the middle bath as a switching parameter for the system.
Another control feature of this setup is that of rectification. This phenomenon consists in an asymmetry of the heat flow through the system when the heat baths at temperature T 1 and T 3 are reversed. The rectification effect is measured with the rectification coefficient: is the heat flow when the temperature gradient is left-right (right-left). Notice that no rectification is observed in the dimer case as the system is symmetrical. To increase the degree of asymmetry in the device with three rotors we set the phases φ i, j to be non-uniform.

VIII. CONCLUSION
To conclude, in this work, we have shown how the local ME modelling of an open quantum system is always consistent with the thermodynamics laws. Our proof is based on the structure of the ME itself, and does not resort on any microscopic model, such as the collisional one. By identifying the relevant contributions to system energy evolution in the Heisenberg picture, we have pinpointed the proper heat currents that satisfy the second law of thermodynamics. Our analysis, therefore, leads us to formulate the correct expression of the second law for the case of the local ME. Furthermore, using an intuitive argument, we have recovered the expression of the quantum heat currents by manipulating the expression of the quantum probability currents, a procedure standardly used in classical stochastic thermodynamics.
We have used our general results to study two and three rotor systems. In particular, we have shown that the trimer behaves, in certain parameters regimes, as an absorption refrigerator. We have characterized its performance in terms of ← − Q (dashed). Other parameters: K 1,2 = −24, K 2,3 = 0, K 1,3 = −20, φ 2,3 = 0, g 1 = g 3 = 1, g 2 = 10 −5 , τ i = τ = 0, T 1 = 1, T 2 = 1.5, T 3 = 2.5. refrigeration power and coefficient of performance. In addition we have shown that such a systems can also be used as thermal control devices able to act both as a switch and a rectifier.
replacing them if they provide an improvement. The full algorithm for this is as follows. First, an initial variable vector U i t is generated, a set of N p d-dimensional vectors , where d is the number of parameters being optimised over. For this work we set N p equal to twice the number of parameters to optimise. In addition, the function F to be optimised is chosen. Then, the algorithm proceeds with the following steps: 1. Generate N p mutant vectors M via the formula: where {k, l, m} i are mutually exclusive integers randomly chosen in the interval [1, N p ] 2. From these N p vectors, a new offspring vector O i is gen-erated: where C r is the crossover factor.
3. The offspring vector is compared to the current parameter vector and is replaced if yileding a better value for the cost function: 4. The process is repeated until convergence.