Wigner entropy production and heat transport in linear quantum lattices

TL;DR

This paper develops a phase-space method to calculate Wigner entropy production in linear quantum networks, revealing how internal interactions influence non-equilibrium steady states and energy transport regimes.

Contribution

It introduces a simple closed-form expression for entropy production in linear quantum networks and analyzes the impact of internal couplings on non-equilibrium steady states.

Findings

Entropy production scales differently in ballistic and diffusive regimes.

Internal couplings are crucial for maintaining non-equilibrium steady states.

The formalism quantifies the entropic cost of diffusive energy transport.

Abstract

When a quantum system is coupled to several heat baths at different temperatures, it eventually reaches a non-equilibrium steady state featuring stationary internal heat currents. These currents imply that entropy is continually being produced in the system at a constant rate. In this paper we apply phase-space techniques to the calculation of the Wigner entropy production on general linear networks of harmonic nodes. Working in the ubiquitous limit of weak internal coupling and weak dissipation, we obtain simple closed-form expressions for the entropic contribution of each individual quasi-probability current. Our analysis highlights the essential role played by the internal unitary interactions (node-node couplings) in maintaining a non-equilibrium steady-state and hence a finite entropy production rate. We also apply this formalism to the paradigmatic problem of energy transfer through a chain of oscillators subject to self-consistent internal baths that can be used to tune the transport from ballistic to diffusive. We find that the entropy production scales with different power law behaviors in the ballistic and diffusive regimes, hence allowing us to quantify what is the "entropic cost of diffusivity."

Equations

$$\frac{\mathrm{d}S(t)}{\mathrm{d}t}=\Pi (t)-\Phi (t),$$

$$\frac{\mathrm{d}\mathbf{\rho }}{\mathrm{d}t}=\mathcal{U}(\mathbf{\rho })+{\sum }_{k=1}^L{\mathcal{L}}_k(\mathbf{\rho })+\text{pazocal}O({\gamma }^2).$$

$$\frac{\mathrm{d}\mathbf{\rho }}{\mathrm{d}t}=-i[\mathbf{H},\mathbf{\rho }]+{\sum }_{k=1}^L{\mathcal{D}}_k(\mathbf{\rho })+\text{pazocal}O(\lambda \gamma ).$$

$$\begin{gathered} {\mathcal{D}}_k(\mathbf{\rho })={\gamma }_k(n_k+1)\left({\mathbf{a}}_k\mathbf{\rho }{\mathbf{a}}_k^{†}-\frac{1}{2}\{{\mathbf{a}}_k^{†}{\mathbf{a}}_k,\mathbf{\rho }{\}}_{+}\right) \\ +{\gamma }_k{n}_k\left({\mathbf{a}}_k^{†}\mathbf{\rho }{\mathbf{a}}_k-\frac{1}{2}\{{\mathbf{a}}_k{\mathbf{a}}_k^{†},\mathbf{\rho }{\}}_{+}\right), \end{gathered}$$

$$W(\bar{\alpha }):=\frac{1}{{\pi }^{2L}}\int \mathrm{d}\bar{\beta }{e}^{-{\sum }_k({\beta }_k{\alpha }_k^{\ast }-{\beta }_k^{\ast }{\alpha }_k)}\mathrm{tr}\left\{\mathbf{\rho }\mathbf{D}(\bar{\beta })\right\},$$

$${\partial }_{t}W=\text{pazocal}U(W)+{\sum }_{k=1}^L\text{pazocal}{D}_k(W),$$

$$\text{pazocal}{D}_k(W)={\partial }_k\text{pazocal}{J}_k(W)+{\partial }_k^{\ast }\text{pazocal}{J}_k^{\ast }(W),$$

$$\text{pazocal}{J}_k(W)=\frac{{\gamma }_k}{2}\left[{\alpha }_{k}W+\left(n_k+1/2\right){\partial }_k^{\ast }W\right].$$

$$W_{\text{eq}}={\prod }_{k=1}^L{W}_{\text{eq}}^{(k)},W_{\text{eq}}^{(k)}=\frac{e^{-|{\alpha }_k{|}^2/(n_k+1/2)}}{\pi (n_k+1/2)}.$$

$$\text{pazocal}U(W)={\sum }_{k=1}^L\left[{\partial }_k\text{pazocal}{A}_k(W)+{\partial }_k^{\ast }\text{pazocal}{A}_k^{\ast }(W)\right],$$

$$\text{pazocal}{A}_k(W)=i{\sum }_{\ell =1}^L{H}_{k\ell }{\alpha }_{\ell }W.$$

$$S(W):=-\int \mathrm{d}\bar{\alpha }W(\bar{\alpha })\ln W(\bar{\alpha }).$$

$$\frac{\mathrm{d}S}{\mathrm{d}t}=-\int \mathrm{d}\bar{\alpha }\left[\text{pazocal}U(W)+{\sum }_k\text{pazocal}{D}_k(W)\right]\ln W.$$

$$\frac{\mathrm{d}S}{\mathrm{d}t}$$

$$-{\sum }_k\int \mathrm{d}\bar{\alpha }\text{pazocal}{D}_k(W)\ln W_{\text{eq}}:=\Pi -\Phi ,$$

$$\Pi ={\sum }_k\int \mathrm{d}\bar{\alpha }\left[\text{pazocal}{J}_k(W){\partial }_k+\text{pazocal}{J}_k^{\ast }(W){\partial }_k^{\ast }\right]\ln (W/W_{\text{eq}}).$$

$$\text{pazocal}{J}_k(W)=\frac{{\gamma }_k}{2}(n_k+1/2)W{\partial }_k^{\ast }\ln (W/W_{\text{eq}}),$$

$${\Pi }_k=\frac{4}{{\gamma }_k(n_k+1/2)}\int \mathrm{d}\bar{\alpha }\frac{{\left|\text{pazocal}{J}_k(W)\right|}^2}{W}.$$

