Non-additive dissipation in open quantum networks out of equilibrium

We theoretically study a simple non-equilibrium quantum network whose dynamics can be expressed and exactly solved in terms of a time-local master equation. Specifically, we consider a pair of coupled fermionic modes, each one locally exchanging energy and particles with an independent, macroscopic thermal reservoir. We show that the generator of the asymptotic master equation is not additive, i.e. it cannot be expressed as a sum of contributions describing the action of each reservoir alone. Instead, we identify an additional interference term that generates coherences in the energy eigenbasis, associated with the current of conserved particles flowing in the steady state. Notably, non-additivity arises even for wide-band reservoirs coupled arbitrarily weakly to the system. Our results shed light on the non-trivial interplay between multiple thermal noise sources in modular open quantum systems.


I. INTRODUCTION
An improved understanding of the dynamics of open quantum networks is desirable for research elds including quantum thermodynamics, 1-4 quantum biology, 5 mesoscopic electronics 6 and the theory of non-equilibrium phase transitions. [7][8][9][10][11][12] The Lindblad master equation 13,14 is a popular and powerful tool for modelling such systems, which approximates the long-time dynamics under the assumption of weak coupling to memoryless environments, i.e. the Born-Markov approximation (BMA). Additional, uncontrolled approximations are typically required in order to obtain a completely positive evolution, such as the secular, 15 singular-couplinglimit 16 or wide-band-limit 17,18 approximations. However, in composite open systems, distinct approximations may lead to di erent and sometimes drastically incorrect predictions.
In particular, standard master equations derived under the BMA describe dissipation in terms of either local processes on a small number of sites or global transitions between energy eigenstates of the entire network. 19 Unfortunately, the local approach appears to violate thermodynamic laws 20,21 for certain kinds of time-independent system-bath interactions, 22,23 while the global approach fails to capture the coherences necessary to properly describe the non-equilibrium steady state (NESS). 24 This poses a particular problem for the ongoing study of quantum thermal machines, where the multifaceted role played by coherence-acting variously as a performance-enhancing resource, 25-32 a useful output [33][34][35][36] or an unavoidable hindrance 37,38 -remains incompletely understood.
Here, we aim to elucidate these issues by studying an exactly solvable model of a non-equilibrium quantum network. More precisely, we derive a time-local master equation describing a pair of coupled, localised fermionic modes. Each mode exchanges particles and energy with a macroscopic reservoir represented by a semi-in nite, uniform tightbinding chain, such that the total Hamiltonian is quadratic in fermionic ladder operators. This represents a prototypical quantum thermal machine, which could be realised by electrons owing through a serial double quantum dot 39 or cold fermionic atoms con ned by an optical lattice. 40,41 a) markTmitchison@gmail.com Note that other authors have already derived and extensively studied exact time-local master equations describing general quadratic Fermi systems. [42][43][44][45] In contrast to these previous approaches, the simplicity and symmetry of our speci c set-up enables compact analytical solutions to be obtained using elementary methods and without needing the Born-Markov, secular, singular-coupling-limit or wide-bandlimit approximations. Instead, we use an alternative, perturbative approximation scheme, valid when the systemenvironment coupling is much smaller than the bandwidth of the reservoir vacuum noise, and explicitly con rm its accuracy using the exact solution. This simpli cation allows us to clearly identify and distinguish the di erent physical processes governing the dynamics.
Our analysis challenges a central tenet of standard opensystems theory: namely, that two initially uncorrelated environments should give rise to two independent, additive contributions to the master equation in the weak-coupling limit (see, for example, Ref. 46). Instead, we distinguish three contributions to the generator of asymptotic time evolution. Two of these describe the individual thermalising e ect of each reservoir, while the third is an interference term arising from the combined action of both reservoirs whenever they are initially out of equilibrium with each other. This interference gives rise to coherence in the energy eigenbasis of the network, a feature which re ects the conserved fermion current owing in the NESS. These ndings connect several previous theoretical studies on composite open quantum systems. In particular, a signi cant body of research has sought to assess the validity of approximate, additive master equations by comparing them to each other or to exact numerical results. [19][20][21]24,[47][48][49][50][51][52][53] Viewed broadly, these investigations indicate that the local and global Lindblad equations are only accurate in limited, complementary parameter regimes, and only for certain observables, even though the requisite conditions for the BMA might hold for each bath individually. On the other hand, a few recent papers have shown that the additivity assumption fails in the presence of multiple strong or structured noise sources. [54][55][56] We extend these results by demonstrating that interference between di erent thermal baths out of equilibrium gives rise to non-additive dynamics, even when the open system couples arbitrarily weakly to spectrally unstructured reservoirs. The interference contribution is not in Lind-N blad form and thus underlies the occurrence of asymptotic non-Markovianity found by Ribeiro et al., 45 where the longtime dynamics cannot be described by a Markovian master equation even when the open system has lost all memory of its initial state. [57][58][59] Moreover, we show that our nonadditive master equation interpolates between the local and global Lindblad equations and recovers them in di erent limits. Since additivity is a necessary consequence of the Born-Markov approximation, 56 this implies the failure of the BMA for describing steady-state transport away from these limits. Nevertheless, we nd that certain observables-such as the steady-state currents owing into the baths-are accurately predicted by an additive master equation across a relatively broad range of parameters. Our work thus helps to clarify the validity of existing Lindblad models, while providing a reference point for future research aimed at moving systematically beyond the standard approximations.
The plan of the paper is as follows. Section II introduces some preliminary concepts and overviews our main results. The microscopic model that forms the core of this work is de ned and exactly solved in Section III. We derive the exact master equation describing the system and develop a weak-coupling approximation scheme in Section IV. Nonadditivity of the asymptotic dynamics is explored in Section V. We discuss our results and conclude in Section VI.

