Non-additive dissipation in open quantum networks out of equilibrium

An improved understanding of the dynamics of open quantum networks is desirable for research fields including quantum thermodynamics [1–4], quantum biology [5], mesoscopic electronics [6] and the theory of non-equilibrium phase transitions [7–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. Additional, uncontrolled approximations are typically required in order to obtain a completely positive evolution, such as the secular [15], singular-coupling-limit [16] or wide-band-limit [17, 18] approximations. However, in composite open systems, distinct approximations may lead to different 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–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-infinite, uniform tight-binding 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 flowing through a serial double quantum dot [39] or cold fermionic atoms confined by an optical lattice [40, 41].

Note that other authors have already derived and extensively studied exact time-local master equations describing general quadratic Fermi systems [42–45]. In contrast to these previous approaches, the simplicity and symmetry of our specific set-up enables compact analytical solutions to be obtained using elementary methods and without needing the Born–Markov, secular, singular-coupling-limit or wide-band-limit approximations. Instead, we use an alternative, perturbative approximation scheme, valid when the system-environment coupling is much smaller than the bandwidth of the reservoir vacuum noise, and explicitly confirm its accuracy using the exact solution. This simplification allows us to clearly identify and distinguish the different physical processes governing the dynamics.

Our analysis challenges a central tenet of standard open-systems 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 [46]). Instead, we distinguish three contributions to the generator of asymptotic time evolution. Two of these describe the individual thermalising effect 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 reflects the conserved fermion current flowing in the NESS.

These findings connect several previous theoretical studies on composite open quantum systems. In particular, a significant 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–21, 24, 47–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–56].