$${\Phi }_k=-\int \mathrm{d}\bar{\alpha }\left[\text{pazocal}{J}_k(W){\partial }_k\ln W_{\text{eq}}+\text{pazocal}{J}_k^{\ast }(W){\partial }_k^{\ast }\ln W_{\text{eq}}\right].$$

$$\int \mathrm{d}\bar{\alpha }{\alpha }_k^{\ast }{\alpha }_{\ell }W=⟨{\mathbf{a}}_k^{†}{\mathbf{a}}_{\ell }⟩+\frac{1}{2}{\delta }_{k,\ell },$$

$${\Phi }_k=\frac{{\gamma }_k}{n_k+1/2}\left(⟨{\mathbf{a}}_k^{†}{\mathbf{a}}_k⟩-n_k\right).$$

$$S(W\Vert W_{\text{eq}}):=\int \mathrm{d}\bar{\alpha }W\ln (W/W_{\text{eq}}),$$

$$\frac{\mathrm{d}}{\mathrm{d}t}S(W\Vert W_{\text{eq}})$$

$$+{\sum }_{k=1}^L\int \mathrm{d}\bar{\alpha }\text{pazocal}{D}_k(W)\ln (W/W_{\text{eq}})$$

$$=-\int \mathrm{d}\bar{\alpha }\text{pazocal}U(W)\ln W_{\text{eq}}-\Pi ,$$

$$\Pi (t\to \infty ):={\Pi }_{\text{NESS}}=-\int \mathrm{d}\bar{\alpha }\text{pazocal}U(W_{\text{NESS}})\ln W_{\text{eq}}.$$

$${\Pi }_{\text{NESS}}$$

$$={\sum }_{k\ne \ell }\frac{2}{n_k+1/2}\text{Im}\{H_{k,\ell }⟨a_k^{†}{a}_{\ell }⟩\}.$$

$$\frac{\mathrm{d}⟨{\mathbf{a}}_k^{†}{\mathbf{a}}_k⟩}{\mathrm{d}t}$$

$$j_{k,\ell }$$

$${\Pi }_{\text{NESS}}=\frac{1}{2}{\sum }_{k\ne \ell }{j}_{k,\ell }\left(\frac{1}{n_k+1/2}-\frac{1}{n_{\ell }+1/2}\right).$$

$$\mathbf{H}=\omega {\sum }_{k=1}^L{\mathbf{a}}_k^{†}{\mathbf{a}}_k+i\lambda {\sum }_{k=1}^{L-1}\left({\mathbf{a}}_k^{†}{\mathbf{a}}_{k+1}-{\mathbf{a}}_{k+1}^{†}{\mathbf{a}}_k\right).$$

$$\frac{\mathrm{d}\mathbf{\rho }}{\mathrm{d}t}=-i[\mathbf{H},\mathbf{\rho }]+{\mathcal{D}}_1(\mathbf{\rho })+{\mathcal{D}}_L(\mathbf{\rho })+{\sum }_{k=1}^L{\tilde{\mathcal{D}}}_k(\mathbf{\rho }).$$

Full text

I Introduction

Entropy plays a fundamental role in both thermodynamics and information theory. Unlike energy, entropy does not satisfy a continuity equation—in addition to the exchange of entropy with the environment, it can also be produced within the system. This additional contribution is known as entropy production and serves to gauge the irreversibility of a physical process Callen (1985); Tomé and Oliveira (2015). Indeed, one can write the following balance equation for the entropy $S(t)$ of a system

[TABLE]

where $\Pi(t)$ is the entropy-production rate and $\Phi(t)$ stands for the entropy flux from the system into the environment. As a consequence of the second law of thermodynamics, $\Pi(t)$ is always non-negative and vanishes if and only if the system is in equilibrium. When the system is allowed to relax in contact with a single heat bath, it will generally reach thermal equilibrium, so that $\,\mathrm{d}S/\,\mathrm{d}t=\Pi=\Phi=0$. However, when it is connected to multiple baths at different temperatures, it will instead reach a non-equilibrium steady-state (NESS) where $\,\mathrm{d}S/\,\mathrm{d}t=0$ but $\Pi=\Phi\geq 0$. The NESS is therefore characterized by a finite entropy production rate $\Pi$, which is continuously converted into an entropy flux $\Phi$ and dumped into the environment.

The theory of entropy production is formulated differently depending on the type of stochastic process at hand. For classical systems, approaches based on Onsager’s theory of chemical kinetics Onsager (1931); Machlup and Onsager (1953); Tisza and Manning (1957), classical master equations Schnakenberg (1976); Tomé and de Oliveira (2012), or Fokker–Planck equations Tomé and de Oliveira (2010); Spinney and Ford (2012); Landi et al. (2013); Van den Broeck and Esposito (2010); Seifert (2012) have been widely used. Conversely, for quantum systems the problem is usually formulated in terms of Gorini–Kossakowski–Lindblad–Sudarshan (GKLS) quantum master equations, Gorini (1976); Lindblad (1976); Spohn (1978); Leggio et al. (2013); Breuer and Petruccione (2007), repeated interactions Scarani et al. (2002); Gennaro et al. (2008); Burgarth and Giovannetti (2007); Bodor et al. (2013); Ziman et al. (2005); Uzdin and Kosloff (2014); Barra (2015); Barra and Lledó (2017); Strasberg et al. (2017); De Chiara et al. (2018), quantum trajectories Elouard et al. (2017), and fluctuation theorems Jarzynski (1997); Crooks (1998), among others.

More recently, a formulation combining quantum phase-space methods and Fokker–Planck equations has been put forward Santos et al. (2017); Brunelli et al. (2018); Santos et al. (2018a, b) and applied to general non-thermal environments, such as squeezed, dephasing and even photon loss reservoirs, that are often used in the description of optical cavities and input-output theory Gardiner and Zoller (2004). This is in fact essential if one wishes to have a complete assessment of entropy production in controlled quantum experiments. In addition, this framework allows to identify irreversible quasi-probability currents in phase space, which are the elementary contributions ultimately responsible for the emergence of irreversibility at the quantum level. However, the efforts so far have been focused exclusively on systems connected to a single reservoir and have not addressed the phenomenology of NESSs. In this paper, we set out to fill this gap.