II. PRELIMINARIES
To motivate the problem at hand, consider an open quantum system S comprising a network of coupled sites. This network is coupled weakly to multiple thermal reservoirs which may exchange particles and energy with S, as illustrated in Fig. 1. Each reservoir B α is characterised by a temperature T α = 1/β α and chemical potential µ α (we work in units where k B = 1 and = 1). LetĤ S andN S respectively be the Hamiltonian and particle number operators on S, where [Ĥ S ,N S ] = 0 unless µ α = 0. We assume for simplicity that H S has a non-degenerate spectrum. The quantum state of the network at time t is denotedρ S (t). We suppose that initially each of the baths is not correlated with the others nor with the open system, and that S asymptotically approaches a unique stationary stateρ ∞ S = lim t→∞ρS (t). We now summarise the behaviour that we expect in general.
One fundamental property of a thermal reservoir is that a small system weakly coupled to it should eventually equilibrate to the same temperature and chemical potential. In addition, the composition of several independent thermal reservoirs in equilibrium itself constitutes a thermal reservoir. Hence, for equal reservoir temperatures and chemical potentials the system should equilibrate, i.e.
On the other hand, if the reservoirs are initially out of equilibrium with each other, the imbalance in temperature or chemical potential sets up a current of energy or particles owing into S from B α . These energy and particle currents are respectively denoted J E α and J P α , while J Q α = J E α − µ α J P α is the corresponding heat current. In the stationary state, the second law of thermodynamics requires that the total rate of entropy production in the reservoirs is non-negative: Assuming that di erent reservoirs couple to di erent regions of the network, basic conservation laws imply that the transfer of energy or particles between the reservoirs can only occur via commensurate energy or particle currents owing within the system. Therefore, there exist one or more current observables on S, here denoted schematically byĴ S , having non-zero expectation value in the NESS: (We denote expectation values at time t by • t , with • ∞ the limiting value as t → ∞.) The existence of such internal currents in a boundary-driven system implies that the nonequilibrium steady state must exhibit coherence in the eigenbasis ofĤ S . Although we defer a detailed demonstration and discussion of this claim to Appendix A, its plausibility can be appreciated by considering the example of a one-dimensional (1D) network with open (i.e. non-periodic) boundary conditions. For this geometry, the eigenstates ofĤ S do not support internal currents at all. 24 This follows because, in the absence of external sources or sinks, any such current would lead to an accumulation of particles or energy in one part of the system, which is incompatible with the fact that energy eigenstates are stationary states of the closed-system dynamics.
A widely used dynamical model of the situation depicted in Fig. 1 is the quantum master equation valid for times t much greater than the environment memory time. In order for Eq. (4) to generate a completely positive and trace-preserving (CPTP) evolution for any stateρ S (t), the dissipator L must be in Lindblad form 13,14 where D[L]ρ S =Lρ SL † − 1 2 {L †L ,ρ S } andL j is a jump operator describing an incoherent transition occurring at a rate γ j .
Since the reservoirs are assumed to initially be statistically independent, the standard construction of the dissipator is a sum of generators L α representing each bath B α , i.e.
This constitutes our de nition of additivity, as studied previously in Refs. 54-56. This is distinct from the concept of additive decoherence rates explored, for example, in Refs. 60 and 61. The additivity assumption permits one to unambiguously identify the particle current J P α (t) and energy current J E α (t) entering the system from B α as Here, L † α is the adjoint generator describing the Heisenbergpicture evolution of observables, de ned by Tr[BL † αÂ ] = Tr[ÂL αB ] for arbitrary operatorsÂ andB.
Regarding the speci c form of the generators L α , various inequivalent approaches are commonly employed. These can be broadly classi ed into two groups according to the xed point of each generator, i.e. the stater α satisfying L αrα = 0. The rst is the "global" approach, where each generator drives the entire network towards the corresponding equilibrium state, i.e.
Note that here we assume the xed pointr α to be unique. Lindblad generators satisfying Eq. (7) can be derived from a microscopic model under the Born-Markov and secular approximations. 15 The latter approximation is assumed to be justi ed in the quantum-optical regime, where the separation between energy levels ofĤ S is much larger than their environment-induced broadening. Eq. (7) directly implies the correct equilibration behaviour (1), and also ensures-via the Spohn inequality 62 -that the second law (2) holds for the heat currents from the baths. However, generators derived under the aforementioned approximations also have the property that they do not couple populations and coherences in the eigenbasis ofĤ S . 15 Adding such generators together according to Eq. (5) therefore generates independent equations of motion for the populations and coherences. Since the evolution is trace-preserving, if the stationary state is unique then it must be diagonal in the eigenbasis ofĤ S , which is generally inconsistent with the condition (3). 24 Alternatively, one can use a "local" approach, where each generator L α acts non-trivially only on one part of the network s α ⊂ S, driving it towards thermal equilibrium while leaving its complements α una ected. That is, whereĤ s α andN s α are respectively the Hamiltonian and particle number operator of s α , whileÔs α is an arbitrary density operator with support on the complements α . Clearly,r α is not unique in this case.
Generators satisfying Eq. (8) may either be derived microscopically using various approximations, 19,24,47,49 or directly postulated on phenomenological grounds. [63][64][65][66] In the local approach,ρ ∞ S is not necessarily diagonal in the energy eigenbasis, and therefore provides a consistent model for the internal current dynamics as required by Eq. (3). However, neither Eq. (1) nor Eq. (2) hold, in general, leading to potential violations of thermodynamic laws. 20,21 In what follows, we consider the simplest case of a twosite network with two independent reservoirs labelled by the index α = L, R. We show that, in the limit of weak systemreservoir coupling, the asymptotic dynamics is governed by a master equation whose dissipator takes the form In the quantum-optical limit, L L and L R determine the currents from the baths according to Eq. (6) and induce thermalisation according to Eq. (7). Outside of the quantum-optical limit, these generators drive S towards the reduction of a global thermal state (i.e. including the baths and the systembath interaction) that accounts for system-reservoir correlations. 67 The interference term L int , which appears whenever B L and B R are not in equilibrium with each other, generates coherence in eigenbasis ofĤ S as required by condition (3). In this way, all three properties (1)-(3) can be satis ed, but only by abandoning the additivity assumption (5). Physically, non-additivity stems from correlations between the two reservoirs. 54 Out of equilibrium, such correlations grow steadily in time due to the particle current owing between the baths. Hence, Eq. (5), which is justi ed by the initial statistical independence of the reservoirs, fails to hold as t → ∞, even if the reservoirs are spectrally unstructured and coupled arbitrarily weakly to the system.