We extend these results by demonstrating that interference between different 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 Lindblad form and thus underlies the occurrence of asymptotic non-Markovianity found by Ribeiro et al [45], where the long-time dynamics cannot be described by a Markovian master equation even when the open system has lost all memory of its initial state [57–59]. Moreover, we show that our non-additive master equation interpolates between the local and global Lindblad equations and recovers them in different limits. Since additivity is a necessary consequence of the BMA [56], this implies the failure of the BMA for describing steady-state transport away from these limits. Nevertheless, we find that certain observables—such as the steady-state currents flowing 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 2 introduces some preliminary concepts and overviews our main results. The microscopic model that forms the core of this work is defined and exactly solved in section 3. We derive the exact master equation describing the system and develop a weak-coupling approximation scheme in section 4. Non-additivity of the asymptotic dynamics is explored in section 5. We discuss our results and conclude in section 6.

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 figure 1. Each reservoir B_α is characterised by a temperature ${T}_{\alpha }=1/{\beta }_{\alpha }$ and chemical potential ${\mu }_{\alpha }$ (we work in units where k_B = 1 and ${\hslash }=1$). Let ${\hat{H}}_{S}$ and ${\hat{N}}_{S}$ respectively be the Hamiltonian and particle number operators on $S$, where $[{\hat{H}}_{S},{\hat{N}}_{S}]=0$ unless ${\mu }_{\alpha }=0$. We assume for simplicity that ${\hat{H}}_{S}$ has a non-degenerate spectrum. The quantum state of the network at time t is denoted ${\hat{\rho }}_{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 ${\hat{\rho }}_{S}^{\infty }={\mathrm{lim}}_{t\to \infty }{\hat{\rho }}_{S}(t)$. We now summarise the behaviour that we expect in general.

Zoom In Zoom Out Reset image size

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. 

$$\begin{eqnarray}&&{\hat{\rho }}_{S}^{\infty }=\displaystyle \frac{{{\rm{e}}}^{-\beta ({\hat{H}}_{S}-\mu {\hat{N}}_{S})}}{{ \mathcal Z }\,(\beta ,{\hat{H}}_{S},\mu ,{\hat{N}}_{S})}\quad \mathrm{if}\ {\beta }_{\alpha }=\beta \ \mathrm{and}\ {\mu }_{\alpha }=\mu .\end{eqnarray} \tag{ 1 }$$

Here, ${ \mathcal Z }\,(\beta ,\hat{H},\mu ,\hat{N})=\mathrm{Tr}[{{\rm{e}}}^{-\beta (\hat{H}-\mu \hat{N})}]$ is the partition function.

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 flowing into $S$ from B_α. These energy and particle currents are respectively denoted ${J}_{\alpha }^{E}$ and ${J}_{\alpha }^{P}$, while ${J}_{\alpha }^{Q}={J}_{\alpha }^{E}-{\mu }_{\alpha }{J}_{\alpha }^{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:

$$\begin{eqnarray}&&-\displaystyle \sum _{\alpha }{\beta }_{\alpha }{J}_{\alpha }^{Q}\geqslant 0.\end{eqnarray} \tag{ 2 }$$

Assuming that different reservoirs couple to different 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 flowing within the system. Therefore, there exist one or more current observables on $S$, here denoted schematically by ${\hat{J}}_{S}$, having non-zero expectation value in the NESS:

$$\begin{eqnarray}&&{\langle {\hat{J}}_{S}\rangle }_{\infty }\ne 0\quad \mathrm{unless}\quad {J}_{\alpha }^{E,P}=0.\end{eqnarray} \tag{ 3 }$$

(We denote expectation values at time t by $\langle \bullet {\rangle }_{t}$, with $\langle \bullet {\rangle }_{\infty }$ the limiting value as $t\to \infty$.) The existence of such internal currents in a boundary-driven system implies that the NESS must exhibit coherence in the eigenbasis of ${\hat{H}}_{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 ${\hat{H}}_{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 figure 1 is the quantum master equation

$$\begin{eqnarray}&&{\partial }_{t}{\hat{\rho }}_{S}(t)=-{\rm{i}}[{\hat{H}}_{S},{\hat{\rho }}_{S}(t)]+{ \mathcal L }{\hat{\rho }}_{S}(t),\end{eqnarray} \tag{ 4 }$$

valid for times t much greater than the environment memory time. In order for equation (4) to generate a completely positive and trace-preserving (CPTP) evolution for any state ${\hat{\rho }}_{S}(t)$, the dissipator ${ \mathcal L }$ must be in Lindblad form [13, 14] ${ \mathcal L }={\sum }_{j}{\gamma }_{j}{ \mathcal D }[{\hat{L}}_{j}]$, where ${ \mathcal D }[\hat{L}]{\hat{\rho }}_{S}=\hat{L}{\hat{\rho }}_{S}{\hat{L}}^{\dagger }-\tfrac{1}{2}\{{\hat{L}}^{\dagger }\hat{L},{\hat{\rho }}_{S}\}$ and ${\hat{L}}_{j}$ is a jump operator describing an incoherent transition occurring at a rate ${\gamma }_{j}$.

Since the reservoirs are assumed to initially be statistically independent, the standard construction of the dissipator is a sum of generators ${{ \mathcal L }}_{\alpha }$ representing each bath B_α, i.e.

$$\begin{eqnarray}&&{ \mathcal L }=\displaystyle \sum _{\alpha }{{ \mathcal L }}_{\alpha }.\end{eqnarray} \tag{ 5 }$$

This constitutes our definition of additivity, as studied previously in [54–56]. This is distinct from the concept of additive decoherence rates explored, for example, in [60, 61]. The additivity assumption permits one to unambiguously identify the particle current ${J}_{\alpha }^{P}(t)$ and energy current ${J}_{\alpha }^{E}(t)$ entering the system from B_α as

$$\begin{eqnarray}&&{J}_{\alpha }^{P}(t)={\langle {{ \mathcal L }}_{\alpha }^{\dagger }{\hat{N}}_{S}\rangle }_{t},\quad {J}_{\alpha }^{E}(t)={\langle {{ \mathcal L }}_{\alpha }^{\dagger }{\hat{H}}_{S}\rangle }_{t}.\end{eqnarray} \tag{ 6 }$$

Here, ${{ \mathcal L }}_{\alpha }^{\dagger }$ is the adjoint generator describing the Heisenberg-picture evolution of observables, defined by $\mathrm{Tr}[\hat{B}{{ \mathcal L }}_{\alpha }^{\dagger }\hat{A}]=\mathrm{Tr}[\hat{A}{{ \mathcal L }}_{\alpha }\hat{B}]$ for arbitrary operators $\hat{A}$ and $\hat{B}$.

Regarding the specific form of the generators ${{ \mathcal L }}_{\alpha }$, various inequivalent approaches are commonly employed. These can be broadly classified into two groups according to the fixed point of each generator, i.e. the state ${\hat{r}}_{\alpha }$ satisfying ${{ \mathcal L }}_{\alpha }{\hat{r}}_{\alpha }=0$. The first is the 'global' approach, where each generator drives the entire network towards the corresponding equilibrium state, i.e. 

$$\begin{eqnarray}&&{{ \mathcal L }}_{\alpha }{\hat{r}}_{\alpha }=0\,\Longleftrightarrow \,{\hat{r}}_{\alpha }=\displaystyle \frac{{{\rm{e}}}^{-{\beta }_{\alpha }({\hat{H}}_{S}-{\mu }_{\alpha }{\hat{N}}_{S})}}{{ \mathcal Z }\,({\beta }_{\alpha },{\hat{H}}_{S},{\mu }_{\alpha },{\hat{N}}_{S})}.\end{eqnarray} \tag{ 7 }$$

Note that here we assume the fixed point ${\hat{r}}_{\alpha }$ to be unique.

Lindblad generators satisfying equation (7) can be derived from a microscopic model under the Born–Markov and secular approximations [15]. The latter approximation is assumed to be justified in the quantum-optical regime, where the separation between energy levels of ${\hat{H}}_{S}$ is much larger than their environment-induced broadening. Equation (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 ${\hat{H}}_{S}$ [15]. Adding such generators together according to equation (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 ${\hat{H}}_{S}$, which is generally inconsistent with the condition (3) [24].

Alternatively, one can use a 'local' approach, where each generator ${{ \mathcal L }}_{\alpha }$ acts non-trivially only on one part of the network ${s}_{\alpha }\subset S$, driving it towards thermal equilibrium while leaving its complement ${\bar{s}}_{\alpha }$ unaffected. That is,

$$\begin{eqnarray}&&{{ \mathcal L }}_{\alpha }{\hat{r}}_{\alpha }=0\,\Longleftrightarrow \,{\hat{r}}_{\alpha }=\displaystyle \frac{{{\rm{e}}}^{-{\beta }_{\alpha }({\hat{H}}_{{s}_{\alpha }}-{\mu }_{\alpha }{\hat{N}}_{{s}_{\alpha }})}}{{ \mathcal Z }\,({\beta }_{\alpha },{\hat{H}}_{{s}_{\alpha }},{\mu }_{\alpha },{\hat{N}}_{{s}_{\alpha }})}{\hat{O}}_{{\bar{s}}_{\alpha }},\end{eqnarray} \tag{ 8 }$$

where ${\hat{H}}_{{s}_{\alpha }}$ and ${\hat{N}}_{{s}_{\alpha }}$ are respectively the Hamiltonian and particle number operator of s_α, while ${\hat{O}}_{{\bar{s}}_{\alpha }}$ is an arbitrary density operator with support on the complement ${\bar{s}}_{\alpha }$. Clearly, ${\hat{r}}_{\alpha }$ is not unique in this case.

Generators satisfying equation (8) may either be derived microscopically using various approximations [19, 24, 47, 49], or directly postulated on phenomenological grounds [63–66]. In the local approach, ${\hat{\rho }}_{S}^{\infty }$ is not necessarily diagonal in the energy eigenbasis, and therefore provides a consistent model for the internal current dynamics as required by equation (3). However, neither equation (1) nor (2) hold, in general, leading to potential violations of thermodynamic laws [20, 21].

In what follows, we consider the simplest case of a two-site network with two independent reservoirs labelled by the index $\alpha =L,R$. We show that, in the limit of weak system-reservoir coupling, the asymptotic dynamics is governed by a master equation whose dissipator takes the form

$$\begin{eqnarray}&&{ \mathcal L }={{ \mathcal L }}_{L}+{{ \mathcal L }}_{R}+{{ \mathcal L }}_{\mathrm{int}}.\end{eqnarray} \tag{ 9 }$$

In the quantum-optical limit, ${{ \mathcal L }}_{L}$ and ${{ \mathcal L }}_{R}$ determine the currents from the baths according to equation (6) and induce thermalisation according to equation (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 system-bath interaction) that accounts for system-reservoir correlations [67]. The interference term ${{ \mathcal L }}_{\mathrm{int}}$, which appears whenever B_L and B_R are not in equilibrium with each other, generates coherence in eigenbasis of ${\hat{H}}_{S}$ as required by condition (3). In this way, all three properties (1)–(3) can be satisfied, 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 flowing between the baths. Hence, equation (5), which is justified by the initial statistical independence of the reservoirs, fails to hold as $t\to \infty$, even if the reservoirs are spectrally unstructured and coupled arbitrarily weakly to the system.

3.1. Description of the model

Let us now detail our specific set-up, which belongs to the well-known family of resonant-level transport models [68, 69]. The total Hamiltonian takes the form $\hat{H}={\hat{H}}_{S}+{\hat{H}}_{B}+{\hat{H}}_{SB}$, describing a central open system $S$ sandwiched between two fermionic particle reservoirs B_L and B_R, as illustrated in figure 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 chemical-potential gradient.

Zoom In Zoom Out Reset image size

Specifically, $S$ comprises two localised fermionic modes with Hamiltonian

$$\begin{eqnarray}&&{\hat{H}}_{S}=\displaystyle \sum _{\alpha =L,R}{h}_{\alpha }{\hat{n}}_{\alpha }-\displaystyle \frac{g}{2}({\hat{c}}_{L}^{\dagger }{\hat{c}}_{R}+{\hat{c}}_{R}^{\dagger }{\hat{c}}_{L}),\end{eqnarray} \tag{ 10 }$$

where ${\hat{c}}_{\alpha }$ annihilates a fermion on site $\alpha =L,R$ and satisfies the anti-commutation relations $\{{\hat{c}}_{\alpha },{\hat{c}}_{\alpha ^{\prime} }^{\dagger }\}={\delta }_{\alpha \alpha ^{\prime} }$ and $\{{\hat{c}}_{\alpha },{\hat{c}}_{\alpha ^{\prime} }\}=0$, while ${\hat{n}}_{\alpha }={\hat{c}}_{\alpha }^{\dagger }{\hat{c}}_{\alpha }$. We parametrise the local energies h_α by their mean $h=\tfrac{1}{2}({h}_{L}+{h}_{R})$ and detuning $\delta ={h}_{L}-{h}_{R}$, with g the tunnel coupling between the sites.

The baths are particle reservoirs described by the Hamiltonian ${\hat{H}}_{B}={\hat{H}}_{L}+{\hat{H}}_{R}$. Each Hamiltonian ${\hat{H}}_{\alpha }$ describes a uniform chain, i.e. a 1D tight-binding model on M sites with hopping amplitude ${\rm{\Omega }}\gt 0$, given explicitly by

$$\begin{eqnarray}&&{\hat{H}}_{\alpha }=-\displaystyle \frac{{\rm{\Omega }}}{2}\displaystyle \sum _{m=1}^{M-1}({\hat{A}}_{m,\alpha }^{\dagger }{\hat{A}}_{m+1,\alpha }+{\hat{A}}_{m+1,\alpha }^{\dagger }{\hat{A}}_{m,\alpha })\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&=\displaystyle \sum _{q}{\omega }_{q}{\hat{a}}_{q,\alpha }^{\dagger }{\hat{a}}_{q,\alpha }.\,\end{eqnarray} \tag{ 12 }$$

Here, ${\hat{A}}_{m,\alpha }$ annihilates a fermion localised on site m of bath B_α and satisfies $\{{\hat{A}}_{m,\alpha },{\hat{A}}_{m^{\prime} ,\alpha ^{\prime} }^{\dagger }\}={\delta }_{{mm}^{\prime} }{\delta }_{\alpha \alpha ^{\prime} }$ and $\{{\hat{A}}_{m,\alpha },{\hat{A}}_{m^{\prime} ,\alpha ^{\prime} }\}=0$. On the second line, the Hamiltonian is diagonalised by the canonical transformation

$$\begin{eqnarray}&&{\hat{a}}_{q,\alpha }=\sqrt{\displaystyle \frac{2}{M+1}}\displaystyle \sum _{m=1}^{M}\sin ({qm}){\hat{A}}_{m,\alpha }.\end{eqnarray} \tag{ 13 }$$

The ladder operators ${\hat{a}}_{q,\alpha }$ describe quasi-free fermionic modes indexed by a dimensionless wave number $q=\pi k/(M+1)$, for $k=1,2,\,\ldots ,\,M$, with dispersion relation ${\omega }_{q}=-{\rm{\Omega }}\cos (q)$.

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 $\sqrt{{\rm{\Gamma }}{\rm{\Omega }}}$, where Γ sets the overall frequency scale of the dissipative dynamics. Explicitly, the interaction Hamiltonian reads as

$$\begin{eqnarray}&&{\hat{H}}_{SB}=\displaystyle \frac{\sqrt{{\rm{\Gamma }}{\rm{\Omega }}}}{2}\displaystyle \sum _{\alpha =L,R}({\hat{c}}_{\alpha }^{\dagger }{\hat{A}}_{1,\alpha }+{\hat{A}}_{1,\alpha }^{\dagger }{\hat{c}}_{\alpha })\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\,=\displaystyle \sum _{\alpha =L,R}\displaystyle \sum _{q}({\tau }_{q}{\hat{c}}_{\alpha }^{\dagger }{\hat{a}}_{q,\alpha }+{\tau }_{q}^{* }{\hat{a}}_{q,\alpha }^{\dagger }{\hat{c}}_{\alpha }),\end{eqnarray} \tag{ 15 }$$

with tunnel couplings ${\tau }_{q}=\sqrt{{\rm{\Gamma }}{\rm{\Omega }}/(2M+2)}\sin (q)$.

It is convenient to analyse the problem in a basis that diagonalises ${\hat{H}}_{S}$. To do this, we collect the ladder operators into column vectors $\hat{{\bf{c}}}={({\hat{c}}_{L},{\hat{c}}_{R})}^{{\mathsf{T}}}$ and ${\hat{{\bf{a}}}}_{q}={({\hat{a}}_{q,L},{\hat{a}}_{q,R})}^{{\mathsf{T}}}$. We then define a new canonical set of ladder operators $\hat{{\bf{d}}}={({\hat{d}}_{1},{\hat{d}}_{2})}^{{\mathsf{T}}}={\mathsf{R}}\hat{{\bf{c}}}$ and ${\hat{{\bf{b}}}}_{q}={({\hat{b}}_{q,1},{\hat{b}}_{q,2})}^{{\mathsf{T}}}={\mathsf{R}}{\hat{{\bf{a}}}}_{q}$, related by the orthogonal rotation matrix

$$\begin{eqnarray}{\mathsf{R}}=\displaystyle \frac{1}{\sqrt{2{\rm{\Delta }}}}\left(\begin{array}{cc}\sqrt{{\rm{\Delta }}+\delta } & -\sqrt{{\rm{\Delta }}-\delta }\\ \sqrt{{\rm{\Delta }}-\delta } & \sqrt{{\rm{\Delta }}+\delta }\end{array}\right),\end{eqnarray} \tag{ 16 }$$

where ${\rm{\Delta }}=\sqrt{{g}^{2}+{\delta }^{2}}$. The Hamiltonian hence splits into two independent pieces $\hat{H}={\hat{H}}_{1}+{\hat{H}}_{2}$, with

$$\begin{eqnarray}&&{\hat{H}}_{j}={E}_{j}{\hat{d}}_{j}^{\dagger }{\hat{d}}_{j}+\displaystyle \sum _{q}({\omega }_{q}{\hat{b}}_{q,j}^{\dagger }{\hat{b}}_{q,j}+{\tau }_{q}{\hat{d}}_{j}^{\dagger }{\hat{b}}_{q,j}+{\tau }_{q}^{* }{\hat{b}}_{q,j}^{\dagger }{\hat{d}}_{j}),\end{eqnarray} \tag{ 17 }$$

where ${E}_{1}=h+\tfrac{1}{2}{\rm{\Delta }}$ and ${E}_{2}=h-\tfrac{1}{2}{\rm{\Delta }}$ are the single-particle energy eigenvalues of ${\hat{H}}_{S}$. For concreteness, we assume that ${E}_{j}\gt 0$, or equivalently that ${h}_{\alpha }\gt 0$ and $g\lt 2\sqrt{{h}_{L}{h}_{R}}$.

The effect of the reservoirs on the system is determined by the spectral density

$$\begin{eqnarray}&&{ \mathcal J }(\omega )=\displaystyle \sum _{q}| {\tau }_{q}{| }^{2}\delta (\omega -{\omega }_{q}).\end{eqnarray} \tag{ 18 }$$

This is a smooth function of ω in the limit $M\to \infty$, where the spacing between adjacent wave vectors ${\rm{\Delta }}q=\pi /(M+1)$ tends to zero and q becomes a continuous variable taking values in the first Brillouin zone $q\in [0,\pi ]$. Using the prescription ${\sum }_{q}{\rm{\Delta }}q\to \int {\rm{d}}q$, we obtain

$$\begin{eqnarray}&&{{ \mathcal J }}_{{\rm{N}}}(\omega )=\displaystyle \frac{{\rm{\Gamma }}}{2\pi }\sqrt{1-\displaystyle \frac{{\omega }^{2}}{{{\rm{\Omega }}}^{2}}}\ {\rm{\Theta }}({\rm{\Omega }}-| \omega | ),\end{eqnarray} \tag{ 19 }$$

where ${\rm{\Theta }}(x)$ is the Heaviside unit step function. We label this spectral density with a subscript ${\rm{N}}$ after Newns, who (to our knowledge) introduced it [70]. According to equation (19), each environment is characterised by a spectral bandwidth Ω, leading to a vacuum correlation time of order ${{\rm{\Omega }}}^{-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 influence.

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–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 ${{\rm{\Omega }}}^{-1}=0$ corresponds to the wide-band-limit approximation [18], which leads to a frequency-independent spectral density. In the following, we compute the solutions for finite Ω, which enables us to retain energy-dependent damping rates even for weak coupling, ${\rm{\Gamma }}\ll {\rm{\Omega }}$, since we need not assume that ${E}_{j}\ll {\rm{\Omega }}$.

3.2. Formal solution

In this section, we provide the exact solution for the open-system density matrix ${\hat{\rho }}_{S}(t)={\mathrm{Tr}}_{B}[\hat{\rho }(t)]$, where $\hat{\rho }(t)$ is the global quantum state at time t. The same solution formally applies to the general scenario depicted in figure 2(a), i.e. any pair of environments described by Hamiltonians of the form (12) and (15), corresponding to a spectral density (18). We first 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 $\hat{\rho }(0)={\hat{\rho }}_{S}(0){\hat{\rho }}_{L}{\hat{\rho }}_{R}$, with reservoir B_α initialised in the Gibbs state

$$\begin{eqnarray}&&{\hat{\rho }}_{\alpha }=\displaystyle \frac{{{\rm{e}}}^{-{\beta }_{\alpha }({\hat{H}}_{\alpha }-{\mu }_{\alpha }{\hat{N}}_{\alpha })}}{{ \mathcal Z }\,({\beta }_{\alpha },{\hat{H}}_{\alpha },{\mu }_{\alpha },{\hat{N}}_{\alpha })}.\end{eqnarray} \tag{ 20 }$$

Here, ${\hat{N}}_{\alpha }={\sum }_{q}{\hat{a}}_{q,\alpha }^{\dagger }{\hat{a}}_{q,\alpha }$ is the number of particles in B_α.

The dynamics preserves the total number of fermions due to the relation $[\hat{H},\hat{N}]=0$, where $\hat{N}={\hat{N}}_{L}+{\hat{N}}_{R}+{\hat{N}}_{S}$ and ${\hat{N}}_{S}={\hat{n}}_{L}+{\hat{n}}_{R}$. For any physical initial state satisfying $[\hat{N},\hat{\rho }(0)]=0$, we thus have $[{\hat{N}}_{S},{\hat{\rho }}_{S}(t)]={\mathrm{Tr}}_{B}[\hat{N},\hat{\rho }(t)]=0$. Therefore, ${\hat{\rho }}_{S}(t)$ is characterised by five independent real numbers, in general. Four of these are encapsulated by the Hermitian correlation matrix

$$\begin{eqnarray}&&{C}_{{jk}}(t)=\mathrm{Tr}[{\hat{\rho }}_{S}(t){\hat{d}}_{j}^{\dagger }{\hat{d}}_{k}],\end{eqnarray} \tag{ 21 }$$

which satisfies $0\leqslant {\mathsf{C}}(t)\leqslant {\mathbb{1}}$. (Here, and throughout this document, the elements of a matrix ${\mathsf{A}}$ are denoted by A_jk.) The fifth degree of freedom is the double-occupancy probability

$$\begin{eqnarray}&&D(t)=\mathrm{Tr}[{\hat{\rho }}_{S}(t){\hat{d}}_{1}^{\dagger }{\hat{d}}_{2}^{\dagger }{\hat{d}}_{2}{\hat{d}}_{1}].\end{eqnarray} \tag{ 22 }$$

We compute equations (21) and (22) by solving the equations of motion for ${\hat{d}}_{j}(t)$ and ${\hat{b}}_{q,j}(t)$ in the Laplace domain. Here, $\hat{O}(t)={{\rm{e}}}^{{\rm{i}}\hat{H}t}\hat{O}{{\rm{e}}}^{-{\rm{i}}\hat{H}t}$ denotes the Heisenberg-picture time dependence of an operator $\hat{O}$. The solution for $t\gt 0$ can be completely expressed in terms of the propagator

$$\begin{eqnarray}&&{G}_{{jk}}(t)=\langle 0| {\hat{d}}_{j}(t){\hat{d}}_{k}^{\dagger }| 0\rangle ,\end{eqnarray} \tag{ 23 }$$

where $| 0\rangle$ is the vacuum state, i.e. $\hat{N}| 0\rangle =0$. According to equation (17), the particle number in each j sector is separately conserved: $[\hat{H},{\hat{N}}_{j}]=0$, with ${\hat{N}}_{j}={\hat{d}}_{j}^{\dagger }{\hat{d}}_{j}+{\sum }_{q}{\hat{b}}_{q,j}^{\dagger }{\hat{b}}_{q,j}$. It thus follows from the definition (23) that ${G}_{{jk}}(t)={\delta }_{{jk}}{G}_{j}(t)$, with

$$\begin{eqnarray}&&{G}_{j}(t)=\int {\rm{d}}\omega \,{{\rm{e}}}^{-{\rm{i}}\omega t}{\varphi }_{j}(\omega ).\end{eqnarray} \tag{ 24 }$$

The function ${\varphi }_{j}(\omega )$ is the probability distribution of excitation energies associated with the state ${\hat{d}}_{j}^{\dagger }| 0\rangle$, i.e.

$$\begin{eqnarray}&&{\varphi }_{j}(\omega )=\displaystyle \sum _{n}| \langle n| {\hat{d}}_{j}^{\dagger }| 0\rangle {| }^{2}\delta (\omega -{\varepsilon }_{n}),\end{eqnarray} \tag{ 25 }$$

where $| n\rangle$ is a single-particle eigenstate of $\hat{H}$ with energy ${\varepsilon }_{n}$, i.e. $\hat{H}| n\rangle ={\varepsilon }_{n}| n\rangle$ (note that ${\varepsilon }_{n}$ may be negative¹ ). Therefore,

$$\begin{eqnarray}&&\int {\rm{d}}\omega \,{\varphi }_{j}(\omega )=1,\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{{\rm{\Gamma }}\to 0}{\varphi }_{j}(\omega )=\delta (\omega -{E}_{j}).\end{eqnarray} \tag{ 27 }$$

(For general spectral densities, the limit ${\rm{\Gamma }}\to 0$ here means that all tunnel couplings ${\tau }_{q}\to 0$.)

The solution can be compactly represented in terms of the matrix ${\mathsf{F}}(\omega )={\mathsf{Rf}}(\omega ){{\mathsf{R}}}^{{\mathsf{T}}}$, where ${\mathsf{f}}(\omega )=\mathrm{diag}[{f}_{L}(\omega ),{f}_{R}(\omega )]$ with ${f}_{\alpha }{(\omega )=({{\rm{e}}}^{{\beta }_{\alpha }(\omega -{\mu }_{\alpha })}+1)}^{-1}$ the Fermi–Dirac function of reservoir B_α, and the noise kernel

$$\begin{eqnarray}&&{\rm{\Phi }}(t)=\int {\rm{d}}\omega \,{{\rm{e}}}^{-{\rm{i}}\omega t}{ \mathcal J }(\omega ){\mathsf{F}}(\omega ).\end{eqnarray} \tag{ 28 }$$

With this notation, the correlation matrix reads as

$$\begin{eqnarray}&&{\mathsf{C}}(t)={{\mathsf{G}}}^{\dagger }(t){\mathsf{C}}(0){\mathsf{G}}(t)+{\mathsf{Z}}(t),\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{\mathsf{Z}}(t)={\int }_{0}^{t}{\rm{d}}t^{\prime} {\int }_{0}^{t}{\rm{d}}t^{\prime\prime} \,{{\mathsf{G}}}^{\dagger }(t^{\prime} ){\rm{\Phi }}(t^{\prime} -t^{\prime\prime} ){\mathsf{G}}(t^{\prime\prime} ).\end{eqnarray} \tag{ 30 }$$

The double-occupancy probability is

$$\begin{eqnarray}&&D(t)=\ {| \det {\mathsf{G}}(t)| }^{2}D(0)+\det {\mathsf{Z}}(t)+\displaystyle \sum _{j\ne k}[| {G}_{j}(t){| }^{2}{Z}_{{kk}}(t){C}_{{jj}}(0)-{G}_{j}^{* }(t){G}_{k}(t){Z}_{{kj}}(t){C}_{{jk}}(0)].\end{eqnarray} \tag{ 31 }$$

We see that, so long as ${\mathrm{lim}}_{t\to \infty }{G}_{j}(t)=0$, the system relaxes to a unique stationary state ${\hat{\rho }}_{S}^{\infty }$ that is independent of the initial conditions. Moreover, this ${\hat{\rho }}_{S}^{\infty }$ is Gaussian [74] because $D(\infty )=\det {\mathsf{C}}(\infty )$. The stationary state is therefore determined completely by its correlation matrix

$$\begin{eqnarray}&&{C}_{{jk}}(\infty )=\int {\rm{d}}\omega \,{ \mathcal J }(\omega ){F}_{{jk}}(\omega )[{\vartheta }_{j}(\omega )+{\rm{i}}\pi {\varphi }_{j}(\omega )][{\vartheta }_{k}(\omega )-{\rm{i}}\pi {\varphi }_{k}(\omega )].\end{eqnarray} \tag{ 32 }$$

Here, we defined the Hilbert transform of ${\varphi }_{j}(\omega )$,

$$\begin{eqnarray}&&{\vartheta }_{j}(\omega )={P}\int {\rm{d}}\omega ^{\prime} \,\displaystyle \frac{{\varphi }_{j}(\omega ^{\prime} )}{\omega -\omega ^{\prime} },\end{eqnarray} \tag{ 33 }$$

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 ${\mathsf{G}}(t)$. However, in general, computing the propagator is a challenging problem.

3.3. Propagator for uniform-chain environments

In order to find 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

$$\begin{eqnarray}&&\sqrt{{E}_{j}^{2}+{\rm{\Gamma }}{\rm{\Omega }}}\lt {\rm{\Omega }}.\end{eqnarray} \tag{ 34 }$$

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

$$\begin{eqnarray}&&{\varphi }_{j}(\omega )=\displaystyle \frac{1}{1-{\rm{\Gamma }}/{\rm{\Omega }}}\displaystyle \frac{{{ \mathcal J }}_{{\rm{N}}}(\omega )}{{(\omega -{E}_{j}^{{\prime} })}^{2}+{{\rm{\Gamma }}}_{j}^{2}/4},\end{eqnarray} \tag{ 35 }$$

where we defined shifted energies and decay rates

$$\begin{eqnarray}&&{E}_{j}^{{\prime} }=\displaystyle \frac{1-{\rm{\Gamma }}/2{\rm{\Omega }}}{1-{\rm{\Gamma }}/{\rm{\Omega }}}{E}_{j},\qquad {{\rm{\Gamma }}}_{j}=\displaystyle \frac{\sqrt{1-{\rm{\Gamma }}/{\rm{\Omega }}-{E}_{j}^{2}/{{\rm{\Omega }}}^{2}}}{1-{\rm{\Gamma }}/{\rm{\Omega }}}{\rm{\Gamma }}.\end{eqnarray} \tag{ 36 }$$

The propagator is then evaluated using equation (24). We also prove in appendix B that

$$\begin{eqnarray}&&{\vartheta }_{j}(\omega )=\left[\left(1-\displaystyle \frac{{\rm{\Gamma }}}{2{\rm{\Omega }}}\right)\omega -{E}_{j}\right]\displaystyle \frac{{\varphi }_{j}(\omega )}{{{ \mathcal J }}_{{\rm{N}}}(\omega )}.\end{eqnarray} \tag{ 37 }$$

3.4. 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. ${\beta }_{\alpha }=\beta$, ${\mu }_{\alpha }=\mu$ and ${f}_{\alpha }(\omega )=f(\omega )$, then

$$\begin{eqnarray}&&{C}_{{jk}}(\infty )={\delta }_{{jk}}\int {\rm{d}}\omega \,{\varphi }_{j}(\omega )f(\omega ),\end{eqnarray} \tag{ 38 }$$

which differs from the strict equilibrium value $f({E}_{j})$ due to the finite width of the energy distribution ${\varphi }_{j}(\omega )$. Using equations (25) and (38), we show in appendix C that, for equilibrium baths, ${\hat{\rho }}_{S}^{\infty }$ is the reduced state of the global Gibbs ensemble

$$\begin{eqnarray}&&{\hat{\rho }}_{S}^{\infty }={\mathrm{Tr}}_{B}\left[\displaystyle \frac{{{\rm{e}}}^{-\beta (\hat{H}-\mu \hat{N})}}{{ \mathcal Z }\,(\beta ,\hat{H},\mu ,\hat{N})}\right].\end{eqnarray} \tag{ 39 }$$

We thus infer that the energy-level broadening represented by equation (38) arises from the system-bath correlations reflected in equation (39), both being related to a finite system-reservoir coupling Γ. With the help of equation (27), we recover thermalisation in the strict weak-coupling sense (1) in the limit ${\rm{\Gamma }}\to 0$.

Out of equilibrium, particle and energy currents flow from each bath B_α into $S$. These currents are defined respectively as ${\hat{J}}_{\alpha }^{P}=-{\rm{i}}[\hat{H},{\hat{N}}_{\alpha }]$ and ${\hat{J}}_{\alpha }^{E}=-{\rm{i}}[\hat{H},{\hat{H}}_{\alpha }]$. Fermion conservation demands that a corresponding particle current ${\hat{J}}_{S}^{P}$ flows between the two sites of the system, which is defined to satisfy the continuity equations, e.g. ${\partial }_{t}{\hat{n}}_{L}(t)={\hat{J}}_{L}^{P}(t)-{\hat{J}}_{S}^{P}(t)$. Explicitly, we have

$$\begin{eqnarray}&&{\hat{J}}_{\alpha }^{P}={\rm{i}}\displaystyle \sum _{q}({\tau }_{q}^{* }{\hat{a}}_{q,\alpha }^{\dagger }{\hat{c}}_{\alpha }-{\tau }_{q}{\hat{c}}_{\alpha }^{\dagger }{\hat{a}}_{q,\alpha }),\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&{\hat{J}}_{\alpha }^{E}={\rm{i}}\displaystyle \sum _{q}{\omega }_{q}({\tau }_{q}^{* }{\hat{a}}_{q,\alpha }^{\dagger }{\hat{c}}_{\alpha }-{\tau }_{q}{\hat{c}}_{\alpha }^{\dagger }{\hat{a}}_{q,\alpha }),\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}&&{\hat{J}}_{S}^{P}=\displaystyle \frac{g}{2{\rm{i}}}({\hat{c}}_{L}^{\dagger }{\hat{c}}_{R}-{\hat{c}}_{R}^{\dagger }{\hat{c}}_{L}).\end{eqnarray} \tag{ 42 }$$

Note that the mean intra-system current $\langle {\hat{J}}_{S}^{P}{\rangle }_{t}=g\,\mathrm{Im}\,[{C}_{12}(t)]$ vanishes identically if $[{\hat{H}}_{S},{\hat{\rho }}_{S}(t)]=0$.

It follows from the definitions that, in the stationary state, the particle currents are homogeneous throughout the system, i.e. ${J}_{L}^{P}={J}_{S}^{P}=-{J}_{R}^{P}$, where ${J}_{\sigma }^{P}=\langle {\hat{J}}_{\sigma }^{P}{\rangle }_{\infty }$, for $\sigma =S,L,R$. Likewise, the asymptotic mean energy currents ${J}_{\alpha }^{E}=\langle {\hat{J}}_{\alpha }^{E}{\rangle }_{\infty }$ satisfy ${J}_{L}^{E}=-{J}_{R}^{E}$.

The asymptotic currents are calculated in appendix D. For the Newns spectral density (19), we obtain the standard Landauer formulae [6]

$$\begin{eqnarray}&&{J}_{L}^{P}=\int {\rm{d}}\omega \,{ \mathcal T }(\omega )[{f}_{L}(\omega )-{f}_{R}(\omega )],\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&{J}_{L}^{E}=\int {\rm{d}}\omega \,\omega { \mathcal T }(\omega )[{f}_{L}(\omega )-{f}_{R}(\omega )],\end{eqnarray} \tag{ 44 }$$

with the transmission function

$$\begin{eqnarray}&&{ \mathcal T }(\omega )=\displaystyle \frac{\pi {g}^{2}}{2}{\varphi }_{1}(\omega ){\varphi }_{2}(\omega ).\end{eqnarray} \tag{ 45 }$$

The transmission probability is proportional to the overlap between the energy distributions of the two eigenmodes of ${\hat{H}}_{S}$ (see equation (25)). Note that the Landauer formulae imply the second law (2) for any transmission function ${ \mathcal T }(\omega )\geqslant 0$ [4].

We plot the transmission function in figure 3. For ${\rm{\Gamma }}\ll {\rm{\Delta }}$, ${ \mathcal T }(\omega )$ is a bimodal distribution that is well approximated by

$$\begin{eqnarray}&&{ \mathcal T }(\omega )\approx \displaystyle \frac{{g}^{2}}{4{{\rm{\Delta }}}^{2}}\displaystyle \sum _{j=1}^{2}{{\rm{\Gamma }}}_{j}{\varphi }_{j}(\omega ).\end{eqnarray} \tag{ 46 }$$

As Γ is increased, the two peaks at $\omega ={E}_{\mathrm{1,2}}$ broaden and ultimately merge into a single maximum for ${\rm{\Gamma }}\gg {\rm{\Delta }}$.

Zoom In Zoom Out Reset image size

4.1. Exact master equation

In order to analyse the non-additive properties of the dynamics, we first write equations (29) and (31) in differential form, corresponding to an exact time-local master equation for the system density operator (see [43, 44] for alternative derivations). All equations presented in this section hold for an arbitrary spectral density (18).

Assuming that ${\mathsf{G}}(t)$ is non-singular, we define Hermitian matrices ${\mathsf{H}}(t)$ and ${\rm{\Gamma }}(t)$ such that ${\rm{i}}{\mathsf{H}}(t)+\tfrac{1}{2}{\rm{\Gamma }}(t)=-{{\mathsf{G}}}^{-1}{\partial }_{t}{\mathsf{G}}$. We also define rate matrices ${{\rm{\Lambda }}}^{\pm }(t)$ by

$$\begin{eqnarray}{{\rm{\Lambda }}}^{+} & = & {\partial }_{t}{\mathsf{Z}}-{\rm{i}}[{\mathsf{H}},{\mathsf{Z}}]+\tfrac{1}{2}\{{\rm{\Gamma }},{\mathsf{Z}}\},\\ {{\rm{\Lambda }}}^{-} & = & {\rm{\Gamma }}-{{\rm{\Lambda }}}^{+},\end{eqnarray} \tag{ 47 }$$

where time arguments are suppressed. The matrix elements of ${{\rm{\Lambda }}}^{\pm }$ correspond to the gain and loss rate coefficients appearing in the master equation, as will be seen shortly. Indeed, upon differentiating equations (29) and (31), some tedious algebra reveals that

$$\begin{eqnarray}&&{\partial }_{t}{\mathsf{C}}={\rm{i}}[{\mathsf{H}},{\mathsf{C}}]-\tfrac{1}{2}\{{\rm{\Gamma }},{\mathsf{C}}\}+{{\rm{\Lambda }}}^{+},\end{eqnarray} \tag{ 48 }$$

$$\begin{eqnarray}&&{\partial }_{t}D=-\mathrm{Tr}[{\rm{\Gamma }}]D+\displaystyle \sum _{j\ne k}({{\rm{\Lambda }}}_{{jj}}^{+}{C}_{{kk}}-{{\rm{\Lambda }}}_{{jk}}^{+}{C}_{{kj}}).\end{eqnarray} \tag{ 49 }$$

Direct comparison confirms that these equations of motion are equivalent to the master equation

$$\begin{eqnarray}&&{\partial }_{t}{\hat{\rho }}_{S}(t)=-{\rm{i}}[{\hat{H}}_{S}^{{\prime} }(t),{\hat{\rho }}_{S}(t)]+{ \mathcal L }(t){\hat{\rho }}_{S}(t),\end{eqnarray} \tag{ 50 }$$

where ${\hat{H}}_{S}^{{\prime} }(t)={\hat{{\bf{d}}}}^{\dagger }{{\mathsf{H}}}^{{\mathsf{T}}}\hat{{\bf{d}}}$ and we defined the dissipator

$$\begin{eqnarray}&&{ \mathcal L }(t){\hat{\rho }}_{S}=\displaystyle \sum _{j,k=1}^{2}{{\rm{\Lambda }}}_{{jk}}^{-}(t)\left({\hat{d}}_{j}{\hat{\rho }}_{S}{\hat{d}}_{k}^{\dagger }-\tfrac{1}{2}\{{\hat{d}}_{k}^{\dagger }{\hat{d}}_{j},{\hat{\rho }}_{S}\}\right)+\displaystyle \sum _{j,k=1}^{2}{{\rm{\Lambda }}}_{{jk}}^{+}(t)\left({\hat{d}}_{k}^{\dagger }{\hat{\rho }}_{S}{\hat{d}}_{j}-\tfrac{1}{2}\{{\hat{d}}_{j}{\hat{d}}_{k}^{\dagger },{\hat{\rho }}_{S}\}\right).\end{eqnarray} \tag{ 51 }$$

After diagonalising the rate matrices ${{\rm{\Lambda }}}^{\pm }(t)$ by a unitary rotation ${[{{\mathsf{U}}}^{\pm }(t)]}^{\dagger }{{\rm{\Lambda }}}^{\pm }(t){{\mathsf{U}}}^{\pm }(t)=\mathrm{diag}[{\lambda }_{1}^{\pm }(t),{\lambda }_{2}^{\pm }(t)]$, we cast the dissipator into canonical form [57]

$$\begin{eqnarray}&&{ \mathcal L }(t)=\displaystyle \sum _{j=1}^{2}\displaystyle \sum _{s=\pm }{\lambda }_{j}^{s}(t){ \mathcal D }[{\hat{L}}_{j}^{s}(t)].\end{eqnarray} \tag{ 52 }$$

Here, the time-dependent jump operators are defined as ${\hat{L}}_{j}^{-}(t)={\sum }_{k}{U}_{{jk}}^{-}(t){\hat{d}}_{k}$ and ${\hat{L}}_{j}^{+}(t)={\sum }_{k}{[{U}_{{jk}}^{+}(t)]}^{* }{\hat{d}}_{k}^{\dagger }$. These describe loss and gain of excitations from and into modes determined by the eigenbases of the matrices ${{\rm{\Lambda }}}^{\pm }(t)$. Note that only if $[{{\rm{\Lambda }}}^{+}(t),{{\rm{\Lambda }}}^{-}(t)]=0$ can we choose ${\hat{L}}_{j}^{+}(t)={[{\hat{L}}_{j}^{-}(t)]}^{\dagger }$, i.e. the loss and gain modes differ, in general.

4.2. Exponential-propagator approximation (EPA)

Throughout the rest of the paper, we specialise to the Newns spectral density (19) and assume that ${\rm{\Gamma }}\ll {\rm{\Omega }}$. Although the latter condition is usually deemed necessary for the BMA to hold, we shall see that it is by no means sufficient.

In order to simplify the subsequent discussion, we introduce an approximation scheme where terms of order $O({\rm{\Gamma }}/{\rm{\Omega }})$ are neglected. This amounts to the replacement

$$\begin{eqnarray}&&{G}_{j}(t)\approx {{\rm{e}}}^{-{\rm{i}}{E}_{j}^{{\prime} }t-{{\rm{\Gamma }}}_{j}t/2}.\end{eqnarray} \tag{ 53 }$$

For the sake of clarity and concision, we henceforth refer to this as the EPA. The EPA is justified 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 ${\rm{\Gamma }}/{\rm{\Omega }}$ and decays to zero as $t\to \infty$. To the same order of approximation, we write ${E}_{j}^{{\prime} }\approx {E}_{j}$ and ${{\rm{\Gamma }}}_{j}=2\pi {{ \mathcal J }}_{{\rm{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 sufficiently small ${\rm{\Gamma }}/{\rm{\Omega }}$, 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 figure 4 for a few example parameters.

Zoom In Zoom Out Reset image size

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 wide-band limit ${\rm{\Omega }}\to \infty$. 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 ${\mathsf{H}}=\mathrm{diag}[{E}_{1},{E}_{2}]$, i.e. ${{\hat{H}}_{S}^{{\prime} }}_{}(t)={\hat{H}}_{S}$, ${\rm{\Gamma }}=\mathrm{diag}[{{\rm{\Gamma }}}_{1},{{\rm{\Gamma }}}_{2}]$ and

$$\begin{eqnarray}&&{{\rm{\Lambda }}}^{+}(t)={\int }_{0}^{t}{\rm{d}}t^{\prime} \,{{\mathsf{G}}}^{\dagger }(t^{\prime} ){\rm{\Phi }}(t^{\prime} )+{\rm{h}}.{\rm{c}}.\end{eqnarray} \tag{ 54 }$$

The loss and gain modes defined by equation (52) can be identified by inspection of equation (54) in two special cases. First, assume that the two baths are in thermal equilibrium with each other, so that ${f}_{L}(\omega )={f}_{R}(\omega )=f(\omega )$. In that case, ${\mathsf{F}}(\omega )=f(\omega ){\mathbb{1}}$ is proportional to the identity, and ${{\rm{\Lambda }}}^{\pm }(t)$ are both diagonal. Hence, ${\hat{L}}_{j}^{-}={({\hat{L}}_{j}^{+})}^{\dagger }={\hat{d}}_{j}$, i.e. excitations are pumped into eigenmodes of ${\hat{H}}_{S}$, driving the system towards a state that is diagonal in the energy eigenbasis of $S$. If the system-environment coupling is sufficiently weak, this ensures the proper thermalisation behaviour (see section 5.1).

Second, consider the white-noise limit, with ${\rm{\Omega }}\to \infty$ and ${\mathsf{F}}(\omega )={\mathsf{F}}$ a constant matrix. We then have that ${\rm{\Phi }}(t)={\rm{\Gamma }}{\mathsf{F}}\delta (t)$ and ${\rm{\Gamma }}={\rm{\Gamma }}{\mathbb{1}}$, and since ${\mathsf{G}}(0)={\mathbb{1}}$ it follows that ${{\rm{\Lambda }}}^{+}(t)={\rm{\Gamma }}{\mathsf{F}}$ and ${{\rm{\Lambda }}}^{-}(t)={\rm{\Gamma }}({\mathbb{1}}-{\mathsf{F}})$. Therefore, the rate matrices ${{\rm{\Lambda }}}^{\pm }$ can be diagonalised by the rotation ${{\mathsf{R}}}^{{\mathsf{T}}}$, leading to Lindblad operators ${\hat{L}}_{j}^{-}={({\hat{L}}_{j}^{+})}^{\dagger }$ with ${\hat{L}}_{1,2}^{-}={\hat{c}}_{L,R}$, i.e. energy is pumped into localised modes. Of course, this only holds exactly in the unphysical scenario of infinite energy density in the baths, corresponding to ${\beta }_{j}\to 0$ or ${\beta }_{j},| {\mu }_{j}| \to \infty$. Nevertheless, the local master equation may be an excellent approximation for sufficiently large chemical-potential bias or temperature.

According to equation (54), the rate matrices are determined by the product of the propagator ${\mathsf{G}}(t)$—which is diagonal in the energy eigenbasis—and the noise kernel ${\rm{\Phi }}(t)$—which is diagonal in the site basis. Hence, in general, the gain and loss modes correspond neither to the energy eigenbasis ${\hat{d}}_{j}$ nor the site basis ${\hat{c}}_{\alpha }$, but instead lie somewhere in between. In particular, whenever the baths are not in equilibrium with each other, the rate matrices ${{\rm{\Lambda }}}^{\pm }$ are non-diagonal and the dissipation generates some coherence in the eigenbasis of ${\hat{H}}_{S}$. This can be seen, for example, in figure 4, where the coherence can be comparable in magnitude to the populations and has a large imaginary part, reflecting the current flowing through the system.

5.1. Non-additivity in the energy eigenbasis

Now we demonstrate that the asymptotic dynamics as $t\to \infty$ is not additive. We write the rate matrices in the limit simply as ${{\rm{\Lambda }}}^{\pm }={\mathrm{lim}}_{t\to \infty }{{\rm{\Lambda }}}^{\pm }(t)$. Explicitly, these have the components

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{{jk}}^{+}={\rm{i}}\int {\rm{d}}\omega \ {{ \mathcal J }}_{{\rm{N}}}(\omega ){F}_{{jk}}(\omega )\left[\displaystyle \frac{1}{{E}_{j}-\omega +{\rm{i}}{{\rm{\Gamma }}}_{j}/2}-\displaystyle \frac{1}{{E}_{k}-\omega -{\rm{i}}{{\rm{\Gamma }}}_{k}/2}\right],\end{eqnarray} \tag{ 55 }$$

and ${{\rm{\Lambda }}}_{{jk}}^{-}={{\rm{\Gamma }}}_{j}{\delta }_{{jk}}-{{\rm{\Lambda }}}_{{jk}}^{+}$. It follows that the asymptotic dissipator ${ \mathcal L }={\mathrm{lim}}_{t\to \infty }{ \mathcal L }(t)$ admits the decomposition

$$\begin{eqnarray}&&{ \mathcal L }={{ \mathcal L }}_{L}+{{ \mathcal L }}_{R}+{{ \mathcal L }}_{\mathrm{int}}.\end{eqnarray} \tag{ 56 }$$

Here, the generators ${{ \mathcal L }}_{\alpha }$ represent the thermalising effect of each individual bath $\alpha =L,R$ on the populations, while the interference term ${{ \mathcal L }}_{\mathrm{int}}$ describes the generation of coherence due to the combined effect of the non-equilibrium baths, as described below.

First, let us examine the generators ${{ \mathcal L }}_{\alpha }$, which take the form

$$\begin{eqnarray}&&{{ \mathcal L }}_{\alpha }=\displaystyle \sum _{j=1}^{2}({\gamma }_{\alpha ,j}^{-}{ \mathcal D }[{\hat{d}}_{j}]+{\gamma }_{\alpha ,j}^{+}{ \mathcal D }[{\hat{d}}_{j}^{\dagger }]).\end{eqnarray} \tag{ 57 }$$

The decay and gain rates are defined by

$$\begin{eqnarray}{\gamma }_{\alpha ,j}^{-} & = & {{\rm{\Gamma }}}_{\alpha ,j}\displaystyle \int {\rm{d}}\omega \,{\varphi }_{j}(\omega )[1-{f}_{\alpha }(\omega )],\\ {\gamma }_{\alpha ,j}^{+} & = & {{\rm{\Gamma }}}_{\alpha ,j}\displaystyle \int {\rm{d}}\omega \,{\varphi }_{j}(\omega ){f}_{\alpha }(\omega ),\end{eqnarray} \tag{ 58 }$$

where ${{\rm{\Gamma }}}_{L,j}={{\rm{\Gamma }}}_{j}{R}_{1j}^{2}$ and ${{\rm{\Gamma }}}_{R,j}={{\rm{\Gamma }}}_{j}{R}_{2j}^{2}$.

One readily verifies that the unique fixed point ${\hat{r}}_{\alpha }$ satisfying ${{ \mathcal L }}_{\alpha }{\hat{r}}_{\alpha }=0$ is Gaussian, with the correlation matrix

$$\begin{eqnarray}&&\mathrm{Tr}[{\hat{r}}_{\alpha }{\hat{d}}_{j}^{\dagger }{\hat{d}}_{k}]={\delta }_{{jk}}\int {\rm{d}}\omega \,{\varphi }_{j}(\omega ){f}_{\alpha }(\omega ).\end{eqnarray} \tag{ 59 }$$

In the quantum-optical limit where ${\rm{\Gamma }}\ll {\rm{\Delta }}$, this describes an approximately thermal distribution, up to corrections due to level broadening as discussed in section 3.4. In fact, using the arguments given in appendix C, one can show that ${\hat{r}}_{\alpha }$ is the reduction of a global Gibbs state in equilibrium with the corresponding bath, i.e.

$$\begin{eqnarray}&&{\hat{r}}_{\alpha }={\mathrm{Tr}}_{B}\left[\displaystyle \frac{{{\rm{e}}}^{-{\beta }_{\alpha }(\hat{H}-{\mu }_{\alpha }\hat{N})}}{{ \mathcal Z }\,({\beta }_{\alpha },\hat{H},{\mu }_{\alpha },\hat{N})}\right].\end{eqnarray} \tag{ 60 }$$

It follows immediately that equation (7) is recovered as ${\rm{\Gamma }}\to 0$.

Another interesting property of the generators ${{ \mathcal L }}_{\alpha }$ in the quantum-optical limit is that they accurately reproduce the steady-state currents according to equation (6). Indeed, a direct computation yields

$$\begin{eqnarray}&&{\langle {{ \mathcal L }}_{L}^{\dagger }{\hat{N}}_{S}\rangle }_{\infty }=\displaystyle \frac{{g}^{2}}{4{{\rm{\Delta }}}^{2}}\int {\rm{d}}\omega \displaystyle \sum _{j=1}^{2}{{\rm{\Gamma }}}_{j}{\varphi }_{j}(\omega )[{f}_{L}(\omega )-{f}_{R}(\omega )],\end{eqnarray} \tag{ 61 }$$

$$\begin{eqnarray}&&{\langle {{ \mathcal L }}_{L}^{\dagger }{\hat{H}}_{S}\rangle }_{\infty }=\displaystyle \frac{{g}^{2}}{4{{\rm{\Delta }}}^{2}}\int {\rm{d}}\omega \displaystyle \sum _{j=1}^{2}{E}_{j}{{\rm{\Gamma }}}_{j}{\varphi }_{j}(\omega )[{f}_{L}(\omega )-{f}_{R}(\omega )],\end{eqnarray} \tag{ 62 }$$

with ${\sum }_{\alpha }\langle {{ \mathcal L }}_{\alpha }^{\dagger }{\hat{N}}_{S}{\rangle }_{\infty }=0={\sum }_{\alpha }\langle {{ \mathcal L }}_{\alpha }^{\dagger }{\hat{H}}_{S}{\rangle }_{\infty }$. These expressions can be shown to be equivalent to the exact Landauer formulae (43) and (44) using equation (46) and writing ${E}_{j}{\varphi }_{j}(\omega )\approx \omega {\varphi }_{j}(\omega )$, which is a valid approximation in the quantum-optical regime.

Since the generators ${{ \mathcal L }}_{\alpha }$ are in Lindblad form, they obey the Spohn inequality [62]

$$\begin{eqnarray}&&\mathrm{Tr}\{[\mathrm{ln}{\hat{r}}_{\alpha }-\mathrm{ln}{\hat{\rho }}_{S}(t)]{{ \mathcal L }}_{\alpha }{\hat{\rho }}_{S}(t)\}\geqslant 0.\end{eqnarray} \tag{ 63 }$$

However, ${{ \mathcal L }}_{\mathrm{int}}$ does not, by itself, generate a positive evolution, and thus does not necessarily satisfy such an inequality [75]. We nonetheless show in appendix F that, in the weak-coupling limit, from the Spohn inequality one may recover the second law of thermodynamics in the form

$$\begin{eqnarray}&&-\mathop{\mathrm{lim}}\limits_{{\rm{\Gamma }}\to 0}\displaystyle \frac{1}{{\rm{\Gamma }}}\displaystyle \sum _{\alpha }{\beta }_{\alpha }{J}_{\alpha }^{Q}\geqslant 0,\end{eqnarray} \tag{ 64 }$$

with the heat current defined by ${J}_{\alpha }^{Q}=\langle {{ \mathcal L }}_{\alpha }^{\dagger }({\hat{H}}_{S}-{\mu }_{\alpha }{\hat{N}}_{S}){\rangle }_{\infty }$. Here, it is necessary to first divide by Γ before taking the limit in order to avoid recovering the trivial equality ${\mathrm{lim}}_{{\rm{\Gamma }}\to 0}{\sum }_{\alpha }{\beta }_{\alpha }{J}_{\alpha }^{Q}=0$ implied by ${\mathrm{lim}}_{{\rm{\Gamma }}\to 0}{J}_{\alpha }^{Q}=0$. Corrections to the left-hand side (LHS) of inequality (64) for finite Γ are quoted explicitly in appendix F.

Now we turn to the interference contribution, defined by

$$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{int}}{\hat{\rho }}_{S}=\displaystyle \sum _{j\ne k}{{\rm{\Lambda }}}_{{jk}}^{+}\left[\left({\hat{d}}_{k}^{\dagger }{\hat{\rho }}_{S}{\hat{d}}_{j}-\tfrac{1}{2}\{{\hat{d}}_{j}{\hat{d}}_{k}^{\dagger },{\hat{\rho }}_{S}\}\right)-\left({\hat{d}}_{j}{\hat{\rho }}_{S}{\hat{d}}_{k}^{\dagger }-\tfrac{1}{2}\{{\hat{d}}_{k}^{\dagger }{\hat{d}}_{j},{\hat{\rho }}_{S}\}\right)\right].\end{eqnarray} \tag{ 65 }$$

Here, ${{\rm{\Lambda }}}_{12}^{+}={({{\rm{\Lambda }}}_{21}^{+})}^{* }=\xi +{\rm{i}}\eta$, with

$$\begin{eqnarray}&&\xi =\displaystyle \frac{g}{2{\rm{\Delta }}}\int {\rm{d}}\omega \,\displaystyle \sum _{j=1}^{2}{{\rm{\Gamma }}}_{j}{\varphi }_{j}(\omega )[{f}_{L}(\omega )-{f}_{R}(\omega )],\end{eqnarray} \tag{ 66 }$$

$$\begin{eqnarray}&&\eta =\displaystyle \frac{g}{2{\rm{\Delta }}}\int {\rm{d}}\omega [({E}_{1}-\omega ){\varphi }_{1}(\omega )-({E}_{2}-\omega ){\varphi }_{2}(\omega )][{f}_{L}(\omega )-{f}_{R}(\omega )].\end{eqnarray} \tag{ 67 }$$

Therefore, ${{ \mathcal L }}_{\mathrm{int}}$ is associated with the difference in distribution functions ${f}_{\alpha }(\omega )$. In particular, ${{ \mathcal L }}_{\mathrm{int}}$ is non-negligible unless ${f}_{L}(\omega )\approx {f}_{R}(\omega )$ over the entire frequency range in which ${\varphi }_{j}(\omega )$ differs appreciably from zero.

The interference term satisfies the properties

$$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{int}}^{\dagger }{\hat{d}}_{j}^{\dagger }{\hat{d}}_{j}=0,\qquad {{ \mathcal L }}_{\mathrm{int}}^{\dagger }{\hat{d}}_{1}^{\dagger }{\hat{d}}_{2}={{\rm{\Lambda }}}_{12}^{+}.\end{eqnarray} \tag{ 68 }$$

Hence, ${{ \mathcal L }}_{\mathrm{int}}^{\dagger }$ generates coherence while leaving the populations unaffected. This is to be expected, since energy-eigenbasis coherence is an intrinsic property of the NESS of a quantum network, as discussed in section 2. It is noteworthy that, in the Schrödinger picture, ${{ \mathcal L }}_{\mathrm{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, ${{ \mathcal L }}_{\mathrm{int}}$ is generally a significant contribution even in the quantum-optical limit where the secular approximation is widely believed to be valid. We note also that ${{ \mathcal L }}_{\mathrm{int}}$ can never be written in Lindblad form. Therefore, departures from Markovian evolution can be used to detect non-additivity, as discussed in section 5.3.

5.2. 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 ${\hat{c}}_{\alpha }$ and decompose the asymptotic dissipator as

$$\begin{eqnarray}&&{ \mathcal L }={\bar{{ \mathcal L }}}_{L}+{\bar{{ \mathcal L }}}_{R}+{\bar{{ \mathcal L }}}_{{LR}}.\end{eqnarray} \tag{ 69 }$$

The local dissipators are Lindblad generators that act only on a single site, given by

$$\begin{eqnarray}&&{\bar{{ \mathcal L }}}_{\alpha }={\bar{\gamma }}_{\alpha }^{-}{ \mathcal D }[{\hat{c}}_{\alpha }]+{\bar{\gamma }}_{\alpha }^{+}{ \mathcal D }[{\hat{c}}_{\alpha }^{\dagger }],\end{eqnarray} \tag{ 70 }$$

where ${\bar{\gamma }}_{\alpha }^{\pm }={\sum }_{j=1}^{2}{\gamma }_{\alpha ,j}^{\pm }$. The remaining contribution to ${ \mathcal L }$ describes delocalised, incoherent processes acting on both sites together:

$$\begin{eqnarray}&&{\bar{{ \mathcal L }}}_{{LR}}{\hat{\rho }}_{S}=\displaystyle \sum _{\alpha \ne \alpha ^{\prime} }\left[{\bar{{\rm{\Lambda }}}}_{\alpha \alpha ^{\prime} }^{+}\left({\hat{c}}_{\alpha ^{\prime} }^{\dagger }{\hat{\rho }}_{S}{\hat{c}}_{\alpha }-\tfrac{1}{2}\{{\hat{c}}_{\alpha }{\hat{c}}_{\alpha ^{\prime} }^{\dagger },{\hat{\rho }}_{S}\}\right)+\ {\bar{{\rm{\Lambda }}}}_{\alpha \alpha ^{\prime} }^{-}\left({\hat{c}}_{\alpha }{\hat{\rho }}_{S}{\hat{c}}_{\alpha ^{\prime} }^{\dagger }-\tfrac{1}{2}\{{\hat{c}}_{\alpha ^{\prime} }^{\dagger }{\hat{c}}_{\alpha },{\hat{\rho }}_{S}\}\right)\right],\end{eqnarray} \tag{ 71 }$$

where ${\bar{{\rm{\Lambda }}}}_{{LR}}^{\pm }={({\bar{{\rm{\Lambda }}}}_{{RL}}^{\pm })}^{* }$, with

$$\begin{eqnarray}\,{\rm{Re}}\,[{\bar{{\rm{\Lambda }}}}_{{LR}}^{+}] & = & \displaystyle \frac{g}{4{\rm{\Delta }}}\displaystyle \int {\rm{d}}\omega [{{\rm{\Gamma }}}_{2}{\varphi }_{2}(\omega )-{{\rm{\Gamma }}}_{1}{\varphi }_{1}(\omega )][{f}_{L}(\omega )+{f}_{R}(\omega )],\\ \,{\rm{Im}}\,[{\bar{{\rm{\Lambda }}}}_{{LR}}^{+}] & = & \displaystyle \frac{g}{2{\rm{\Delta }}}\displaystyle \int {\rm{d}}\omega [({E}_{1}-\omega ){\varphi }_{1}(\omega )-({E}_{2}-\omega ){\varphi }_{2}(\omega )][{f}_{L}(\omega )-{f}_{R}(\omega )],\\ \,{\rm{Re}}\,[{\bar{{\rm{\Lambda }}}}_{{LR}}^{-}] & = & \displaystyle \frac{g}{4{\rm{\Delta }}}\displaystyle \int {\rm{d}}\omega [{{\rm{\Gamma }}}_{2}{\varphi }_{2}(\omega )-{{\rm{\Gamma }}}_{1}{\varphi }_{1}(\omega )][2-{f}_{L}(\omega )-{f}_{R}(\omega )],\\ \,{\rm{Im}}\,[{\bar{{\rm{\Lambda }}}}_{{LR}}^{-}] & = & -\,{\rm{Im}}\,[{\bar{{\rm{\Lambda }}}}_{{LR}}^{+}].\end{eqnarray} \tag{ 72 }$$

The cross-term ${\bar{{ \mathcal L }}}_{{LR}}$ is associated with the difference between the frequency distributions ${\varphi }_{j}(\omega )$ of the open system's energy levels. That is, ${\bar{{ \mathcal L }}}_{{LR}}$ reflects the extent to which the occupation numbers ${f}_{\alpha }(\omega )$ of reservoir states differ between the distinct frequency ranges sampled by the two distributions ${\varphi }_{j}(\omega )$. In particular, ${\bar{{ \mathcal L }}}_{{LR}}$ is negligible only in the white-noise limit, where ${\rm{\Delta }}\ll {\rm{\Omega }}$ and the distribution functions ${f}_{L,R}(\omega )$ are essentially constant over the range in which ${\varphi }_{j}(\omega )$ are non-zero.

Note that here the local generators ${\bar{{ \mathcal L }}}_{\alpha }$ 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, ${\bar{{ \mathcal L }}}_{L,R}$ do not obey the properties (6) and (8) expected of additive, thermal dissipators unless $\delta \approx 0$, as discussed below.

Let us first examine the validity of equation (8). The fixed point ${\hat{r}}_{\alpha }$ satisfying ${\bar{{ \mathcal L }}}_{\alpha }{\hat{r}}_{\alpha }=0$ is of the form

$$\begin{eqnarray}&&{\hat{r}}_{\alpha }=\displaystyle \frac{{{\rm{e}}}^{-{\beta }_{\alpha }({\hat{H}}_{{s}_{\alpha }}-{\mu }_{\alpha }{\hat{n}}_{\alpha })}}{{ \mathcal Z }\,({\beta }_{\alpha },{\hat{H}}_{{s}_{\alpha }},{\mu }_{\alpha },{\hat{n}}_{\alpha })}{\hat{O}}_{\alpha ^{\prime} },\end{eqnarray} \tag{ 73 }$$

where, by analogy with equation (8), we defined an effective local Hamiltonian acting on site α,

$$\begin{eqnarray}&&{\hat{H}}_{{s}_{\alpha }}=\left[{T}_{\alpha }\mathrm{ln}\left(\displaystyle \frac{{\bar{\gamma }}_{\alpha }^{-}}{{\bar{\gamma }}_{\alpha }^{+}}\right)+{\mu }_{\alpha }\right]{\hat{n}}_{\alpha },\end{eqnarray} \tag{ 74 }$$

while ${\hat{O}}_{\alpha ^{\prime} }$ is an arbitrary density operator with support on the other site $\alpha ^{\prime} \ne \alpha$. Using equation (27), we find

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{{\rm{\Gamma }}\to 0}\displaystyle \frac{{\bar{\gamma }}_{\alpha }^{-}}{{\bar{\gamma }}_{\alpha }^{+}}=\left[\displaystyle \frac{{{\rm{e}}}^{{\beta }_{\alpha }(h-{\mu }_{\alpha })}+{\rm{\cosh }}({\beta }_{\alpha }{\rm{\Delta }}/2)}{1+{{\rm{e}}}^{{\beta }_{\alpha }(h-{\mu }_{\alpha })}{\rm{\cosh }}({\beta }_{\alpha }{\rm{\Delta }}/2)}\right]{{\rm{e}}}^{{\beta }_{\alpha }(h-{\mu }_{\alpha })}.\end{eqnarray} \tag{ 75 }$$

This only has the detailed-balance form required for thermalisation if ${\beta }_{\alpha }{\rm{\Delta }}\ll 1$ or if $| {\beta }_{\alpha }(h-{\mu }_{\alpha })| \gg 1$, which corresponds to the white-noise limit. In that case, we have that ${\mathrm{lim}}_{{\rm{\Gamma }}\to 0}\mathrm{ln}({\bar{\gamma }}_{\alpha }^{-}/{\bar{\gamma }}_{\alpha }^{+})\approx {\beta }_{\alpha }(h-{\mu }_{\alpha })$, and the weak-coupling fixed point (73) can be approximated by equation (8) with ${\hat{H}}_{{s}_{\alpha }}=h{\hat{n}}_{\alpha }$ and ${\hat{N}}_{\alpha }={\hat{n}}_{\alpha }$. However, this only corresponds to the 'true' site Hamiltonian ${h}_{\alpha }{\hat{n}}_{\alpha }$ (see equation (10)) if $\delta =0$.

Second, let us discuss the definition of the currents via equation (6). For the particle current, we have, for example,

$$\begin{eqnarray}&&{\langle {\bar{{ \mathcal L }}}_{L}^{\dagger }{\hat{N}}_{S}\rangle }_{\infty }={\langle {{ \mathcal L }}^{\dagger }{\hat{n}}_{L}-{\bar{{ \mathcal L }}}_{{LR}}^{\dagger }{\hat{n}}_{L}\rangle }_{\infty }={J}_{L}^{P}-{\langle {\bar{{ \mathcal L }}}_{{LR}}^{\dagger }{\hat{n}}_{L}\rangle }_{\infty },\end{eqnarray} \tag{ 76 }$$

where we identified $\langle {{ \mathcal L }}^{\dagger }{\hat{n}}_{L}{\rangle }_{\infty }={J}_{L}^{P}$ using the exact master equation (50). Hence, equation (76) differs from the true particle current by an amount

$$\begin{eqnarray}&&-{\langle {\bar{{ \mathcal L }}}_{{LR}}^{\dagger }{\hat{n}}_{L}\rangle }_{\infty }=\displaystyle \frac{g}{2{\rm{\Delta }}}({{\rm{\Gamma }}}_{2}-{{\rm{\Gamma }}}_{1})\,\mathrm{Re}\,{\langle {\hat{c}}_{L}^{\dagger }{\hat{c}}_{R}\rangle }_{\infty }.\end{eqnarray} \tag{ 77 }$$

Clearly, this correction is negligible if ${\bar{{ \mathcal L }}}_{{LR}}\approx 0$ or if ${\rm{\Delta }}/{\rm{\Omega }}\to 0$ so that ${{\rm{\Gamma }}}_{1}\approx {{\rm{\Gamma }}}_{2}$. In addition, equation (77) vanishes as $\delta \to 0$, since one can easily show that $\,\mathrm{Re}\,\langle {\hat{c}}_{L}^{\dagger }{\hat{c}}_{R}{\rangle }_{\infty }\propto \delta /{\rm{\Delta }}$.

On the other hand, we find for the energy current

$$\begin{eqnarray}&&{\langle {\bar{{ \mathcal L }}}_{L}^{\dagger }{\hat{H}}_{S}\rangle }_{\infty }={h}_{L}{\langle {\bar{{ \mathcal L }}}_{L}^{\dagger }{\hat{N}}_{S}\rangle }_{\infty }-\displaystyle \frac{g}{2}\displaystyle \sum _{j=1}^{2}{R}_{1j}^{2}{{\rm{\Gamma }}}_{j}\,\mathrm{Re}\,{\langle {\hat{c}}_{L}^{\dagger }{\hat{c}}_{R}\rangle }_{\infty }.\end{eqnarray} \tag{ 78 }$$

This expression is only correct when ${\bar{{ \mathcal L }}}_{{LR}}\approx 0$ and $\delta \ll h,g$, such that $\langle {\bar{{ \mathcal L }}}_{L}^{\dagger }{\hat{H}}_{S}{\rangle }_{\infty }\approx {{hJ}}_{L}^{P}$. This approximately agrees with the exact Landauer formulae (43) and (44) in the white-noise limit, where ${\rm{\Omega }}\to \infty$, ${ \mathcal T }(h+\omega )={ \mathcal T }(h-\omega )$ and ${f}_{\alpha }(\omega )\approx \mathrm{const}.$

5.3. Non-Markovianity witnesses non-additivity

Since neither ${{ \mathcal L }}_{\mathrm{int}}$ nor ${\bar{{ \mathcal L }}}_{{LR}}$ are in Lindblad form, it may not be possible to cast the overall dissipator ${ \mathcal L }$ into Lindblad form for certain values of the parameters. Therefore, its canonical representation (52) may exhibit one or more negative rates ${\lambda }_{j}^{\pm }(t)$ even as $t\to \infty$. Negativity of the rates ${\lambda }_{j}^{s}(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 [76]. In this framework, the degree of instantaneous non-Markovianity is quantified by [57]

$$\begin{eqnarray}&&\nu (t)=\displaystyle \sum _{s=\pm }\displaystyle \sum _{j=1,2}\max [0,-{\lambda }_{j}^{s}(t)].\end{eqnarray} \tag{ 79 }$$

Asymptotic non-Markovianity (ANM) corresponds to ${\mathrm{lim}}_{t\to \infty }\nu (t)={\nu }_{\infty }\gt 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 ${{\rm{\Omega }}}^{-1}=0$. ANM has also been recently identified in several other open-system scenarios [57–59].

We plot ${\nu }_{\infty }$ in figure 5, as calculated using the EPA rate matrices defined by equation (55). At low temperatures, we find ${\nu }_{\infty }\gt 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 effect is significantly 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 difference, without any engineering or microscopic control of the environment Hamiltonian. Indeed, the ANM effect survives even in the wide-band limit ${\rm{\Omega }}\to \infty$ where each environment is completely unstructured. In figure 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 different temperatures.

Zoom In Zoom Out Reset image size

Since the generators ${{ \mathcal L }}_{L,R}$ (${\bar{{ \mathcal L }}}_{L,R}$) are in Lindblad form, the appearance of asymptotic non-Markovianity is associated with the cross-term ${{ \mathcal L }}_{\mathrm{int}}$ (${\bar{{ \mathcal L }}}_{{LR}}$). Hence, ANM witnesses the violation of additivity in either the global or local pictures of dissipation. Indeed, in figure 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.

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 ${{ \mathcal L }}_{\alpha }$ associated to each bath, which drives the open network towards thermal equilibrium according to equation (7) (or, more generally, equation (60)). Furthermore, ${{ \mathcal L }}_{\alpha }$ determines the currents via equation (6) in the quantum-optical limit. This generator can therefore be interpreted as an individual contribution describing the effect of bath B_α acting in isolation.

Nevertheless, additivity is violated away from equilibrium due to an additional interference term ${{ \mathcal L }}_{\mathrm{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 ${{ \mathcal L }}_{\mathrm{int}}$ leaves the energy-eigenbasis populations invariant and therefore the global, additive Lindblad model obtained by setting ${{ \mathcal L }}_{\mathrm{int}}=0$ still leads to accurate approximations for the boundary currents in the quantum-optical regime. Such an additive model coincides in the limit ${\rm{\Gamma }}\to 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 identification of ${\bar{{ \mathcal L }}}_{\alpha }$ as an 'individual' bath contribution is not fully justified, since it does not accurately represent the effect of a single bath acting in isolation. Nevertheless, an additive, local Lindblad model is recovered when the modes are near-resonant, i.e. $\delta \ll h,g$, and when the noise correlation time is much smaller than the inverse frequency splitting ${{\rm{\Delta }}}^{-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 second-order 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 system-reservoir density matrix factorises, it is more accurately described as a projection ${ \mathcal P }$ onto product states of the form

$$\begin{eqnarray}&&{ \mathcal P }\hat{\rho }(t)={\hat{\rho }}_{S}(t){\hat{\rho }}_{B},\end{eqnarray} \tag{ 80 }$$

where ${\hat{\rho }}_{B}$ is the initial state of the environment, in combination with a low-order perturbation expansion in the system-reservoir interaction [19].

The breakdown of additivity for arbitrarily weak system-bath coupling Γ indicates that equation (80) 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 flowing between them. As time increases, these correlations cause the joint state $\hat{\rho }(t)$ to move arbitrarily far away from the subspace spanned by states of the form (80) with ${\hat{\rho }}_{B}={\hat{\rho }}_{L}{\hat{\rho }}_{R}$. A very interesting open question is whether a projection of the type (80) can still facilitate a good perturbative description of the long-time dynamics, but with a correlated reservoir state ${\hat{\rho }}_{B}$.

Although the evidence presented here pertains only to the specific 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 [77, 78]). We note also a very recent study of a different bosonic model of coupled mechanical oscillators that displays coherences in the NESS [79], 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 affect the thermodynamic power and efficiency 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 fluctuations of the energy and particle currents. Finally, the framework presented herein appears to be an ideal setting to explore strong-coupling effects in thermodynamics [80–85], since system-environment correlations in both equilibrium and far-from-equilibrium states may be taken into account.

We gratefully acknowledge the stimulating comments of and conversations with: Jonatan Brask, Nicolas Brunner, Luis Correa, John Goold, Géraldine Haack, Patrick Hofer, Susana Huelga, 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.

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 is

$$\begin{eqnarray}&&\hat{H}={\hat{H}}_{S}+\displaystyle \sum _{\alpha }({\hat{H}}_{{B}_{\alpha }}+{\hat{H}}_{S{B}_{\alpha }}),\end{eqnarray} \tag{ A1 }$$

where ${\hat{H}}_{S}$ and ${\hat{H}}_{{B}_{\alpha }}$ are respectively the Hamiltonians of system and baths, while ${\hat{H}}_{S{B}_{\alpha }}$ 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 figure A1 for an illustration.

Zoom In Zoom Out Reset image size

Suppose that there exists a conserved quantity

$$\begin{eqnarray}&&\hat{X}={\hat{X}}_{S}+\displaystyle \sum _{\alpha }{\hat{X}}_{{B}_{\alpha }},\end{eqnarray} \tag{ A2 }$$

such that $[\hat{H},\hat{X}]=0$, which is also locally conserved in the sense that $[{\hat{H}}_{S},{\hat{X}}_{S}]=[{\hat{H}}_{{B}_{\alpha }},{\hat{X}}_{{B}_{\alpha }}]=0$. It follows that

$$\begin{eqnarray}&&{\partial }_{t}{\hat{X}}_{{B}_{\alpha }}={\rm{i}}[{\hat{H}}_{S{B}_{\alpha }},{\hat{X}}_{{B}_{\alpha }}]=-{\rm{i}}[{\hat{H}}_{S{B}_{\alpha }},{\hat{X}}_{S}].\end{eqnarray} \tag{ A3 }$$

For simplicity, we assume that ${\hat{X}}_{S}$ is a one-body observable of the form

$$\begin{eqnarray}&&{\hat{X}}_{S}=\displaystyle \sum _{k}{\hat{x}}_{k},\end{eqnarray} \tag{ A4 }$$

where the local 'charge' ${\hat{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 for ${\hat{x}}_{k}$ takes the form of a continuity equation:

$$\begin{eqnarray}&&{\partial }_{t}{\hat{x}}_{k}={\hat{J}}_{{B}_{\alpha }\to k}^{X}+{\hat{J}}_{S\to k}^{X}.\end{eqnarray} \tag{ A5 }$$

The current operators defined by ${\hat{J}}_{{B}_{\alpha }\to k}^{X}={\rm{i}}[{\hat{H}}_{S{B}_{\alpha }},{\hat{x}}_{k}]$ and ${\hat{J}}_{S\to k}^{X}={\rm{i}}[{\hat{H}}_{S},{\hat{x}}_{k}]$ describe the flow of charge from B_α to site k and from the other sites of $S$ to site k, respectively. In particular, the operators ${\hat{J}}_{S\to k}^{X}$ correspond to the internal current observables denoted schematically by ${\hat{J}}_{S}$ in section 2. By definition, these observables have zero mean unless the quantum state has some coherence in the eigenbasis of ${\hat{H}}_{S}$. This follows because, by the cyclic invariance of the trace,

$$\begin{eqnarray}&&{\langle {\hat{J}}_{S\to k}^{X}\rangle }_{t}=\mathrm{Tr}\{{\rm{i}}[{\hat{\rho }}_{S}(t),{\hat{H}}_{S}]{\hat{x}}_{k}\},\end{eqnarray} \tag{ A6 }$$

which clearly vanishes if $[{\hat{H}}_{S},{\hat{\rho }}_{S}(t)]=0$. However, these expectation values do not vanish whenever the NESS supports a current due to the influx of charge from the external reservoirs. Indeed, summing the expectation value of equation (A5) over all sites in V_α, we obtain

$$\begin{eqnarray}&&{\langle {\partial }_{t}{\hat{X}}_{{B}_{\alpha }}\rangle }_{\infty }=\displaystyle \sum _{k\in {V}_{\alpha }}{\langle {\hat{J}}_{S\to k}^{X}\rangle }_{\infty }.\end{eqnarray} \tag{ A7 }$$

Here, we have used equation (A3) and the fact that expectation values of system observables are time-independent in the NESS. Identifying ${J}_{\alpha }^{X}=-\langle {\partial }_{t}{\hat{X}}_{{B}_{\alpha }}{\rangle }_{\infty }$ as the current entering the system from B_α, the assertion (3) is confirmed. We conclude that boundary-driven currents imply energy-eigenbasis coherence in the NESS.

For clarity, it is useful to consider the specific case of a lattice with two-body interactions described by the generic Hamiltonian

$$\begin{eqnarray}&&{\hat{H}}_{S}=\tfrac{1}{2}\displaystyle \sum _{k,l}{\hat{h}}_{{kl}},\end{eqnarray} \tag{ A8 }$$

where ${\hat{h}}_{{kl}}={\hat{h}}_{{lk}}={\hat{h}}_{{kl}}^{\dagger }$ is an operator with support on sites k and l. In this case, we have

$$\begin{eqnarray}&&{\hat{J}}_{S\to k}^{X}=\displaystyle \sum _{l\in {A}_{k}}{\hat{J}}_{l\to k}^{X},\end{eqnarray} \tag{ A9 }$$

where A_k denotes the the neighbourhood of site k, i.e. all sites l such that ${\hat{h}}_{{kl}}\ne 0$, while ${\hat{J}}_{l\to k}^{X}={\rm{i}}[{\hat{h}}_{{lk}},{\hat{x}}_{k}]$ denotes the current flowing from site l to site k. According to equation (A6), eigenstates $| E\rangle$ of ${\hat{H}}_{S}$ obey the constraint

$$\begin{eqnarray}&&\displaystyle \sum _{l\in {A}_{k}}\langle E| {\hat{J}}_{l\to k}^{X}| E\rangle =0.\end{eqnarray} \tag{ A10 }$$

This constraint allows for non-zero net currents around loops in the network (see, for example, [86]), such that the local charge distribution is stationary (red arrows in figure A1). 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 figure A1).

The simplest example of this principle is that of a 1D chain with open (i.e. non-periodic) boundary conditions. Since the boundary sites of such a chain have only a single neighbour, equation (A10) reduces to the identity

$$\begin{eqnarray}&&\langle E| {\hat{J}}_{k\to k+1}^{X}| E\rangle =0.\end{eqnarray} \tag{ A11 }$$

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 ${\hat{X}}_{S}=\tfrac{1}{2}{\sum }_{k,l}{\hat{x}}_{{kl}}$. For a network with two-body interactions, this class includes the Hamiltonian itself, which is of course always conserved.

Here, we provide details of the exact solution presented in section 3. The equations of motion read as

$$\begin{eqnarray}&&{\rm{i}}{\partial }_{t}{\hat{d}}_{j}(t)={E}_{j}{\hat{d}}_{j}(t)+\displaystyle \sum _{q}{\tau }_{q}{\hat{b}}_{q,j}(t),\end{eqnarray} \tag{ B1 }$$

$$\begin{eqnarray}&&{\rm{i}}{\partial }_{t}{\hat{b}}_{q,j}(t)={\omega }_{q}{\hat{b}}_{q,j}(t)+{\tau }_{q}^{* }{\hat{d}}_{j}(t).\end{eqnarray} \tag{ B2 }$$

These are readily solved by transforming to Laplace space [46], e.g. ${\tilde{d}}_{j}(z)={\int }_{0}^{\infty }{\rm{d}}t\ {{\rm{e}}}^{-{zt}}{\hat{d}}_{j}(t)$. After rearranging the resulting linear, algebraic system of equations, we transform back to the time domain to obtain

$$\begin{eqnarray}&&{\hat{d}}_{j}(t)={G}_{j}(t){\hat{d}}_{j}+\displaystyle \sum _{q}{\tau }_{q}{K}_{j}({\omega }_{q},t){\hat{b}}_{q,j},\end{eqnarray} \tag{ B3 }$$

$$\begin{eqnarray}&&{\hat{b}}_{q,j}(t)={{\rm{e}}}^{-{\rm{i}}{\omega }_{q}t}{\hat{b}}_{q,j}+{\tau }_{q}^{* }{K}_{j}({\omega }_{q},t){\hat{d}}_{j}+{\tau }_{q}^{* }\displaystyle \sum _{p}{\tau }_{p}{I}_{j}({\omega }_{q},{\omega }_{p},t){\hat{b}}_{p,j}.\end{eqnarray} \tag{ B4 }$$

Here, G_j(t) denotes the inverse Laplace transform of ${\tilde{G}}_{j}(z)$, where

$$\begin{eqnarray}&&{\tilde{G}}_{j}(z)={[z+{\rm{i}}{E}_{j}+\tilde{W}(z)]}^{-1},\end{eqnarray} \tag{ B5 }$$

$$\begin{eqnarray}&&\tilde{W}(z)=\int {\rm{d}}\omega \ \displaystyle \frac{{ \mathcal J }(\omega )}{z+{\rm{i}}\omega }.\end{eqnarray} \tag{ B6 }$$

We also introduced the functions

$$\begin{eqnarray}&&{K}_{j}(\omega ,t)=-{\rm{i}}{\int }_{0}^{t}{\rm{d}}t^{\prime} \,{{\rm{e}}}^{-{\rm{i}}\omega (t-t^{\prime} )}{G}_{j}(t^{\prime} ),\end{eqnarray} \tag{ B7 }$$

$$\begin{eqnarray}&&{I}_{j}(\omega ,\omega ^{\prime} ,t)=-{\rm{i}}{\int }_{0}^{t}{\rm{d}}t^{\prime} \,{{\rm{e}}}^{-{\rm{i}}\omega (t-t^{\prime} )}{K}_{j}(\omega ^{\prime} ,t^{\prime} ).\end{eqnarray} \tag{ B8 }$$

By inspection of equation (B3) one sees that ${\mathsf{G}}(t)$ is indeed given by equation (23). Expressions (29) and (31) for the correlation matrix ${\mathsf{C}}(t)$ and double occupancy D(t) follow directly upon plugging equation (B3) into the definitions (21) and (22). Note that the above equations hold for an arbitrary spectral density ${ \mathcal J }(\omega )$, which determines the propagator via equations (B5) and (B6).

Let us now compute the propagator for the specific problem at hand, with spectral density (19). The integral (B6) evaluates to

$$\begin{eqnarray}&&\tilde{W}(z)=\displaystyle \frac{{\rm{\Gamma }}}{2{\rm{\Omega }}}(\sqrt{{z}^{2}+{{\rm{\Omega }}}^{2}}-z),\end{eqnarray} \tag{ B9 }$$

where the positive (negative) sign for the square root is used for $\,\mathrm{Re}\,z\gt 0$ ($\,\mathrm{Re}\,z\lt 0$). This implies that any poles of ${\tilde{G}}_{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 that ${\tilde{G}}_{j}(z)$ has no poles. As a result, ${\tilde{G}}_{j}(z)$ is holomorphic everywhere in the complex z plane except for along a cut connecting the branch points at $z=\pm {\rm{i}}{\rm{\Omega }}$. Therefore, ${G}_{j}(t)\to 0$ as $t\to \infty$ 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\gt 0$,

$$\begin{eqnarray}&&{G}_{j}(t)=\displaystyle \frac{1}{2\pi {\rm{i}}}{\int }_{{{ \mathcal C }}_{1}}{\rm{d}}z\,{{\rm{e}}}^{{zt}}{\tilde{G}}_{j}(z)=-\displaystyle \frac{1}{2\pi {\rm{i}}}{\int }_{{{ \mathcal C }}_{2}}{\rm{d}}z\,{{\rm{e}}}^{{zt}}{\tilde{G}}_{j}(z).\end{eqnarray} \tag{ B10 }$$

Here, ${{ \mathcal C }}_{1}$ denotes a Bromwich contour, while ${{ \mathcal 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 figure B1). After a change of variables to $\omega ={\rm{i}}z$, the integral along ${{ \mathcal C }}_{2}$ takes the form (24), with ${\varphi }_{j}(\omega )$ given by equation (35).

Zoom In Zoom Out Reset image size

It remains to justify equations (32) and (37). By taking the Laplace transform of equation (24), we obtain

$$\begin{eqnarray}&&{\tilde{G}}_{j}(z)=\int {\rm{d}}\omega \,\displaystyle \frac{{\varphi }_{j}(\omega )}{z+{\rm{i}}\omega }.\end{eqnarray} \tag{ B11 }$$

We therefore deduce that, in general,

$$\begin{eqnarray}\mathop{\mathrm{lim}}\limits_{\epsilon \searrow 0}{\tilde{G}}_{j}(\epsilon -{\rm{i}}\omega ) & = & {\displaystyle \int }_{0}^{\infty }{\rm{d}}t\,{{\rm{e}}}^{-{\rm{i}}\omega t}{G}_{j}(t)\\ & = & \pi {\varphi }_{j}(\omega )+{\rm{i}}{\vartheta }_{j}(\omega ).\end{eqnarray} \tag{ B12 }$$

Using the first line, we obtain equation (32) from (29), while comparison of the second line with equations (B5) and (B9) yields equation (37) for the uniform-chain environment model.

In this appendix, we prove equation (39). That is, when ${\beta }_{\alpha }=\beta$, ${\mu }_{\alpha }=\mu$ and ${f}_{\alpha }(\omega )=f(\omega )$, the open system equilibrates to a state ${\hat{\rho }}_{S}^{\infty }$ given by the reduction of a global Gibbs distribution. The following arguments can also be used to deduce equation (60) from (59).

Since the marginal of a Gaussian state is itself Gaussian [87], it suffices to check that equation (39) yields the correct correlation matrix (38). Due to the symmetry $[\hat{H},{\hat{N}}_{j}]=0$ of the global Hamiltonian, it follows that its equilibrium correlation matrix is diagonal, i.e. $\mathrm{Tr}[{{\rm{e}}}^{-\beta (\hat{H}-\mu \hat{N})}{\hat{d}}_{j}^{\dagger }{\hat{d}}_{k}]\propto {\delta }_{{jk}}$. We explicitly compute the diagonal elements using a basis transformation to the eigenmodes of $\hat{H}$, given by [18] ${\hat{d}}_{j}={\sum }_{n}\langle 0| {\hat{d}}_{j}| n\rangle {\hat{e}}_{n}$, where $| n\rangle ={\hat{e}}_{n}^{\dagger }| 0\rangle$ is a single-particle eigenstate of $\hat{H}$ (see equation (25)) and the ladder operators ${\hat{e}}_{n}$ diagonalise the total Hamiltonian as $\hat{H}={\sum }_{n}{\varepsilon }_{n}{\hat{e}}_{n}^{\dagger }{\hat{e}}_{n}$. Substituting the above and using equation (25), straightforward manipulations lead to

$$\begin{eqnarray}&&\displaystyle \frac{\mathrm{Tr}[{{\rm{e}}}^{-\beta (\hat{H}-\mu \hat{N})}{\hat{d}}_{j}^{\dagger }{\hat{d}}_{j}]}{{ \mathcal Z }\,(\beta ,\hat{H},\mu ,\hat{N})}=\int {\rm{d}}\omega \,{\varphi }_{j}(\omega )f(\omega ),\end{eqnarray} \tag{ C1 }$$

in agreement with equation (39), which completes the proof.

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 ${\langle {\hat{J}}_{S}^{P}\rangle }_{t}=g\,\mathrm{Im}\,[{C}_{12}(t)]$. In the limit $t\to \infty$, the particle current can thus be found directly from equation (32):

$$\begin{eqnarray}&&{J}_{S}^{P}=\displaystyle \frac{\pi {g}^{2}}{2{\rm{\Delta }}}\int {\rm{d}}\omega \,{ \mathcal J }(\omega )[{\varphi }_{1}(\omega ){\vartheta }_{2}(\omega )-{\varphi }_{2}(\omega ){\vartheta }_{1}(\omega )][{f}_{L}(\omega )-{f}_{R}(\omega )].\end{eqnarray} \tag{ D1 }$$

To compute the expectation value of equations (40) and (41), we make use of equations (B3) and (B4) and obtain, for example

$$\begin{eqnarray}&&{\langle {\hat{J}}_{L}^{P}\rangle }_{t}=2\,\mathrm{Im}\,{\left[{{\mathsf{R}}}^{{\mathsf{T}}}\int {\rm{d}}\omega { \mathcal J }(\omega ){{\mathsf{K}}}^{\dagger }(\omega ,t){\mathsf{F}}(\omega )\left({{\rm{e}}}^{-{\rm{i}}\omega t}{\mathbb{1}}+\int {\rm{d}}\omega ^{\prime} { \mathcal J }(\omega ^{\prime} ){\mathsf{I}}(\omega ^{\prime} ,\omega ,t)\right){\mathsf{R}}\right]}_{11},\end{eqnarray} \tag{ D2 }$$

$$\begin{eqnarray}&&{\langle {\hat{J}}_{L}^{E}\rangle }_{t}=2\,\mathrm{Im}\,{\left[{{\mathsf{R}}}^{{\mathsf{T}}}\int {\rm{d}}\omega { \mathcal J }(\omega ){{\mathsf{K}}}^{\dagger }(\omega ,t){\mathsf{F}}(\omega )\left(\omega {{\rm{e}}}^{-{\rm{i}}\omega t}{\mathbb{1}}+\int {\rm{d}}\omega ^{\prime} \omega ^{\prime} { \mathcal J }(\omega ^{\prime} ){\mathsf{I}}(\omega ^{\prime} ,\omega ,t)\right){\mathsf{R}}\right]}_{11},\end{eqnarray} \tag{ D3 }$$

where ${\mathsf{K}}(\omega ,t)=\mathrm{diag}[{K}_{1}(\omega ,t),{K}_{2}(\omega ,t)]$ and ${\mathsf{I}}(\omega ^{\prime} ,\omega ,t)=\mathrm{diag}[{I}_{1}(\omega ^{\prime} ,\omega ,t),{I}_{2}(\omega ^{\prime} ,\omega ,t)]$. In the limit $t\to \infty$, we find that

$$\begin{eqnarray}&&{J}_{L}^{P}=2\pi \displaystyle \sum _{k,l}{R}_{k1}{R}_{l1}\int {\rm{d}}\omega \,{ \mathcal J }(\omega ){F}_{{kl}}(\omega )\{{\varphi }_{k}(\omega )-{ \mathcal J }(\omega )[{\vartheta }_{k}(\omega ){\vartheta }_{l}(\omega )+{\pi }^{2}{\varphi }_{k}(\omega ){\varphi }_{l}(\omega )]\},\end{eqnarray} \tag{ D4 }$$

$$\begin{eqnarray}&&{J}_{L}^{E}=2\pi \displaystyle \sum _{k,l}{R}_{k1}{R}_{l1}\int {\rm{d}}\omega \,\omega { \mathcal J }(\omega ){F}_{{kl}}(\omega )\{{\varphi }_{k}(\omega )-{ \mathcal J }(\omega )[{\vartheta }_{k}(\omega ){\vartheta }_{l}(\omega )+{\pi }^{2}{\varphi }_{k}(\omega ){\varphi }_{l}(\omega )]\}.\end{eqnarray} \tag{ D5 }$$

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 equation (37).

In this appendix, we justify the EPA (53) and prove that the error is of order $O({\rm{\Gamma }}/{\rm{\Omega }})$. The integral (24), with ${\varphi }_{j}(\omega )$ given by equation (35), can be estimated with high accuracy for ${\rm{\Gamma }}\ll {\rm{\Omega }}$ by analytically continuing the integrand into the complex ω plane. Choosing (arbitrarily) the branch of the square root function whose real part is positive for $\,\mathrm{Im}\,\omega \lt 0$, we integrate along a closed contour encircling the pole at $\omega ={E}_{j}^{{\prime} }-{\rm{i}}{{\rm{\Gamma }}}_{j}/2$ in the anti-clockwise sense (see figure E1). This integral differs from the exact G_j(t) by the contribution along the semi-circular arc ${{ \mathcal C }}_{\mathrm{err}}$. We obtain

$$\begin{eqnarray}&&{G}_{j}(t)={p}_{j}{{\rm{e}}}^{-{\rm{i}}{E}_{j}^{{\prime} }t-{{\rm{\Gamma }}}_{j}t}+{\mathrm{err}}_{j}(t).\end{eqnarray} \tag{ E1 }$$

Here, p_j is the residue at the pole,

$$\begin{eqnarray}&&{p}_{j}=\sqrt{\displaystyle \frac{{{\rm{\Omega }}}^{2}-{({E}_{j}^{{\prime} }-{\rm{i}}{{\rm{\Gamma }}}_{j}/2)}^{2}}{{{\rm{\Omega }}}^{2}-{\rm{\Gamma }}{\rm{\Omega }}-{E}_{j}^{2}}},\end{eqnarray} \tag{ E2 }$$

while the error term is given by

$$\begin{eqnarray}&&{\mathrm{err}}_{j}(t)=-{\int }_{{{ \mathcal C }}_{\mathrm{err}}}{\rm{d}}\omega \ {{\rm{e}}}^{-{\rm{i}}\omega t}{\varphi }_{j}(\omega ).\end{eqnarray} \tag{ E3 }$$

A trivial rearrangement of terms leads to

$$\begin{eqnarray}&&{G}_{j}(t)={{\rm{e}}}^{-{\rm{i}}{E}_{j}^{{\prime} }t-{{\rm{\Gamma }}}_{j}t}+({p}_{j}-1){{\rm{e}}}^{-{\rm{i}}{E}_{j}^{{\prime} }t-{{\rm{\Gamma }}}_{j}t}+{\mathrm{err}}_{j}(t).\end{eqnarray} \tag{ E4 }$$

The first term above is the EPA contribution (53). The second term is a small relative correction to the leading-order result, which is clearly of first order in ${\rm{\Gamma }}/{\rm{\Omega }}$ and decays exponentially in time. Finally, the remaining absolute error $| {\mathrm{err}}_{j}(t)|$ can be bounded as

$$\begin{eqnarray}| {\mathrm{err}}_{j}(t)| & \leqslant & \displaystyle \frac{{\rm{\Omega }}}{1-{\rm{\Gamma }}/{\rm{\Omega }}}{\displaystyle \int }_{0}^{\pi }{\rm{d}}\theta \ \left|\displaystyle \frac{{{\rm{e}}}^{-{\rm{i}}(\theta +{\rm{\Omega }}t{{\rm{e}}}^{-{\rm{i}}\theta })}{{ \mathcal J }}_{{\rm{N}}}({\rm{\Omega }}{{\rm{e}}}^{-{\rm{i}}\theta })}{{({\rm{\Omega }}{{\rm{e}}}^{-{\rm{i}}\theta }-{E}_{j}^{{\prime} })}^{2}+{{\rm{\Gamma }}}_{j}^{2}/4}\right|\\ & = & \displaystyle \frac{{\rm{\Gamma }}{{\rm{\Omega }}}^{2}}{\sqrt{2}\pi ({\rm{\Omega }}-{\rm{\Gamma }})}{\displaystyle \int }_{0}^{\pi }{\rm{d}}\theta \,\displaystyle \frac{{{\rm{e}}}^{-{\rm{\Omega }}t{\rm{\sin }}\theta }{{\rm{\sin }}}^{1/2}\theta }{| {({\rm{\Omega }}{{\rm{e}}}^{-{\rm{i}}\theta }-{E}_{j}^{{\prime} })}^{2}+{{\rm{\Gamma }}}_{j}^{2}/4| }\\ & & \lt \displaystyle \frac{{\rm{\Gamma }}}{{\rm{\Omega }}}\times \displaystyle \frac{1}{(1-{\rm{\Gamma }}/{\rm{\Omega }})}\displaystyle \frac{{{\rm{\Omega }}}^{2}}{{({\rm{\Omega }}-{E}_{j}^{{\prime} })}^{2}+{{\rm{\Gamma }}}_{j}^{2}/4}u({\rm{\Omega }}t).\end{eqnarray} \tag{ E5 }$$

On the first line above, we changed variables to $\omega ={\rm{\Omega }}{{\rm{e}}}^{-{\rm{i}}\theta }$ and used the fact that $| \int {\rm{d}}x\,f(x)| \leqslant \int {\rm{d}}x\,| f(x)|$, while the third line follows from setting the denominator of the integrand to its minimum value and defining the dimensionless function

$$\begin{eqnarray}u(\tau ) & = & \displaystyle \frac{1}{\sqrt{2}\pi }{\displaystyle \int }_{0}^{\pi }{\rm{d}}\theta \,{{\rm{e}}}^{-\tau \sin \theta }{\sin }^{1/2}\theta \\ & = & \displaystyle \frac{\sqrt{\pi \tau }}{2}\left({I}_{-\tfrac{3}{4}}(\tau /2){I}_{\tfrac{1}{4}}(\tau /2)-{I}_{-\tfrac{1}{4}}(\tau /2){I}_{\tfrac{3}{4}}(\tau /2)\right),\end{eqnarray} \tag{ E6 }$$

where I_s(x) is a modified Bessel function of the first kind. The function $u(\tau )$ is positive, monotonically decreasing and vanishes as $\tau \to \infty$ with the limiting behaviour ${\mathrm{lim}}_{\tau \to \infty }{\tau }^{3/2}u{(\tau )=(2\pi )}^{-1/2}$. Hence, by replacing $u({\rm{\Omega }}t)$ in equation (E5) by $u(0)\approx 0.54$ one obtains a time-independent upper bound on $| {\mathrm{err}}_{j}(t)|$ that is of order $O({\rm{\Gamma }}/{\rm{\Omega }})$, as claimed. Note, however, that the magnitude of the error also depends on ${E}_{j}/{\rm{\Omega }}$. In particular, the error is larger when E_j is closer to the band edge at $\pm {\rm{\Omega }}$.

Zoom In Zoom Out Reset image size

This appendix demonstrates the entropy-production inequality (64), starting from the Spohn inequality (63). We first demonstrate that ${{ \mathcal L }}_{\mathrm{int}}$ does not generate a positive evolution, i.e. the map ${{\rm{e}}}^{{{ \mathcal L }}_{\mathrm{int}}t}$ is not positive. We write ${{ \mathcal L }}_{\mathrm{int}}$ in diagonal form

$$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{int}}=\displaystyle \sum _{j=1}^{4}{\lambda }_{j}{ \mathcal D }[{\hat{L}}_{j}],\end{eqnarray} \tag{ F1 }$$

where ${\lambda }_{1}=-{\lambda }_{2}=-{\lambda }_{3}={\lambda }_{4}=| {{\rm{\Lambda }}}_{12}^{+}|$, ${\hat{L}}_{1}={\hat{L}}_{2}^{\dagger }=\tfrac{1}{\sqrt{2}}({\hat{d}}_{1}-{{\rm{e}}}^{-{\rm{i}}\theta }{\hat{d}}_{2})$, ${\hat{L}}_{3}={\hat{L}}_{4}^{\dagger }=\tfrac{1}{\sqrt{2}}({\hat{d}}_{1}+{{\rm{e}}}^{-{\rm{i}}\theta }{\hat{d}}_{2})$, and $\theta =\arg ({{\rm{\Lambda }}}_{12}^{+})$. Note that these Lindblad operators satisfy the canonical anti-commutation relations $\{{\hat{L}}_{j},{\hat{L}}_{k}^{\dagger }\}={\delta }_{{jk}}$, and are thus linearly independent, i.e. $\mathrm{Tr}[{\hat{L}}_{j}^{\dagger }{\hat{L}}_{k}]=2{\delta }_{{jk}}$. Therefore, ${{ \mathcal L }}_{\mathrm{int}}$ generates a positive evolution if and only if

$$\begin{eqnarray}&&\displaystyle \sum _{j=1}^{4}{\lambda }_{j}{| \langle a| {\hat{L}}_{j}| b\rangle | }^{2}\geqslant 0,\end{eqnarray} \tag{ F2 }$$

for every orthonormal pair of states $| a\rangle$ and $| b\rangle$ [88]. However, it is straightforward to find a counterexample. Consider, for instance, the choice $| b\rangle ={\hat{L}}_{3}^{\dagger }| 0\rangle$ and $| a\rangle =| 0\rangle$. Then $\langle a| {\hat{L}}_{j}| b\rangle ={\delta }_{j3}$ and, since ${\lambda }_{3}\lt 0$, inequality (F2) does not hold. Hence, the Spohn inequality cannot be directly applied to ${{ \mathcal L }}_{\mathrm{int}}$ [75].

Nevertheless, the Lindblad generators ${{ \mathcal L }}_{\alpha }$ do satisfy the Spohn inequality (63). Hence, taking ${\hat{\rho }}_{S}(t)={\hat{\rho }}_{S}^{\infty }$ and summing over the baths, we obtain

$$\begin{eqnarray}&&\mathrm{Tr}[\mathrm{ln}{\hat{\rho }}_{S}^{\infty }{{ \mathcal L }}_{\mathrm{int}}{\hat{\rho }}_{S}^{\infty }]-\displaystyle \sum _{\alpha =L,R}{\beta }_{\alpha }{\langle {{ \mathcal L }}_{\alpha }^{\dagger }({\hat{H}}_{S,\alpha }^{* }-{\mu }_{\alpha }{\hat{N}}_{S})\rangle }_{\infty }\geqslant 0.\end{eqnarray} \tag{ F3 }$$

To obtain the first term on the LHS, we have used the stationary property ${\sum }_{\alpha }{{ \mathcal L }}_{\alpha }{\hat{\rho }}_{S}^{\infty }=-{{ \mathcal L }}_{\mathrm{int}}{\hat{\rho }}_{S}^{\infty }$. In the second term, we introduced the quantum Hamiltonian of mean force [89]

$$\begin{eqnarray}&&{\hat{H}}_{S,\alpha }^{* }=-{T}_{\alpha }\mathrm{ln}\left\{\displaystyle \frac{{\mathrm{Tr}}_{B}[{{\rm{e}}}^{-{\beta }_{\alpha }(\hat{H}-{\mu }_{\alpha }\hat{N})}]}{{ \mathcal Z }\,({\beta }_{\alpha },{\hat{H}}_{B},{\mu }_{\alpha },{\hat{N}}_{B})}\right\}+{\mu }_{\alpha }{\hat{N}}_{S},\end{eqnarray} \tag{ F4 }$$

which is defined so that (see equation (60))

$$\begin{eqnarray}&&{\hat{r}}_{\alpha }=\displaystyle \frac{{{\rm{e}}}^{-{\beta }_{\alpha }({\hat{H}}_{S,\alpha }^{* }-{\mu }_{\alpha }{\hat{N}}_{S})}}{{ \mathcal Z }\,({\beta }_{\alpha },{\hat{H}}_{S,\alpha }^{* },{\mu }_{\alpha },{\hat{N}}_{S})}.\end{eqnarray} \tag{ F5 }$$

Let us examine the first 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 ${{ \mathcal L }}_{\mathrm{int}}$. This term is non-negative, which can be seen by the following geometrical argument, illustrated in figure F1. According to equation (68), ${{ \mathcal L }}_{\mathrm{int}}$ changes only the coherences. Hence, it suffices to consider the action of ${{ \mathcal L }}_{\mathrm{int}}$ in the two-dimensional subspace spanned by the single-particle states $| {E}_{j}\rangle ={\hat{d}}_{j}^{\dagger }| 0\rangle$. The restriction of the density matrix ${\hat{\rho }}_{S}^{\infty }$ to the single-particle subspace is represented by a Bloch vector, ${\bf{v}}$, which we choose to lie perpendicular to the y axis. Specifically, the projection of ${\bf{v}}$ along the z axis is given by $\tfrac{1}{2}+{C}_{11}-{C}_{22}$, while the projection along the x axis is C₁₂, where all quantities are evaluated in the limit $t\to \infty$.

Zoom In Zoom Out Reset image size

Consider now the flow generated by ${{ \mathcal L }}_{\mathrm{int}}$ over a small time interval $\delta t$, i.e. the vector ${\bf{u}}$ such that the state ${{\rm{e}}}^{{{ \mathcal L }}_{\mathrm{int}}\delta t}{\hat{\rho }}_{S}^{\infty }$ corresponds to the shifted Bloch vector ${\bf{v}}+{\bf{u}}\delta t+O(\delta {t}^{2})$. According to equation (68), ${\bf{u}}$ lies in the plane perpendicular to the z axis and points in a direction determined by the complex argument of ${{\rm{\Lambda }}}_{12}^{+}$. Using equation (48), we deduce that ${{\rm{\Lambda }}}_{12}^{+}=(\mathrm{Tr}[{\rm{\Gamma }}]/2-{\rm{i}}{\rm{\Delta }}){C}_{12}$. Thus, the angle subtended by ${\bf{u}}$ from the x axis is $\phi =\arctan (2{\rm{\Delta }}/\mathrm{Tr}[{\rm{\Gamma }}])$, such that $0\leqslant \phi \leqslant \pi /2$. We conclude that the flow generated by ${{ \mathcal L }}_{\mathrm{int}}$ never moves the Bloch vector away the surface of the Bloch sphere, i.e. the length of ${\bf{v}}$ does not decrease. Since the von Neumann entropy is a monotonically decreasing function of the length of ${\bf{v}}$, we conclude that it is a non-increasing function under this flow.

Now let us demonstrate that the derivative of the von Neumann entropy generated by ${{ \mathcal L }}_{\mathrm{int}}$ behaves as $o({\rm{\Gamma }})$ in the limit ${\rm{\Gamma }}\to 0$. For this it is convenient to use an explicit representation of the Gaussian NESS,

$$\begin{eqnarray}&&{\hat{\rho }}_{S}^{\infty }=\displaystyle \frac{\exp (-{\hat{{\bf{d}}}}^{\dagger }{{\mathsf{P}}}^{{\mathsf{T}}}\hat{{\bf{d}}})}{\det [1+{{\rm{e}}}^{-{\mathsf{P}}}]},\end{eqnarray} \tag{ F6 }$$

where ${\mathsf{P}}$ is a positive semi-definite matrix that satisfies ${\mathsf{C}}={[{{\rm{e}}}^{{\mathsf{P}}}+1]}^{-1}$. Then, making use of equation (68), we find the explicit representation

$$\begin{eqnarray}&&\mathrm{Tr}[\mathrm{ln}{\hat{\rho }}_{S}^{\infty }{{ \mathcal L }}_{\mathrm{int}}{\hat{\rho }}_{S}^{\infty }]=-2\,\mathrm{Re}\,[{{\rm{\Lambda }}}_{12}^{+}{P}_{21}].\end{eqnarray} \tag{ F7 }$$

Clearly, ${{\rm{\Lambda }}}_{12}^{+}=O({\rm{\Gamma }})$ by its definition (55). Furthermore, since the coherence in the stationary state ${C}_{12}={{\rm{\Lambda }}}_{12}^{+}/[\mathrm{Tr}[{\rm{\Gamma }}]/2-{\rm{i}}{\rm{\Delta }}]$ also vanishes as ${\rm{\Gamma }}\to 0$, we conclude that ${P}_{21}\to 0$ as ${\rm{\Gamma }}\to 0$. Overall, the expression (F7) is therefore of order $o({\rm{\Gamma }})$. Hence, this term vanishes upon dividing the relation (F3) by Γ and taking the limit ${\rm{\Gamma }}\to 0$. Taking into account the fact that ${\mathrm{lim}}_{{\rm{\Gamma }}\to 0}{\hat{H}}_{S,\alpha }^{* }={\hat{H}}_{S}$, the inequality (64) is thus recovered.