In particular, we study generic lattice systems of linearly coupled harmonic nodes connected to various heat baths. To describe the ensuing dissipative dynamics, we adopt a local master equation Wichterich et al. (2007); Levy and Kosloff (2014)—which is accurate to lowest order in the inter-node interaction strength Trushechkin and Volovich (2016)—and exploit its simple structure. We obtain closed-form analytical expressions for the steady-state irreversible entropy production and entropy flux, broken down into the elementary contributions corresponding to individual dissipation channels. In addition, we also obtain an expression which neatly illustrates the essential role played by the unitary dynamics in sustaining the irreversible entropy production, by enabling energy transport through the lattice.

As an application, we study diffusive heat transport through a harmonic chain connected to two reservoirs at each end Žnidarič (2010); Manzano et al. (2012); Asadian et al. (2013); Santos and Landi (2016); Karevski and Platini (2009). In addition, in order to switch between ballistic and diffusive behavior, we introduce auxiliary self-consistent reservoirs Bolsterli et al. (1970). That is, additional baths that inject noise in the chain without affecting the energy transport, but nonetheless causing decoherence. As is well known, this mechanism turns an a priori ballistic transport, typical of harmonic lattices, into a diffusive one Bolsterli et al. (1970); Bonetto et al. (2004); Pereira and Falcao (2004); Delfini et al. (2006); Dhar and Roy (2006). Self-consistent baths lead to a NESS current similar to that of dephasing baths, as studied for instance in Žnidarič (2010); Manzano et al. (2012); Asadian et al. (2013); Nicacio et al. (2015). However, they have the advantage that they lead to a Gaussian NESS (which is not the case for dephasing). Using our results, we are then able to split the total entropy production rate into a contribution stemming from the physical reservoirs at the boundaries and a second one, solely due to the self-consistent baths. We observe that each contribution scales polynomially with a distinct exponent. This then allows us to unambiguously associate the later with the entropic cost of maintaining a diffusive profile in the NESS.

This paper is structured as follows: We begin by describing our weakly interacting harmonic lattice systems and the phase-space representation of the underlying (local) master equation in Sec. II. We then introduce the central object of our analysis—the Wigner entropy production—and carry out its explicit calculation and discuss the rich physics underlying its neat decomposition into elementary contributions (cf. Sec. III). Finally, Sec. IV discusses the application to estimating the entropic cost of diffusivity. Conclusions are summarized in Sec. V and, in order to make this paper self-contained, we also include a series of Appendices (A–E) providing further details on the calculations outlined in Sec. IV.

II Quantum phase-space methods

II.1 The model and the quantum master equation

We consider here a linear network of $L$ harmonic oscillators with Hamiltonian $\boldsymbol{H}\coloneqq\boldsymbol{H}_{0}+\boldsymbol{H}_{\text{int}}$, where the local part is simply $\boldsymbol{H}_{0}\coloneqq\sum\nolimits_{k=1}^{L}H_{kk}\,\boldsymbol{a}_{k}^{\dagger}\boldsymbol{a}_{k}$ and $\boldsymbol{H}_{\text{int}}=\sum_{k\neq\ell}H_{k,\ell}\,\boldsymbol{a}_{k}^{\dagger}\boldsymbol{a}_{\ell}$ stands for the (linear) inter-node couplings. Here, $\boldsymbol{a}_{i}$ denotes the bosonic annihilation operator of the $i$-th mode and, in what follows, $\hbar=k_{B}=1$. Otherwise, we impose no restrictions on the structure of the network, shaped by the non-zero off-diagonal matrix elements $H_{\ell,k}=H_{k,\ell}^{*}$ ($H_{kk}=\omega_{k}$).

If we let each mode interact weakly with a local heat bath, the resulting dissipative dynamics can generally be described by a Gorini, Kossakowski, Lindblad and Sudarshan (GKLS) quantum master equation of the form Breuer and Petruccione (2007)

[TABLE]

Here, $\mathscr{U}(\boldsymbol{\rho})\coloneqq-i[\boldsymbol{H},\boldsymbol{\rho}]$ and $\mathscr{L}_{k}$ stands for a dissipation super-operator in the standard GKLS form Gorini (1976); Lindblad (1976). Moreover, $\gamma\coloneqq\max\{\gamma_{1},\cdots,\gamma_{L}\}$ equals the largest of the node–bath coupling strengths $\gamma_{k}$, i.e., it carries the order of magnitude of the dissipative interactions. Individually, each super-operator $\mathscr{L}_{k}$ is $\pazocal{O}(\gamma_{k})$. In general, any of the “dissipators” $\mathscr{L}_{k}$ can act globally on all nodes of the network. One may, however, make the additional assumption of weak internal coupling between the nodes, which leads to the local GKLS quantum master equation Trushechkin and Volovich (2016)

[TABLE]

where the local dissipators are given by

[TABLE]

with $n_{k}\coloneqq(e^{\omega_{k}/T_{k}}-1)^{-1}$ being the Bose–Einstein thermal occupation of at temperature $T_{k}$ and frequency $\omega_{k}$, and where $\{\cdot,\cdot\}_{+}$ stands for anti-commutator.

Care must be taken when using Eq. (3) to describe a multipartite open quantum system, since it is well known that going beyond its range of validity might lead to thermodynamic inconsistencies Levy and Kosloff (2014); Stockburger and Motz (2016); Trushechkin and Volovich (2016). The parameter range in which the many approximations underlying both global and local master equations are satisfied have been critically (and extensively) discussed in the literature Rivas et al. (2010); Wichterich et al. (2007); Karevski and Platini (2009); Purkayastha et al. (2016); González et al. (2017); Hofer et al. (2017); Seah et al. (2018); Kołodyński et al. (2018); Mitchison and Plenio (2018); Naseem et al. (2018); Seah and Nimmrichter (2018). Interestingly, local master equations also arise naturally when considering certain repeated-interaction models Karevski and Platini (2009); Barra (2015); Seah and Nimmrichter (2018). This picture enables a thermodynamically consistent bookkeeping of all energy exchanges occurring in dissipative processes exactly described by Eq. (3) Barra and Lledó (2017); De Chiara et al. (2018). In this paper, we remain on the safe side by working well within the range of applicability of Eq. (3), understood as a mere perturbative expansion of Eq. (2) to lowest order in the internal couplings.