III. EXACTLY SOLVABLE MODEL
A. Description of the model Let us now detail our speci c set-up, which belongs to the well-known family of resonant-level transport models. 68,69 The total Hamiltonian takes the formĤ =Ĥ S +Ĥ B +Ĥ SB , describing a central open system S sandwiched between two fermionic particle reservoirs B L and B R , as illustrated in Fig. 2. This models a thermoelectric tunnel junction 4 or entangler, 33,34 which channels a current of fermions between two conducting leads due to a temperature or chemicalpotential gradient.
The system and bath interact via tunnelling of fermions between the terminal site of each reservoir and the adjacent site of S. We parametrise the tunnelling energy as √ ΓΩ, where Γ sets the overall frequency scale of the dissipative dynamics. Explicitly, the interaction Hamiltonian reads aŝ with tunnel couplings τ q = √ ΓΩ/(2M + 2) sin(q).
It is convenient to analyse the problem in a basis that di-agonalisesĤ S . To do this, we collect the ladder operators into column vectorsĉ = (ĉ L ,ĉ R ) T andâ q = (â q,L ,â q,R ) T . We then de ne a new canonical set of ladder operatorŝ where ∆ = g 2 + δ 2 . The Hamiltonian hence splits into two independent piecesĤ =Ĥ 1 +Ĥ 2 , witĥ where ∆ are the single-particle energy eigenvalues ofĤ S . For concreteness, we assume that E j > 0, or equivalently that h α > 0 and g < 2 √ h L h R . The e ect of the reservoirs on the system is determined by the spectral density This is a smooth function of ω in the limit M → ∞, where the spacing between adjacent wave vectors ∆q = π/(M +1) tends to zero and q becomes a continuous variable taking values in the rst Brillouin zone q ∈ [0, π]. Using the prescription q ∆q → dq, we obtain where Θ(x) is the Heaviside unit step function. We label this spectral density with a subscript N after Newns, who (to our knowledge) introduced it. 70 According to Eq. (19), each environment is characterised by a spectral bandwidth Ω, leading to a vacuum correlation time of order Ω −1 . This is rather intuitive, since Ω sets the rate at which an excitation created at the boundary of B α propagates irreversibly along the chain and away from the central system's domain of in uence. Note that choosing a 1D geometry for each environment is not as restrictive an assumption as it may appear. This is because an environment with arbitrary geometry can be mapped onto a 1D tight-binding model coupled to the system at a single boundary site, so long as the reservoir and interaction Hamiltonians take the generic forms (12) and (15). [71][72][73] In general, the resulting 1D chain is described by inhomogeneous inter-site couplings and local site energies. This leads to scattering of excitations back towards the system, potentially giving rise to recurrences or other non-Markovian effects. In contrast, such backscattering is absent in the uniform chain considered here, which is characterised completely by just two frequencies Γ and Ω.
Directly setting Ω −1 = 0 corresponds to the wideband-limit approximation, 18 which leads to a frequencyindependent spectral density. In the following, we compute the solutions for nite Ω, which enables us to retain energydependent damping rates even for weak coupling, Γ Ω, since we need not assume that E j Ω.

B. Formal solution
In this section, we provide the exact solution for the opensystem density matrixρ S (t) = Tr B [ρ(t)], whereρ(t) is the global quantum state at time t. The same solution formally applies to the general scenario depicted in Fig. 2 (a), i.e. any pair of environments described by Hamiltonians of the form (12) and (15), corresponding to a spectral density (18). We rst describe the formal solution for this general case, before specialising to the Newns spectral density (19) of a uniform chain in later sections. Details of the calculation are presented in Appendix B.
We consider factorised initial conditions of the form Here,N α = qâ † q,αâq,α is the number of particles in B α . The dynamics preserves the total number of fermions due to the relation Therefore,ρ S (t) is characterised by ve independent real numbers, in general. Four of these are encapsulated by the Hermitian correlation matrix which satis es 0 ≤ C(t) ≤ 1. (Here, and throughout this document, the elements of a matrix A are denoted by A jk .) The fth degree of freedom is the double-occupancy probability We compute Eqs. (21) and (22) by solving the equations of motion ford j (t) andb q, j (t) in the Laplace domain. Here, O(t) = e iĤtÔ e −iĤt denotes the Heisenberg-picture time dependence of an operatorÔ. The solution for t > 0 can be completely expressed in terms of the propagator where |0 is the vacuum state, i.e.N |0 = 0. According to Eq. (17), the particle number in each j sector is separately conserved: The function ϕ j (ω) is the probability distribution of excitation energies associated with the stated † j |0 , i.e.
where |n is a single-particle eigenstate ofĤ with energy ε n , i.e.Ĥ |n = ε n |n (note that ε n may be negative 74 ). Therefore, (For general spectral densities, the limit Γ → 0 here means that all tunnel couplings τ q → 0.) The solution can be compactly represented in terms of the matrix With this notation, the correlation matrix reads as The double-occupancy probability is We see that, so long as lim t→∞ G j (t) = 0, the system relaxes to a unique stationary stateρ ∞ S that is independent of the initial conditions. Moreover, thisρ ∞ S is Gaussian 75 because D(∞) = det C(∞). The stationary state is therefore determined completely by its correlation matrix Here, we de ned the Hilbert transform of ϕ j (ω), where P denotes the Cauchy principal value. The foregoing equations, which hold for any spectral density (18), constitute the complete formal solution given knowledge of the propagator G(t). However, in general, computing the propagator is a challenging problem.

C. Propagator for uniform-chain environments
In order to nd an explicit expression for the propagator, we now specialise to environments modelled by uniform chains with spectral density (19). Relaxation to a unique steady state is guaranteed in this case if we make the additional, technical assumption that N Physically, this inequality requires that the energy levels of the open system lie well within the reservoirs' energy bands, so that any initial excitation can eventually be absorbed. Under this condition, we show in Appendix B that where we de ned shifted energies and decay rates The propagator is then evaluated using Eq. (24). We also prove in Appendix B that

D. Steady-state observables for uniform-chain environments
We now present solutions for some interesting observables in the asymptotic stationary state, using the explicit formulae quoted in the previous section. If the baths are initially in equilibrium, i.e. β α = β, µ α = µ and f α (ω) = f (ω), then which di ers from the strict equilibrium value f (E j ) due to the nite width of the energy distribution ϕ j (ω). Using Eqs. (25) and (38), we show in Appendix C that, for equilibrium baths,ρ ∞ S is the reduced state of the global Gibbs ensembleρ We thus infer that the energy-level broadening represented by Eq. (38) arises from the system-bath correlations re ected in Eq. (39), both being related to a nite system-reservoir coupling Γ. With the help of Eq. (27), we recover thermalisation in the strict weak-coupling sense (1) in the limit Γ → 0. Out of equilibrium, particle and energy currents ow from each bath B α into S. These currents are de ned respectively Fermion conservation demands that a corresponding particle currentĴ P S ows between the two sites of the system, which is de ned to satisfy the continuity equations, e.g. ∂ tnL (t) =Ĵ P L (t) −Ĵ P S (t). Explicitly, we havê  Note that the mean intra-system current Ĵ It follows from the de nitions that, in the stationary state, the particle currents are homogeneous throughout the system, i.e. J P The asymptotic currents are calculated in Appendix D. For the Newns spectral density (19), we obtain the standard Landauer formulae 6 with the transmission function The transmission probability is proportional to the overlap between the energy distributions of the two eigenmodes of H S [see Eq. (25)]. Note that the Landauer formulae imply the second law (2) for any transmission function T (ω) ≥ 0. 4 We plot the transmission function in Fig. 3. For Γ ∆, T (ω) is a bimodal distribution that is well approximated by As Γ is increased, the two peaks at ω = E 1,2 broaden and ultimately merge into a single maximum for Γ ∆.

A. Exact master equation
In order to analyse the non-additive properties of the dynamics, we rst write Eqs. (29) and (31) in di erential form, M P 7 corresponding to an exact time-local master equation for the system density operator (see Refs. 43 and 44 for alternative derivations). All equations presented in this section hold for an arbitrary spectral density (18). Assuming that G(t) is non-singular, we de ne Hermitian matrices H(t) and Γ(t) such that iH(t) + 1 2 Γ(t) = −G −1 ∂ t G. We also de ne rate matrices Λ ± (t) by where time arguments are suppressed. The matrix elements of Λ ± correspond to the gain and loss rate coe cients appearing in the master equation, as will be seen shortly. Indeed, upon di erentiating Eqs. (29) and (31), some tedious algebra reveals that Direct comparison con rms that these equations of motion are equivalent to the master equation whereĤ S (t) =d † H Td and we de ned the dissipator After diagonalising the rate matrices Λ ± (t) by a unitary rotation , we cast the dissipator into canonical form 57 Here, the time-dependent jump operators are de ned aŝ These describe loss and gain of excitations from and into modes determined by the eigenbases of the matrices Λ ± (t). Note that i.e. the loss and gain modes di er, in general.

B. Exponential-propagator approximation
Throughout the rest of the paper, we specialise to the Newns spectral density (19) and assume that Γ Ω. Although the latter condition is usually deemed necessary for the Born-Markov approximation to hold, we shall see that it is by no means su cient.  In order to simplify the subsequent discussion, we introduce an approximation scheme where terms of order O(Γ/Ω) are neglected. This amounts to the replacement For the sake of clarity and concision, we henceforth refer to this as the exponential-propagator approximation (EPA).
The EPA is justi ed in Appendix E, where we give an explicit expression for the error in G j (t) thus incurred. We also derive a rigorous upper bound on the magnitude of this error that is proportional to Γ/Ω and decays to zero as t → ∞.
To the same order of approximation, we write E j ≈ E j and Γ j = 2πJ N (E j ). Unfortunately, we have not been able to derive a bound on the error induced by calculating general expectation values such as the correlation matrix (29) within the EPA. Nevertheless, a direct numerical comparison shows that, for suciently small Γ/Ω, the EPA gives an excellent approximation to both the transient and steady-state dynamics, even if E j is comparable to Ω. We illustrate the agreement between the approximation and the exact solution in Fig. 4 for a few example parameters.
We emphasise that the EPA only requires that the coupling Γ is weak in comparison to the environment's energy scale Ω, so that the approximation (53) becomes exact in the wideband limit Ω → ∞. On the other hand, the relation between Γ and the system energy scales E j is not restricted.
The master equation (50) takes a simple form under the EPA. In particular, we have that The loss and gain modes de ned by Eq. (52) can be identi ed by inspection of Eq. (54) in two special cases. First, assume that the two baths are in thermal equilibrium with each other, so that f L (ω) = f R (ω) = f (ω). In that case, F(ω) = f (ω)1 is proportional to the identity, and Λ ± (t) are both diagonal. Hence,L − j = (L + j ) † =d j , i.e. excitations are pumped into eigenmodes ofĤ S , driving the system towards N a state that is diagonal in the energy eigenbasis of S. If the system-environment coupling is su ciently weak, this ensures the proper thermalisation behaviour (see Section V A).
Second, consider the white-noise limit, with Ω → ∞ and F(ω) = F a constant matrix. We then have that Φ(t) = ΓFδ(t) and Γ = Γ1, and since G(0) = 1 it follows that Λ + (t) = ΓF and Λ − (t) = Γ(1 − F). Therefore, the rate matrices Λ ± can be diagonalised by the rotation R T , leading to Lindblad opera-torsL − j = (L + j ) † withL − 1,2 =ĉ L,R , i.e. energy is pumped into localised modes. Of course, this only holds exactly in the unphysical scenario of in nite energy density in the baths, corresponding to β j → 0 or β j , |µ j | → ∞. Nevertheless, the local master equation may be an excellent approximation for sufciently large chemical-potential bias or temperature. According to Eq. (54), the rate matrices are determined by the product of the propagator G(t)-which is diagonal in the energy eigenbasis-and the noise kernel Φ(t)-which is diagonal in the site basis. Hence, in general, the gain and loss modes correspond neither to the energy eigenbasisd j nor the site basisĉ α , but instead lie somewhere in between. In particular, whenever the baths are not in equilibrium with each other, the rate matrices Λ ± are non-diagonal and the dissipation generates some coherence in the eigenbasis ofĤ S . This can be seen, for example, in Fig. 4, where the coherence can be comparable in magnitude to the populations and has a large imaginary part, re ecting the current owing through the system.