II.2 The Fokker–Planck equation

We now move from Hilbert space to the quantum phase space. Introducing the displacement operator \boldsymbol{D}(\bar{\beta})\coloneqq\exp{\big{(}\sum\nolimits_{k=1}^{L}\beta_{k}\boldsymbol{a}_{k}^{\dagger}-\beta_{k}^{*}\boldsymbol{a}_{k}\big{)}}, the Wigner function of our $L$-mode Gaussian state $\boldsymbol{\rho}$ can be written as Gardiner and Zoller (2004)

[TABLE]

where $\bar{\alpha}$ and $\bar{\beta}$ are $2L$-dimensional complex vectors containing phase-space variables, i.e., $\bar{\alpha}\coloneqq(\alpha_{1},\alpha_{1}^{*},\ldots,\alpha_{L},\alpha_{L}^{*})^{\top}$. The reason for choosing the Wigner representation lies in the Gaussianity of $\boldsymbol{\rho}$, which brings at least three advantages. First, it is always strictly positive, which allows us to identify the Wigner function as a quasi-probability distribution. The second, it will depend only on the first and second order moments, which greatly simplifies the analysis (cf. Appendix A for details). Finally, the corresponding Wigner entropy, to be introduced below, can be related to the Rényi-2 entropy and also satisfies the strong subadditivity inequality Adesso et al. (2012), hence giving it a more physical appeal.

Using standard correspondence tables (see, e.g., Gardiner and Zoller (2004)) one can turn Eq. (3) into the following Fokker–Planck equation

[TABLE]

where the phase-space super-operator $\pazocal{U}(W)$ represents the unitary part of Eq. (3), while $\pazocal{D}_{k}(W)$ stands for the dissipative contributions. In particular, the latter may be written as divergences in the complex plane; namely

[TABLE]

where $\partial_{k}\coloneqq\partial/\partial\alpha_{k}$ ($\partial_{k}^{*}\coloneqq\partial/\partial\alpha_{k}^{*}$) and

[TABLE]

Eq. (6) can be interpreted as a continuity equation for $W$ with $\pazocal{J}_{k}(W)$ playing the role of irreversible quasi-probability currents in phase space. In fact, these currents are identically zero if and only if each oscillator is in local thermal equilibrium 111This follows directly from the fact that $\boldsymbol{\rho}_{\text{eq}}=\bigotimes_{k}\,\boldsymbol{\rho}_{\text{eq}}^{(k)}$, with \boldsymbol{\rho}_{\text{eq}}^{(k)}=\exp{\big{(}-\omega_{k}\,\boldsymbol{a}_{k}^{\dagger}\boldsymbol{a}_{k}/T_{k}\big{)}}, is the unique fixed point of the dissipative part of the local master equation (3). Notice, however, that $\,\mathrm{d}\boldsymbol{\rho}_{\text{eq}}/\,\mathrm{d}t\neq 0$; not even within the range of validity of the local approximation. Specifically, $\,\mathrm{d}\boldsymbol{\rho}_{\text{eq}}/\,\mathrm{d}t=-i[\boldsymbol{H},\boldsymbol{\rho}_{\text{eq}}]=-i[\boldsymbol{H}_{\text{int}},\boldsymbol{\rho}_{\text{eq}}]=\pazocal{O}(\lambda)\gg\pazocal{O}(\lambda\,\gamma)$.; i.e., $\pazocal{J}_{k}(W_{\text{eq}})=0~{}\forall k$ with

[TABLE]

Since it is only in $W_{\text{eq}}$ that all the individual currents $\pazocal{J}_{k}(W)$ vanish exactly, we shall adopt it as the reference state for quantifying entropy-production rates and fluxes Santos et al. (2017). In this limited sense, $W_{\text{eq}}$ would correspond to the thermal equilibrium state in the standard formulation.

Similarly, the unitary part $\pazocal{U}$ in Eq. (6) can be cast as

[TABLE]

where, for our choice of Hamiltonian, the reversible quasi-probability currents $\pazocal{A}_{k}(W)$ are simply

[TABLE]

III Wigner entropy-production rate

III.1 Individual dissipative contributions

We shall now decompose $\,\mathrm{d}S/\,\mathrm{d}t$ as in Eq. (1) into an entropy-production rate and the entropy flux. We base our analysis in the Wigner entropy defined as

[TABLE]

As shown in Refs. Santos et al. (2017); Brunelli et al. (2018); Santos et al. (2018a, b), for Gaussian systems this is entirely equivalent to the standard approach based on von Neumann entropy. Besides, the Wigner entropy offers several advantages, as we shall see below.

Differentiating Eq. (12) with respect to time and using Eq. (6) yields

[TABLE]

Adding and subtracting $\sum\nolimits_{k}\int\,\mathrm{d}\bar{\alpha}\,\pazocal{D}_{k}(W)\ln{W_{\text{eq}}}$, so as to introduce our fiducial state $W_{\text{eq}}$, gives

[TABLE]

where we have used the fact that the unitary dynamics $\pazocal{U}$ does not change the Wigner entropy. The identification of the two terms in this equation as an entropy production rate and an entropy flux rate will be better justified once they are evaluated explicitly, as we shall now do.

Inserting Eq. (7) in $\Pi$, we can integrate by parts in each variable. Noticing that the corresponding boundary terms vanish due to the Gaussianity of $W$ and $W_{\text{eq}}$, we get

[TABLE]

Next, we rewrite the currents $\pazocal{J}_{k}(W)$ in Eq. (8) as

[TABLE]

which allows us to express $\partial_{k}^{*}\ln{(W/W_{\text{eq}})}$ in terms of $\pazocal{J}_{k}(W)$ and $W$. This leads to the decomposition $\Pi=\sum\nolimits_{k}\,\Pi_{k}$, where

[TABLE]

As shown in Appendix C, further manipulations lead to a compact expression in terms of the covariances of the system. Similarly, the entropy flux decomposes as $\Phi=\sum\nolimits_{k}\Phi_{k}$ with

[TABLE]

Making use of the identity

[TABLE]

we finally get the compact expression

[TABLE]

Eqs. (17) and (20) constitute our first main result. Their identification as entropy-production rate and entropy flux, respectively, is based on several supporting arguments and has been extensively debated in the past, both in the quantum Santos et al. (2017) and classical Seifert (2012); Tomé and de Oliveira (2010) contexts. First and foremost, $\Pi$ is clearly non-negative and zero if and only if the currents themselves vanish, which only occurs in the reference state $W_{\text{eq}}$. Second, the proposed entropy-production rate is an even function of the irreversible currents whereas the entropy flux is odd, in agreement with other studies based on Fokker–Planck equations Tomé and de Oliveira (2010). Finally and most remarkably, it can be shown that, within the framework of stochastic trajectories, this expression satisfies integral fluctuation theorems Santos et al. (2017); Spinney and Ford (2012).

Rather intuitively, we can also see that Eq. (20) is proportional to the difference between the actual occupation of the $k$-th mode and the equilibrium occupation at the temperature $T_{k}$ of the local bath. Namely, if the mode “looks hotter” than its local environment, one has $\Phi_{k}>0$, i.e., entropy flows from the system into the bath, as expected. Conversely, if $\langle\boldsymbol{a}_{k}^{\dagger}\boldsymbol{a}_{k}\rangle<n_{k}$, $\Phi_{k}<0$ and entropy flows into the system. Note that $\Phi_{k}$ is, therefore, an inherently observable quantity Santos et al. (2017).

For us, the most relevant feature of Eqs. (17) and (20) is precisely that they provide the individual contribution of each dissipation channel to the the total entropy-production rate and entropy flux. In particular, from a classical viewpoint Seifert (2012) one may understand Eq. (17) as an average of the “phase-space velocity” $\pazocal{J}_{k}(W)/W$ of each individual quasi-probability current $\pazocal{J}_{k}(W)$. Hence, the entropy production can be thought-of as weighted average of such mean velocities.

III.2 Role of the unitary dynamics in maintaining a NESS

Let us introduce the Wigner relative entropy (or Kullback–Leibler divergence), defined as

[TABLE]

and calculate its rate of change using Eq. (6). This yields

[TABLE]

where we have used Eq. (III.1). In the long-time limit, the stationarity condition $\,\mathrm{d}S(W_{\text{NESS}}\;\|\;W_{\text{eq}})/\,\mathrm{d}t=0$ entails

[TABLE]

This is our second main result: it illustrates the essential role of the unitary internal dynamics in sustaining stationary heat currents across the network which, in turn, translates into a finite rate of steady-state entropy production.

Inserting Eq. (10) for the unitary currents and using (19) we also find that the last term in Eq. (23) may be written as

[TABLE]

We point to the fact that $\Pi_{\text{NESS}}$ in Eq. (24) contains exclusively off-diagonal elements $H_{k\neq\ell}$ and none of the local components $H_{kk}\,\boldsymbol{a}_{k}^{\dagger}\boldsymbol{a}_{k}$; it is thus entirely due to the node–node couplings. Introducing the energy currents $j_{k,\ell}$ defined by Karevski and Platini (2009)

[TABLE]

allows to cast Eq. (24) in the symmetric form

[TABLE]

This result can be connected with Onsager’s theory of irreversible thermodynamics Callen (1985); Onsager (1931). Within this framework, the entropy production is defined as the product of ‘fluxes’ and ‘affinities’ (also called ‘generalized forces’). For instance, the current of energy is related to the affinity $1/T$ so that, in a classical scenario, the Onsager entropy production between two bodies kept at temperatures $T_{A}$ and $T_{B}$ is given by \Pi=j_{AB}\,\big{(}1/T_{A}-1/T_{B}\big{)}, where $j_{AB}$ stands for the energy current from $B$ to $A$. We see that Eq. (26) has the exact same mathematical structure, which is yet another consistency test for our analysis. Moreover, notice that, due to the fact that we are using Wigner entropies, the thermodynamic affinity related to the current is not the inverse temperature, but rather the inverse Bose–Einstein occupation $n+1/2$. For high temperatures $n+1/2\propto T$ so that both frameworks coincide. However, our results hold true even in the limit of vanishingly low temperature222Always provided that the Markov approximation behind the master equation (3) continues to hold. This is often seen as a high-temperature limit.. We also remark that Onsager’s formula is valid only close to equilibrium (linear response theory), whereas Eq. (26) holds true for arbitrary non-equilibrium states. This is a consequence of the Gaussianity of the problem.

Finally note that, as $\Pi_{\text{NESS}}\geq 0$, Eq. (26) implies that $j_{k,\ell}$ has the same sign as $1/(n_{k}+1/2)-1/(n_{\ell}+1/2)$; i.e., when node $k$ is hotter than node $\ell$, the local current $j_{k,\ell}$ flows from $k$ to $\ell$.

IV Application: The entropic cost of diffusivity

In order to illustrate the usefulness of this framework, we now apply it to a typical quantum transport problem. In particular, we show how to exploit it to quantify the irreversibility associated with sustaining diffusive heat transfer through a chain of harmonic oscillators.

Before describing the model further, let us briefly discuss the notions of ‘ballistic’ and ‘diffusive’ transport. Consider the NESS established when a system (of “length” $L$) is placed between two reservoirs at different temperatures. This exhibits a stationary heat current which, for small temperature gradients $\Delta T$, may be fitted by $j=\kappa\,\Delta T/L^{\alpha}$, where $\alpha$ is an exponent that characterizes the phenomenological heat transport law. Whenever $\alpha=1$ (i.e., Fourier’s heat conduction law) one speaks of diffusive behavior. However, low-dimensional—and usually integrable—quantum and classical systems can display a ballistic heat flow Rieder et al. (1967); Landi and de Oliveira (2013); Dhar (2008), characterized by $\alpha=0$; that is, a heat current that is independent of the system size.