V. ASYMPTOTIC NON-ADDITIVITY
A. Non-additivity in the energy eigenbasis Now we demonstrate that the asymptotic dynamics as t → ∞ is not additive. We write the rate matrices in the limit simply as Λ ± = lim t→∞ Λ ± (t). Explicitly, these have the components and Λ − jk = Γ j δ jk − Λ + jk . It follows that the asymptotic dissipator L = lim t→∞ L(t) admits the decomposition Here, the generators L α represent the thermalising e ect of each individual bath α = L, R on the populations, while the interference term L int describes the generation of coherence due to the combined e ect of the non-equilibrium baths, as described below.
First, let us examine the generators L α , which take the form The decay and gain rates are de ned by where Γ L, j = Γ j R 2 1 j and Γ R, j = Γ j R 2 2 j . One readily veri es that the unique xed pointr α satisfying L αrα = 0 is Gaussian, with the correlation matrix In the quantum-optical limit where Γ ∆, this describes an approximately thermal distribution, up to corrections due to level broadening as discussed in Section III D. In fact, using the arguments given in Appendix C, one can show thatr α is the reduction of a global Gibbs state in equilibrium with the corresponding bath, i.e.
Another interesting property of the generators L α in the quantum-optical limit is that they accurately reproduce the steady-state currents according to Eq. (6). Indeed, a direct computation yields with α L † αNS ∞ = 0 = α L † αĤS ∞ . These expressions can be shown to be equivalent to the exact Landauer formulae (43) and (44) using Eq. (46) and writing E j ϕ j (ω) ≈ ωϕ j (ω), which is a valid approximation in the quantum-optical regime.
Since the generators L α are in Lindblad form, they obey the Spohn inequality 62 However, L int does not, by itself, generate a positive evolution, and thus does not necessarily satisfy such an inequality. 76 We nonetheless show in Appendix F that, in the weakcoupling limit, from the Spohn inequality one may recover the second law of thermodynamics in the form with the heat current de ned by J Q α = L † α (Ĥ S − µ αNS ) ∞ . Here, it is necessary to rst divide by Γ before taking the limit in order to avoid recovering the trivial equality M P 9 lim Γ→0 α β α J Q α = 0 implied by lim Γ→0 J Q α = 0. Corrections to the LHS of inequality (64) for nite Γ are quoted explicitly in Appendix F. Now we turn to the interference contribution, de ned by Here, Λ + 12 = Λ + Therefore, L int is associated with the di erence in distribution functions f α (ω). In particular, L int is non-negligible unless f L (ω) ≈ f R (ω) over the entire frequency range in which ϕ j (ω) di ers appreciably from zero. The interference term satis es the properties Hence, L † int generates coherence while leaving the populations una ected. This is to be expected, since energyeigenbasis coherence is an intrinsic property of the NESS of a quantum network, as discussed in Section II. It is noteworthy that, in the Schrödinger picture, L int couples populations and coherences. Such a contribution is therefore precluded in the standard global Lindblad approach due to the secular approximation, which enforces decoupling of populations and coherences. 15 Nevertheless, L int is generally a signi cant contribution even in the quantum-optical limit where the secular approximation is widely believed to be valid. We note also that L int can never be written in Lindblad form. Therefore, departures from Markovian evolution can be used to detect non-additivity, as discussed in Section V C.

B. Non-additivity in the local basis
It is also possible to investigate the violation of additivity in a local, rather than global, picture of dissipation. We transform to the site basisĉ α and decompose the asymptotic dissipator as The local dissipators are Lindblad generators that act only on a single site, given bȳ whereγ ± α = 2 j=1 γ ± α, j . The remaining contribution to L describes delocalised, incoherent processes acting on both sites together: The cross-termL LR is associated with the di erence between the frequency distributions ϕ j (ω) of the open system's energy levels. That is,L LR re ects the extent to which the occupation numbers f α (ω) of reservoir states di er between the distinct frequency ranges sampled by the two distributions ϕ j (ω). In particular,L LR is negligible only in the white-noise limit, where ∆ Ω and the distribution functions f L,R (ω) are essentially constant over the range in which ϕ j (ω) are non-zero.
Note that here the local generatorsL α can be meaningfully associated with bath B α only in the sense that each one depends only on the variables of B α and acts non-trivially only on site α of the system. Even in the weak-coupling limit,L L,R do not obey the properties (6) and (8) expected of additive, thermal dissipators unless δ ≈ 0, as discussed below.
Let us rst examine the validity of Eq. (8). The xed point r α satisfyingL αrα = 0 is of the form where, by analogy with Eq. (8), we de ned an e ective local Hamiltonian acting on site α, whileÔ α is an arbitrary density operator with support on the other site α α. Using Eq. (27), we nd This only has the detailed-balance form required for thermalisation if β α ∆ 1 or if |β α (h − µ α )| 1, which corresponds to the white-noise limit. In that case, we have that Second, let us discuss the de nition of the currents via Eq. (6). For the particle current, we have, for example, where we identi ed L †n L ∞ = J P L using the exact master equation (50). Hence, Eq. (79) di ers from the true particle current by an amount Clearly, this correction is negligible ifL LR ≈ 0 or if ∆/Ω → 0 so that Γ 1 ≈ Γ 2 . In addition, Eq. (80) vanishes as δ → 0, since one can easily show that Re ĉ † Lĉ R ∞ ∝ δ/∆. On the other hand, we nd for the energy current This expression is only correct whenL LR ≈ 0 and δ h, g, such that L † LĤ S ∞ ≈ hJ P L . This approximately agrees with the exact Landauer formulae (43) and (44) in the white-noise limit, where Ω → ∞, T (h + ω) = T (h − ω) and f α (ω) ≈ const.