Here, we aim at understanding how the entropy-production rate in the steady state relates to the heat transfer law; specifically, considering the divide between diffusive and ballistic heat conduction. To do so, we must adopt a model allowing to switch between these two regimes. In particular, our one-dimensional harmonic chain coupled to two different reservoirs at each end would result in a ballistic NESS, provided that its dissipative dynamics is well described by Eq. (3). Obtaining diffusive behavior in general requires complex interactions, which are seldom treatable analytically. An alternative, widely used in the literature on quantum systems, is to use dephasing baths acting locally on each site Žnidarič (2010); Manzano et al. (2012); Asadian et al. (2013). These baths introduce some noise while not compromising the energy balance. As a consequence, they always lead to diffusive transport.

The NESS for a dephasing model, however, is not Gaussian, so that the results presented here would not apply. Instead, we fix this problem by using the notion of ‘self consistent reservoirs’ Bolsterli et al. (1970); Žnidarič (2010); Manzano et al. (2012); Asadian et al. (2013); Nicacio et al. (2015) (also called Büttiker probes in the condensed matter literature). As the name suggests, self-consistent reservoirs are baths that act on all sites, but whose temperature are chosen self-consistently so that, in the steady state, they do not exchange any energy with the system, but only inject noise to cause decoherence Bolsterli et al. (1970). Consequently, heat only flows between the two physical reservoirs at the ends of the chain (see Fig. 1). As shown in Appendix D, the resulting equation of motion for the covariances is the same for self-consistent and dephasing models although the corresponding NESSs are different (i.e., they only agree up to the second-order moments).

IV.1 The model and its stationary solution

In what follows, we consider $L$ resonant harmonic oscillators in a one-dimensional lattice interacting weakly via the linear Hamiltonian

[TABLE]

As described above, the system is coupled to $L+2$ baths, out of which two are physical and the rest, auxiliary self-consistent reservoirs. The master equation (3) then reads

[TABLE]

Here all dissipators have the GKLS structure of Eq. (4). Those of the physical baths have parameters $\gamma_{1}=\gamma_{L}=\gamma$ and the corresponding occupations at frequency $\omega$ are $n_{1}$ and $n_{L}$. On the other hand, the auxiliary baths all have dissipation rate $\tilde{\gamma}_{k}=\Gamma$ and occupations $\tilde{n}_{k}$, calculated self-consistently so as to satisfy $\tilde{n}_{k}=\operatorname{tr}\{\boldsymbol{a}_{k}^{\dagger}\boldsymbol{a}_{k}\,\boldsymbol{\rho}_{\text{NESS}}\}$. All objects denoted with tilde refer to the auxiliary baths. In particular, notice the difference between $\tilde{\mathscr{D}}_{1,L}(\boldsymbol{\rho})$, or $\tilde{n}_{1,L}$, and $\mathscr{D}_{1,L}(\boldsymbol{\rho})$ and $n_{1,L}$.

The stationary state of Eq. (28) (i.e., $\boldsymbol{\rho}_{\text{NESS}}$) can be found analytically with the methods developed in Refs. Žnidarič (2010); Manzano et al. (2012); Asadian et al. (2013); Nicacio et al. (2015). For completeness, we give details of the calculation in Appendix D. We find that the only non-zero stationary covariances in our chain are

[TABLE]

Eqs. (29) have a rich physical interpretation. First, setting $\Gamma=0$ amounts to decoupling the self-consistent reservoirs from the chain. As already anticipated, in that case the occupations $\langle\boldsymbol{a}_{k}^{\dagger}\boldsymbol{a}_{k}\rangle$ become independent of the index $k$ (except for $k=1$ and $k=L$), which is the hallmark of a ballistic heat transfer law [see Fig. 2(a)]. Conversely, if we choose $\Gamma>0$ and sufficiently large $L$, the occupations converge to a linear profile, which is expected for diffusive heat conduction in one dimension. This is illustrated in Figs. 2(b)–(d).

Furthermore, combining Eqs. (25b) and (29b) we can compute the current between neighboring sites

[TABLE]

which are all independent of $k$. Here we see once again that if $\Gamma=0$, the energy current becomes independent of $L$, implying a ballistic transport. On the other hand, $\Gamma>0$ yields a diffusive heat current, scaling as $\sim 1/L$. We remark that discussing the divide between ballistic and diffusive conduction is only meaningful in the thermodynamic limit (i.e., $L\gg 1$), where our model suffers an abrupt transition from a ballistic profile to a diffusive one as soon as $\Gamma$ becomes non-zero.

IV.2 Entropy-production rate in the NESS

The total (Wigner) entropy-production rate can be readily computed from Eqs. (26) and (30). Since the currents $j_{k,k+1}$ are translationally invariant, all inner terms in the sum of Eq. (26) vanish and we are left with

[TABLE]

Note that this is always non-negative and zero only if $n_{1}=n_{L}$, in which case there is no current flowing through the chain. We also remark that neither $j$ nor $\Pi_{\text{NESS}}$ depend on $\tilde{n}_{k}$. This is an important consistency check, as our goal is precisely to ensure that the auxiliary reservoirs do not modify the rate of heat flow across the chain. Notice as well that $\Pi_{\text{NESS}}$ scales in size just like the rate of heat exchange $j$.

Importantly, exploiting Eq. (17) we are also able to separate the contributions of each dissipation channel to the stationary entropy-production rate. More specifically, we can split $\Pi_{\text{NESS}}$ into one component stemming from the real reservoirs at the boundaries and another one coming from the auxiliary self-consistent reservoirs; that is, $\Pi_{\text{NESS}}\coloneqq\Pi_{\text{r}}+\Pi_{\text{sc}}$. Crucially, due to our approach to emulate the anharmonicity leading to the diffusive profile, we can meaningfully identify the self-consistent contribution $\Pi_{\text{sc}}$ with the entropic cost of steady-state diffusivity. Specifically, one has

[TABLE]

while the entropy production from the self-consistent reservoirs reads

[TABLE]

where we already used the fact that the occupations of the self-consistent reservoirs are $\tilde{n}_{k}=\langle\boldsymbol{a}_{k}^{\dagger}\boldsymbol{a}_{k}\rangle$. The sum in Eq. (32) comprises the terms $k=1$ and $k=L$ only, while that of Eq. (33a) runs over all indices $k\in\{1,\cdots,L\}$. As discussed in Appendix E, Eqs. (32) and (33a) can be readily computed from the covariances of the chain for arbitrarily large system size, owing to the Gaussianity of our NESS.

In Fig. 3 we plot $\Pi_{\text{r}}$ and $\Pi_{\text{sc}}$ as a function of $L$ in a logarithmic scale. Interestingly, we see that the irreversibility associated with the physical reservoirs is dominant only for small system size. As the chain is scaled up $\Pi_{\text{sc}}$ quickly surpasses $\Pi_{\text{r}}$ and adopts a distinct power-law-like decay $\Pi_{\text{sc}}\sim 1/L$ for large $L$. On the other hand—and also in the thermodynamic limit—we observe that the contribution from the real baths to the total steady-state irreversibility decays as $\Pi_{\text{r}}\sim 1/L^{2}$. Looking back at Eq. (31), it is also noteworthy that, fixing the temperature gradient, the ballistic irreversible entropy-production rate $\Pi_{\text{NESS}}(\gamma,\Gamma=0)$ is larger than the diffusive one $\Pi_{\text{NESS}}(\gamma,\Gamma>0)$, meaning that $\Pi_{\text{r}}(\gamma,\Gamma=0)>\Pi_{\text{r}}(\gamma,\Gamma\neq 0)+\Pi_{\text{sc}}(\gamma,\Gamma\neq 0)$. The dependence on both $\gamma$ and $\Gamma$, although not explicit in Eqs. (32) and (33), arises from $W_{\text{NESS}}$.

We have thus seen that combining the ease of calculation of the Wigner entropy production for Gaussian states with our simple harmonic model to mimic diffusive heat conduction, provides valuable insights into the relative weight of the individual irreversible processes in the thermodynamic limit.

V Conclusions

We have addressed the calculation of the irreversible entropy-production rate using a quantum phase-space approach based on the Wigner entropy. In particular, we focused on networks of weakly interacting harmonic nodes with arbitrary connectivity, and coupled to various reservoirs at different temperatures. For this wide class of systems, we were able to obtain simple and useful closed-form expressions for both the Wigner entropy-production rate and the entropy flux, solely in terms of the second-order moments of the system. This is possible since we work with an overall harmonic Hamiltonian, which preserves Gaussianity. In addition, we could split the entropy-production rate and flux into contributions stemming from individual quasi-probability currents associated with each open decay channel, which enables one to identify the irreversibility generated by a single bath in the lattice. We also discussed how the internal coherent dynamics plays a central role in generating steady-state irreversibility, as it is the leading mechanism enabling energy transport across the network. All of the above should be understood as the leading (lowest-order) contribution in a perturbative expansion in the node-to-node interaction strengths (i.e., the limit of weak internal coupling), as it relies on the simple local structure of the underlying quantum master equation. Finally, we used our framework to better understand the interplay between the various sources of irreversibility at play in diffusive heat conduction through a harmonic chain. We mimicked the anharmonicity required to establish the desired diffusive profile by adding auxiliary (self-consistent) reservoirs to our model. In turn, this allowed us to break down the total steady-state irreversible entropy production into a contribution due to the heat transport across the chain, and another one which can be interpreted as the entropic cost of maintaining a stationary diffusive transport.

As we have shown, this approach offers several advantages, both practical and fundamental. From a practical standpoint, the fact that Gaussian states can be fully characterized by their covariance matrix allows one to readily access any entropic quantifier, even for very large system sizes. Secondly, from a fundamental point of view, our approach allows for a microscopic description of the problem of entropy production, enabling the identification of irreversible quasi-probability currents in phase space which are ultimately responsible for the emergence of irreversibility. Moreover, the description in terms of the Wigner function also allows us to take into account both thermal and quantum fluctuations, remaining valid even in the limit of very low temperatures, where the standard von Neumann approach becomes problematic. More generally, our framework can be readily adapted to quantify irreversibility in non-trivial energy-conversion processes (e.g., refrigeration) implemented on small-scale quantum heat devices. Indeed, networks of periodically driven and weakly interacting oscillators are a suitable platform to generate the required non-equilibrium states supporting useful continuous quantum-thermodynamic cycles.

Acknowledgments

The authors are grateful to T. Tufarelli for useful comments. WTBM and GTL acknowledge the financial support from the São Paulo research foundation (FAPESP), under grants 2016/08721-7, 2017/50304-7, 2017/06323-7, and 2018/08211-4. LAC acknowledges funding from the European Research Council (StG GQCOP No. 637352). GTL acknowledges the financial support from the Brazilian funding agency CNPq (grant number INCT-IQ 246569/2014-0). MP acknowledges the H2020 Collaborative Project TEQ (grant agreement nr. 766900), the SFI-DfE Investigator Programme grant QuNaNet (grant 15/IA/2864), and the Leverhulme Trust project UltraQuTe. GTL and MP are funded by the FAPESP-QUB SPRINT program.

Appendix A The covariance matrix formalism

From $\bar{\boldsymbol{R}}=(\boldsymbol{a}_{1},\boldsymbol{a}_{1}^{\dagger},\cdots,\boldsymbol{a}_{L},\boldsymbol{a}_{L}^{\dagger})^{\mathsf{T}}$, we may define the $2L\times 2L$ covariance matrix (CM) $\mathsf{\mathsf{\Theta}}$ of our $L$-node network as

[TABLE]

Hence, in our convention, the CM of system with $L=2$ is

[TABLE]

This ordering turns out to be more convenient for the problem at hand. Indeed, $\mathsf{\Theta}$ can be readily decomposed in terms of the reduced covariance matrices