C. Non-Markovianity witnesses non-additivity
Since neither L int norL LR are in Lindblad form, it may not be possible to cast the overall dissipator L into Lindblad form for certain values of the parameters. Therefore, its canonical representation (52) may exhibit one or more negative rates λ ± j (t) even as t → ∞. Negativity of the rates λ s j (t) signals the onset of non-Markovian evolution, according to the measure of non-Markovianity based on completely-positive divisible maps proposed by Rivas et al. 77 In this framework, the degree of instantaneous non-Markovianity is quanti ed by 57 Asymptotic non-Markovianity (ANM) corresponds to lim t→∞ ν(t) = ν ∞ > 0. In a previous work, Ribeiro et al. 45 demonstrated that ANM arises in a class of spin and fermionic transport models, which includes the set-up considered here within the wide-band-limit approximation Ω −1 = 0. ANM has also been recently identi ed in several other open-system scenarios. [57][58][59] We plot ν ∞ in Fig. 5, as calculated using the EPA rate matrices de ned by Eq. (55). At low temperatures, we nd ν ∞ > 0 whenever the chemical potential of one reservoir lies in between the two energy eigenlevels of the system while that of the other reservoir lies outside. This e ect is signi cantly reduced once the thermal energy grows comparable to the splitting ∆. It is remarkable that a highly non-Markovian evolution can be induced merely by tuning a macroscopic potential di erence, without any engineering or microscopic control of the environment Hamiltonian. Indeed, the ANM e ect survives even in the wide-band limit Ω → ∞ where each environment is completely unstructured. In Fig. 5(c) we also show that some non-Markovianity may be generated by purely thermal driving, i.e. where the reservoirs have identical chemical potentials and di erent temperatures.
Since the generators L L,R (L L,R ) are in Lindblad form, the appearance of asymptotic non-Markovianity is associated with the cross-term L int (L LR ). Hence, ANM witnesses the violation of additivity in either the global or local pictures of dissipation. Indeed, in Fig. 5 we observe that non-Markovianity is associated with parameters for which the master equation is not additive, in neither the local nor energy eigenbases. As discussed in previous sections, this occurs when the baths are far from equilibrium, yet the chemical-potential bias and the temperature are not so large that the white-noise limit is recovered. Since ANM is, in principle, experimentally accessible via quantum process tomography, it represents a measurable signature of non-additivity. M P 11

VI. CONCLUSION
In this work, we have presented a study of an exactly solvable microscopic system-environment model, featuring a two-site fermionic network driven out of equilibrium by two independent thermal baths. We have derived the exact master equation describing the open-system density matrix and developed a simplifying approximation scheme valid when the system-reservoir coupling Γ is much smaller than the reservoir bandwidth Ω.
Our results demonstrate that the asymptotic master equation is not additive, i.e. the dissipator cannot be written as a sum of the form (5) where each term pertains to a single bath, even if the system-reservoir coupling is vanishingly small in comparison to every other relevant energy scale. Examining the master equation in the energy eigenbasis of the open system, we identify a generator L α associated to each bath, which drives the open network towards thermal equilibrium according to Eq. (7) [or, more generally, Eq. (60)]. Furthermore, L α determines the currents via Eq. (6) in the quantumoptical limit. This generator can therefore be interpreted as an individual contribution describing the e ect of bath B α acting in isolation.
Nevertheless, additivity is violated away from equilibrium due to an additional interference term L int that generates coherence, as required for a current-carrying NESS in an extended system. In extreme cases, this non-additivity can lead to an asymptotically non-Markovian master equation. Despite this, the interference term L int leaves the energyeigenbasis populations invariant and therefore the global, additive Lindblad model obtained by setting L int = 0 still leads to accurate approximations for the boundary currents in the quantum-optical regime. Such an additive model coincides in the limit Γ → 0 with the standard master equation derived under the BMA and secular approximation.
We have also examined the violation of additivity in the local basis. However, here it is not always possible to identify thermal generators satisfying the properties (6) and (8), even in the weak-coupling limit. Thus, the identi cation ofL α as an "individual" bath contribution is not fully justi ed, since it does not accurately represent the e ect of a single bath acting in isolation. Nevertheless, an additive, local Lindblad model is recovered when the modes are near-resonant, i.e. δ h, g, and when the noise correlation time is much smaller than the inverse frequency splitting ∆ −1 . This conclusion is entirely consistent with recent studies showing that local additive Lindblad models may perform well in non-equilibrium scenarios when the network nodes are nearly degenerate, 52,53 but fail for large detunings. 20 It has been shown 56 that additive open-system dynamics is obtained whenever the BMA (or, more generally, the secondorder time-convolutionless projection operator method 15 ) holds. Hence, our results indicate a breakdown of these assumptions at asymptotically large times when multiple thermal reservoirs act in competition. Note that, while the BMA is commonly cast in terms of an assumption that the systemreservoir density matrix factorises, it is more accurately de-scribed as a projection P onto product states of the form whereρ B is the initial state of the environment, in combination with a low-order perturbation expansion in the systemreservoir interaction. 19 The breakdown of additivity for arbitrarily weak systembath coupling Γ indicates that Eq. (83) does not serve as a good reference point for a perturbative expansion in powers of Γ. This can be understood as a consequence of the correlations that build up between the two parts of the environment due to the current owing between them. As time increases, these correlations cause the joint stateρ(t) to move arbitrarily far away from the subspace spanned by states of the form (83) withρ B =ρ LρR . A very interesting open question is whether a projection of the type (83) can still facilitate a good perturbative description of the long-time dynamics, but with a correlated reservoir stateρ B .
Although the evidence presented here pertains only to the speci c system considered, we expect the qualitative conclusions regarding non-additivity to hold for other extended open systems, where currents and energy-eigenbasis coherence are inextricably linked. Indeed, previous results demonstrating asymptotic non-Markovianity in larger 1D fermionic and spin networks driven out of equilibrium 45 already support this conclusion, because ANM witnesses non-additivity, as we have shown. Furthermore, our model can be readily generalised to a bosonic setting where the ladder operators obey commutation, rather than anti-commutation, relations. Since the Heisenberg equations are essentially identical in this case, one expects the same conclusions regarding additivity to hold, although we leave a detailed analysis of this problem to future work (see also Refs. 78 and 79). We note also a very recent study of a di erent bosonic model of coupled mechanical oscillators that displays coherences in the NESS, 80 which are naturally explained by the interference between non-equilibrium baths.
Looking ahead, it will be interesting to investigate the consequences of the results presented here for the thermodynamics of small heat machines running between multiple thermal reservoirs. For example, it remains to be seen how non-additive noise may a ect the thermodynamic power and e ciency of such machines, or their ability to generate quantum resources such as coherence and entanglement. One may also ask whether the interference between the reservoirs is manifested in the uctuations of the energy and particle currents. Finally, the framework presented herein appears to be an ideal setting to explore strong-coupling effects in thermodynamics, 81-86 since system-environment correlations in both equilibrium and far-from-equilibrium states may be taken into account. N Andreas Lemmer, Fabio Mascherpa, Martí Perarnau-Llobet, Ralph Silva, Andrea Smirne and Philipp Strasberg. This work was funded by the ERC Synergy grant BioQ and the EU project QUCHIP.

APPENDICES Appendix A: Connection between conserved currents and coherence
In this appendix, we elucidate the relationship between conserved currents and coherence in the energy eigenbasis. Consider an open quantum network S connected to several baths B α , such that the total Hamiltonian iŝ whereĤ S andĤ B α are respectively the Hamiltonians of system and baths, whileĤ SB α describes the coupling between S and B α . Crucially, we assume that each interaction term couples to a distinct region on the boundary of S. We denote the corresponding interaction region by V α . Each V α is assumed to form a proper subset of the network S. See Fig. 6 for an illustration.
Suppose that there exists a conserved quantitŷ such that [Ĥ,X] = 0, which is also locally conserved in the sense that [Ĥ S ,X S ] = [Ĥ B α ,X B α ] = 0. It follows that For simplicity, we assume thatX S is a one-body observable of the formX where the local "charge"x k has support only on site k of S. This scenario could describe, for example, a lattice system of conserved particles or a spin network with conserved magnetisation.
Consider now a particular site on the boundary of the system, which lies in the region V α that is coupled to bath B α . On this site, the equation of motion forx k takes the form of a continuity equation: The current operators de ned byĴ X B α →k = i[Ĥ SB α ,x k ] and J X S→k = i[Ĥ S ,x k ] describe the ow of charge from B α to site k and from the other sites of S to site k, respectively. In particular, the operatorsĴ X S→k correspond to the internal current observables denoted schematically byĴ S in Section II. By definition, these observables have zero mean unless the quantum state has some coherence in the eigenbasis ofĤ S . This follows because, by the cyclic invariance of the trace, which clearly vanishes if [Ĥ S ,ρ S (t)] = 0. However, these expectation values do not vanish whenever the NESS supports a current due to the in ux of charge from the external reservoirs. Indeed, summing the expectation value of Eq. A5 over all sites in V α , we obtain Here, we have used Eq. (A3) and the fact that expectation values of system observables are time-independent in the NESS.
Identifying J X α = − ∂ tXB α ∞ as the current entering the system from B α , the assertion (3) is con rmed. We conclude that boundary-driven currents imply energy-eigenbasis coherence in the NESS.
For clarity, it is useful to consider the speci c case of a lattice with two-body interactions described by the generic HamiltonianĤ whereĥ kl =ĥ lk =ĥ † kl is an operator with support on sites k and l. In this case, we havê where A k denotes the the neighbourhood of site k, i.e. all sites l such thatĥ kl 0, whileĴ X l→k = i[ĥ lk ,x k ] denotes the current owing from site l to site k. According to Eq. (A6), eigenstates |E ofĤ S obey the constraint N This constraint allows for non-zero net currents around loops in the network (see, for example, Ref. 87), such that the local charge distribution is stationary (red arrows in Fig. 6). However, currents induced by external sources violate the constraint (A10) at the boundaries of the system and therefore some coherence in the energy eigenbasis is necessary to represent them (blue arrows in Fig. 6).
The simplest example of this principle is that of a onedimensional chain with open (i.e. non-periodic) boundary conditions. Since the boundary sites of such a chain have only a single neighbour, Eq. (A10) reduces to the identity That is, the energy eigenstates of a chain do not support conserved currents.
In conclusion, the existence of energy-eigenbasis coherence in a boundary-driven NESS of an extended open system is a rather general consequence of conservation laws and the locality of Hamiltonian interactions. We note that one may generalise the above argument to account for more general conserved quantities, such as two-body operators of the form X S = 1 dt e −iω(t−t ) K j (ω , t ).
(B8) (a) (b) Figure 7. Integration contours used to invert the Laplace trans-formG j (z), assuming that the inequality (34) holds. (a) A Bromwich contour lies to the right of the branch cut (thick red line) and extends parallel to the entire imaginary axis. (b) A closed integration contour that avoids enclosing any singularities. Taking the limits R → ∞ and → 0, the integrals along C 3 , C 4 and C 5 sum to zero, leaving only the contributions from C 1 and C 2 .
By inspection of Eq. (B3) one sees that G(t) is indeed given by Eq. (23). Expressions (29) and (31) for the correlation matrix C(t) and double occupancy D(t) follow directly upon plugging Eq. (B3) into the de nitions (21) and (22). Note that the above equations hold for an arbitrary spectral density J(ω), which determines the propagator via Eqs. (B5) and (B6).
Let us now compute the propagator for the speci c problem at hand, with spectral density (19). The integral (B6) evaluates toW where the positive (negative) sign for the square root is used for Re z > 0 ( Re z < 0). This implies that any poles ofG j (z) lie on the imaginary axis, which would correspond to undamped oscillations in the time domain. In order to avoid this behaviour, we assume that inequality (34) holds so thatG j (z) has no poles. As a result,G j (z) is holomorphic everywhere in the complex z plane except for along a cut connecting the branch points at z = ±iΩ. Therefore, G j (t) → 0 as t → ∞ and S relaxes to a unique stationary state (see Appendix E). The propagator is now obtained by carrying out the inverse Laplace transform, for t > 0, Here, C 1 denotes a Bromwich contour, while C 2 is a closed contour encircling the branch cut in the clockwise sense. The two integration contours can be shown to give equal and opposite contributions using Cauchy's theorem (see Fig. 7). After a change of variables to ω = iz, the integral along C 2 takes the form (24), with ϕ j (ω) given by Eq. (35). It remains to justify Eqs. (32) and (37). By taking the Laplace transform of Eq. (24), we obtaiñ We therefore deduce that, in general, Using the rst line, we obtain Eq. (32) from Eq. (29), while comparison of the second line with Eqs. (B5) and (B9) yields Eq. (37) for the uniform-chain environment model.
Appendix C: Thermalisation at finite system-bath coupling In this appendix, we prove Eq. (39). That is, when β α = β, µ α = µ and f α (ω) = f (ω), the open system equilibrates to a stateρ ∞ S given by the reduction of a global Gibbs distribution. The following arguments can also be used to deduce Eq. (60) from Eq. (59).
Since the marginal of a Gaussian state is itself Gaussian, 88 it su ces to check that Eq. (39) yields the correct correlation matrix (38). Due to the symmetry [Ĥ,N j ] = 0 of the global Hamiltonian, it follows that its equilibrium correlation matrix is diagonal, i.e. Tr[e −β(Ĥ−µN)d † jd k ] ∝ δ jk . We explicitly compute the diagonal elements using a basis transformation to the eigenmodes ofĤ, given by 18 d j = n 0|d j |n ê n , where |n =ê † n |0 is a single-particle eigenstate ofĤ [cf. Eq. (25)] and the ladder operatorsê n diagonalise the total Hamiltonian asĤ = n ε nê † nên . Substituting the above and using Eq. (25), straightforward manipulations lead to in agreement with Eq. (39), which completes the proof.