[TABLE]

One may then readily verify that

[TABLE]

where $\openone_{2L}$ is the $2L$-dimensional identity matrix and $\sigma_{\pm}$ are raising and lowering spin-$1/2$ matrices. In the NESS of the problem discussed in Sec. IV, $S$ is identically zero, so that one only needs the reduced matrix $\mathsf{C}$.

If the state of the system is Gaussian, the Wigner function (5) is completely determined by the CM and the vector of means $\bar{\mu}=\langle\bar{\boldsymbol{R}}\rangle$, as

[TABLE]

In passing, it is also convenient to note that the inverse CM $\mathsf{\Theta}^{-1}$ can be written as

[TABLE]

where \mathsf{B}=\big{[}\mathsf{C}-\mathsf{S}\,\big{(}\mathsf{C}^{-1}\big{)}^{\top}\mathsf{S}^{*}\big{]}^{-1} and $\mathsf{P}=-\mathsf{C}^{-1}\,\mathsf{S}\,\mathsf{B}^{\top}$. In particular, if $\mathsf{S}=0$, as in Sec. IV, one simply gets $\mathsf{B}=\mathsf{C}^{-1}$.

Appendix B Lyapunov equation

Since the dynamics is Gaussian preserving, the state can be characterized at all times using only the CM $\mathsf{\Theta}$ and the mean vector $\bar{\mu}$, when starting from a Gaussian initial condition. In particular, the dynamics of the CM can be cast as the Lyapunov equation

[TABLE]

where $\mathsf{W}$ and $\mathsf{F}$ are $2L\times 2L$ matrices, which follow from the underlying master equation. In the case of Eq. (3)

[TABLE]

where $\mathsf{\Gamma}=\text{diag}\,(\gamma_{1},\gamma_{1},\ldots,\gamma_{L},\gamma_{L})$ is a $2L\times 2L$ diagonal matrix with the specified diagonal elements and $\mathsf{H}$ is the matrix with elements $H_{k,\ell}$ which determine the structure of the network. In turn, the matrix $\mathsf{F}$ can be cast as $\mathsf{F}\coloneqq\mathsf{f}+\mathsf{\Gamma}/2$, where

[TABLE]

Since the master equation (3) does not spontaneously generate squeezing, it is possible to convert Eq. (40) into two separate equations for the reduced CMs $C$ and $S$ [cf. Eqs. (36)]. To accomplish this, one simply needs to exploit the tensor structure of Eqs. (41) and (42) together with that of Eq. (37). This gives

[TABLE]

Interestingly, we see that, while the equation for $\mathsf{C}$ acquires a clean structure, the equation for $\mathsf{S}$ becomes dependent on whether or not the matrix $\mathsf{H}$ is real or complex. In this respect, we note that the choice of phase in Eq. (27) is particularly convenient, as it entails $\mathsf{H}^{*}=-\mathsf{H}$ and Eq. (43b) acquires the same structure as Eq. (43a). We also remark that, while Eq. (43a) has an inhomogeneous term $\mathsf{f}$, this is not present in Eq. (43b). As a result, in the NESS we get $\mathsf{S}=0$ whereas $\mathsf{C}$ is the solution to

[TABLE]

Appendix C Entropy-production rate

The entropy production rate (17) can be written in terms of the entries of the CM. We begin by substituting explicitly the current (16) into Eq. (17), which yields

[TABLE]

Next, we substitute the explicit formula (38) for the Gaussian Wigner function in the logarithm and carry out the remaining integrals to obtain

[TABLE]

Although this is not as simple as the expression for the entropy flux in Eq. (20), this closed formula for $\Pi_{k}$ is equally useful.

Appendix D NESS from the model in Sec. IV

In this Appendix, we obtain the steady-state solution of Eq. (28) from Sec. IV.1. As discussed in Appendix B, this translates into solving an equation with the structure (44). In our specific case, the matrix $\mathsf{\Gamma}$ takes the form

[TABLE]

where the subindex “r” refers to quantities of the physical reservoirs. Similarly, the elements of the matrix $\mathsf{f}$ in Eq. (42) becomes

[TABLE]

where $\mathsf{C}_{\text{d}}$ is a matrix containing only the diagonal elements of $\mathsf{C}$, with all other entries being zero. We thus see that the self-consistent reservoirs introduce the elements of $\mathsf{C}$ directly into $\mathsf{f}$. Eq. (44) then becomes

[TABLE]

Interestingly, this is precisely the same equation that appears when one uses local dephasing to enforce diffusive heat conduction (cf., e.g., Ref. Žnidarič (2010); Manzano et al. (2012); Asadian et al. (2013)). We emphasize, however, that the NESS is different in both models since dephasing is not Gaussianity-preserving, i.e., they only coincide up to the second-order moments.

The solution of Eq. (49) was discussed in Žnidarič (2010); Manzano et al. (2012); Asadian et al. (2013) and yields a tridiagonal matrix of the form

[TABLE]

where $x\coloneqq\langle\boldsymbol{a}_{k}^{\dagger}\boldsymbol{a}_{k+1}\rangle$ is given in Eq. (29b) and the occupations are those of Eq. (29a).

Appendix E Calculation of $\Pi_{\text{r}}$ and $\Pi_{\text{sc}}$ for the model in Sec. IV

Finally, we discuss how the above tools can be used to facilitate the computation of the individual entropy production rates for the real and the self-consistent reservoirs; namely Eqs. (32) and (33). This is readily accomplished using Eq. (46). For the contribution to the real reservoir we get

[TABLE]

where the entropy flux $\Phi_{k}$ is given by Eq. (20). Similarly, the entropy production rate due to the self-consistent reservoirs evaluates to

[TABLE]

were we have used the fact that the occupations of the self-consistent reservoirs are $\tilde{n}_{k}=\langle\boldsymbol{a}_{k}^{\dagger}\boldsymbol{a}_{k}\rangle$, so that the corresponding entropy fluxes $\tilde{\Phi}_{k}$ are identically zero.