Appendix D: Computing the currents
Due to the homogeneity of the steady-state currents, the particle current may be computed from the expectation value of either the intrasystem current (42) or the boundary current (40). The expected value of the intrasystem current is simply Ĵ P S t = g Im [C 12 (t)]. In the limit t → ∞, the particle current can thus be found directly from Eq. (32): To compute the expectation value of Eqs. (40) and (41), we make use of Eqs. (B3) and (B4) and obtain, for example where K(ω, t) = diag[K 1 (ω, t), K 2 (ω, t)] and I(ω , ω, t) = diag[I 1 (ω , ω, t), I 2 (ω , ω, t)]. In the limit t → ∞, we nd that The above equations hold for arbitrary spectral densities. The explicit expressions (43)-(45) for the Newns spectral density (19) can be derived by straightforward algebra with the help of Eq. (37). while the error term is given by A trivial rearrangement of terms leads to The rst term above is the EPA contribution (53). The second term is a small relative correction to the leading-order result, which is clearly of rst order in Γ/Ω and decays exponentially in time. Finally, the remaining absolute error |err j (t)| can be bounded as where I s (x) is a modi ed Bessel function of the rst kind. The function u(τ) is positive, monotonically decreasing and vanishes as τ → ∞ with the limiting behaviour lim τ→∞ τ 3/2 u(τ) = (2π) −1/2 . Hence, by replacing u(Ωt) in Eq. (E5) by u(0) ≈ 0.54 one obtains a time-independent upper bound on |err j (t)| that is of order O(Γ/Ω), as claimed. Note, however, that the magnitude of the error also depends on E j /Ω. In particular, the error is larger when E j is closer to the band edge at ±Ω.

Appendix F: Entropy-production inequality
This appendix demonstrates the entropy-production inequality (64), starting from the Spohn inequality (63). We rst demonstrate that L int does not generate a positive evolution, i.e. the map e L int t is not positive. We write L int in diagonal form where λ 1 = −λ 2 = −λ 3 = λ 4 = |Λ + 12 |,L 1 =L † 2 = for every orthonormal pair of states |a and |b . 89 However, it is straightforward to nd a counterexample. Consider, for instance, the choice |b =L † 3 |0 and |a = |0 . Then a|L j |b = δ j3 and, since λ 3 < 0, inequality (F2) does not hold. Hence, the Spohn inequality cannot be directly applied to L int . 76 Nevertheless, the Lindblad generators L α do satisfy the Spohn inequality (63). Hence, takingρ S (t) =ρ ∞ S and summing over the baths, we obtain Let us examine the rst term on the LHS of inequality (F3) more closely. This term represents the negative time derivative of the system's von Neumann entropy generated by L int . Figure 9. Depiction of the NESS restricted to the single-particle subspace. The state, described by a Bloch vector v, is displaced by a ow u generated by L int , which always points towards the surface of the Bloch sphere and thus decreases the entropy.
This term is non-negative, which can be seen by the following geometrical argument, illustrated in Fig. 9. According to Eq. (68), L int changes only the coherences. Hence, it su ces to consider the action of L int in the two-dimensional subspace spanned by the single-particle states |E j =d † j |0 . The restriction of the density matrixρ ∞ S to the single-particle subspace is represented by a Bloch vector, v, which we choose to lie perpendicular to the y axis. Speci cally, the projection of v along the z axis is given by 1 2 + C 11 − C 22 , while the projection along the x axis is C 12 , where all quantities are evaluated in the limit t → ∞.
Consider now the ow generated by L int over a small time interval δt, i.e. the vector u such that the state e L int δtρ∞ S corresponds to the shifted Bloch vector v + uδt + O(δt 2 ). According to Eq. (68), u lies in the plane perpendicular to the z axis and points in a direction determined by the complex argument of Λ + 12 . Using Eq. (48), we deduce that Λ + 12 = (Tr[Γ]/2 − i∆)C 12 . Thus, the angle subtended by u from the x axis is φ = arctan(2∆/Tr[Γ]), such that 0 ≤ φ ≤ π/2. It follows that the ow generated by L int never moves the Bloch vector away the surface of the Bloch sphere, i.e. the length of v does not decrease. Since the von Neumann entropy is a monotonically decreasing function of the length of v, we conclude that it is a non-increasing function under this ow. Now let us demonstrate that the derivative of the von Neumann entropy generated by L int behaves as o(Γ) in the limit Γ → 0. For this it is convenient to use an explicit representation of the Gaussian NESS, where P is a positive semi-de nite matrix that satis es C = [e P + 1] −1 . Then, making use of Eq. (68), we nd the ex-plicit representation Tr lnρ ∞ S L intρ ∞ S = −2 Re Λ + 12 P 21 .