Nonequilibrium Thermodynamics and Steady State Density Matrix for Quantum Open Systems

TL;DR

This paper investigates the nonequilibrium thermodynamics of a finite quantum electron system connected to reservoirs, deriving the steady state density matrix and analyzing entropy production, confirming its positivity and expanding theoretical understanding.

Contribution

It introduces a detailed analysis of the nonequilibrium density matrix for quantum open systems driven by reservoir differences, extending previous theoretical frameworks.

Findings

Entropy production rate is always positive.

Different expressions for entropy production are consistent.

Results expand and corroborate earlier theoretical work.

Abstract

We consider the generic model of a finite-size quantum electron system connected to two (temperature and particle) reservoirs. The quantum open system is driven out of equilibrium by the presence of both a temperature and a chemical potential differences between the two reservoirs. The nonequilibrium (NE) thermodynamical properties of such a quantum open system are studied for the steady state regime. In such a regime, the corresponding NE density matrix is built on the so-called generalised Gibbs ensembles. From different expressions of the NE density matrix, we can identify the terms related to the entropy production in the system. We show, for a simple model, that the entropy production rate is always a positive quantity. Alternative expressions for the entropy production are also obtained from the Gibbs-von Neumann conventional formula and discussed in detail. Our results corroborate and expand earlier works found in the literature.

Equations

$$\begin{array}{r}[H_{\alpha },H_{\beta }]=0,[H_{\alpha },N_{\beta }]=0,[N_{\alpha },N_{\beta }]=0\end{array}$$

$$\begin{array}{r}{\rho }_{\alpha }=\frac{1}{Z_{\alpha }}{e}^{-{\beta }_{\alpha }(H_{\alpha }-{\mu }_{\alpha }{N}_{\alpha })},\end{array}$$

$$\begin{array}{r}{\rho }_C={\rho }_C(H_C)=Z_C^{-1}\sum _n|C_n⟩{\delta }_{\Delta }({ϵ}_n-E_C)⟨C_n|\end{array}$$

$$\begin{array}{r}⟨A{⟩}^{\mathrm{NE}}=\mathrm{Tr}[{\rho }^{\mathrm{NE}}A],\end{array}$$

$$\begin{array}{r}[X_{\alpha }^{(+)},H]={\Omega }^{(+)}[X_{\alpha },H_0]{\Omega }^{(+)-1}=0.\end{array}$$

$$\begin{array}{rl}{\rho }^{\mathrm{NE}}& ={\Omega }^{(+)}{\rho }_0{\Omega }^{(+)-1}={\Omega }^{(+)}{\rho }_L{\Omega }^{(+)-1}{\Omega }^{(+)}{\rho }_C{\Omega }^{(+)-1}{\Omega }^{(+)}{\rho }_R{\Omega }^{(+)-1}\\ & =Z_L^{-1}{Z}_R^{-1}{e}^{-{\beta }_L(H_L^{(+)}-{\mu }_L{N}_L^{(+)})}{\rho }_C(H_C^{(+)})e^{-{\beta }_R(H_R^{(+)}-{\mu }_R{N}_R^{(+)})}\\ & =Z^{-1}\exp \left\{-{\beta }_L(H_L^{(+)}-{\mu }_L{N}_L^{(+)})-{\beta }_R(H_R^{(+)}-{\mu }_R{N}_R^{(+)})\right\}{\rho }_C(H_C^{(+)}),\end{array}$$

$$\begin{array}{r}{\rho }^{\mathrm{NE}}=Z^{-1}\exp \left\{-\bar{\beta }(H-Y^Q+Y^E)\right\}{\rho }_C(H_C^{(+)})e^{+\bar{\beta }{H}_C^{(+)}},\end{array}$$

$$\begin{array}{r}{Y}^Q=({\beta }_L{\mu }_L{N}_L^{(+)}+{\beta }_R{\mu }_R{N}_R^{(+)})/\bar{\beta },\mathrm{and}{Y}^E=({\beta }_L-{\beta }_R)\frac{1}{2}(H_L^{(+)}-H_R^{(+)})/\bar{\beta }.\end{array}$$

$$\begin{array}{rl}{Y}^Q& =\bar{\mu }N+{\Delta }_{\mu }{Q}^{(+)}/\bar{\beta }\end{array}$$

$$\begin{array}{rl}{\rho }^{\mathrm{NE}}& =Z^{-1}\exp \left\{-\bar{\beta }(H_{L+R}^{(+)}-\bar{\mu }{N}_{L+R}^{(+)})+{\Delta }_{\mu }{Q}^{(+)}-({\beta }_L-{\beta }_R)E^{(+)}\right\}{\rho }_C(H_C^{(+)})\\ & =Z^{-1}\exp \left\{-\bar{\beta }(H-\bar{\mu }N)+{\Delta }_{\mu }{Q}^{(+)}-({\beta }_L-{\beta }_R)E^{(+)}\right\}{\rho }_C(H_C^{(+)})e^{+\bar{\beta }(H_C^{(+)}-\bar{\mu }{N}_C^{(+)})},\end{array}$$

$$\begin{array}{r}{\rho }^{\mathrm{NESO}}=\frac{1}{Z}\exp \left\{-\sum _{\alpha }{\beta }_{\alpha }\left(H_{\alpha }-{\mu }_{\alpha }{N}_{\alpha }\right)+{\int }_{-\infty }^0\mathrm{d}u{e}^{\eta u}{J}_S(u)\right\}.\end{array}$$

$$\begin{array}{r}{J}_S(u)=\sum _{\alpha }{\beta }_{\alpha }{J}_{S,\alpha }(u)\mathrm{where}{J}_{S,\alpha }(u)=\frac{d}{du}(H_{\alpha }(u)-{\mu }_{\alpha }{N}_{\alpha }(u))\end{array}$$

$$\begin{array}{r}Y=\sum _{n=0}^{\infty }{Y}_{n,I}\mathrm{and}{\partial }_t{Y}_{n+1,I}(t)=-\frac{i}{\hslash }[W_I(t),Y_{n,I}(t)],\end{array}$$

$$\begin{array}{r}\Delta S^{\mathrm{NE}}/k_B={\Delta }_{\mu }⟨Q^{(+)}⟩-({\beta }_L-{\beta }_R)⟨E^{(+)}⟩=\bar{\beta }(⟨Y^Q⟩+⟨Y^E⟩-\bar{\mu }⟨N⟩)=\int \mathrm{d}\tau \Delta {\dot{S}}^{\mathrm{NE}}(\tau )/k_B,\end{array}$$

$$\Delta {\dot{S}}^{\mathrm{NE}}/k_B=({\beta }_L{\mu }_L-{\beta }_R{\mu }_R)⟨j_Q⟩/e-({\beta }_L-{\beta }_R)⟨j_E⟩.$$

$$\begin{array}{r}{M}_n=\frac{1}{\hslash }\int \frac{\mathrm{d}\omega }{2\pi }{\omega }^{n}T(\omega )(f_L(\omega )-f_R(\omega )),\end{array}$$

$$\begin{array}{r}\mathrm{Tr}[A^{(+)}f(A^{(+)})]=\mathrm{Tr}[A^{(+)}{\Omega }^{(+)}f(A){\Omega }^{(+)-1}]=\mathrm{Tr}[{\Omega }^{(+)}Af(A){\Omega }^{(+)-1}]=\mathrm{Tr}[Af(A)].\end{array}$$

$$S^{\mathrm{NE}}=-k_B\mathrm{Tr}[{\rho }^{\mathrm{NE}}\ln {\rho }_0]=S_L^{\mathrm{NE}}+S_C^{\mathrm{NE}}+S_R^{\mathrm{NE}}\mathrm{where}{S}_{\alpha }^{\mathrm{NE}}=-k_B{\mathrm{Tr}}_{(\alpha )}[{\rho }_{\mathrm{red},\alpha }^{\mathrm{NE}}\ln {\rho }_{\alpha }].$$

$$\begin{array}{rl}{\rho }^{\mathrm{NE}}& \approx Z^{-1}\exp \left\{-\beta (H_{L+R}^{(+)}-\mu N_{L+R}^{(+)})\right\}\left(1+{\Delta }_{\mu }{Q}^{(+)}\right)\left(1-\Delta \beta E^{(+)}\right){\rho }_C(H_C^{(+)})\\ & \approx Z^{-1}\exp \left\{-\beta (H_{L+R}^{(+)}-\mu N_{L+R}^{(+)})\right\}{\rho }_C(H_C^{(+)})\left(1+{\Delta }_{\mu }{Q}^{(+)}-\Delta \beta E^{(+)}\right),\end{array}$$

$${\rho }_{\mathrm{red},L}^{\mathrm{NE}}={\mathrm{Tr}}_{(C,R)}[{\rho }^{\mathrm{NE}}]\approx {\rho }_L{\mathrm{Tr}}_{(C,R)}\left[1+{\Delta }_{\mu }{Q}^{(+)}-\Delta \beta E^{(+)}\right],$$

$$\begin{array}{rl}{S}_L^{\mathrm{NE}}& =-k_B{\mathrm{Tr}}_{(L)}\left[{\rho }_{\mathrm{red},L}^{\mathrm{NE}}\ln {\rho }_L\right]\\ & \approx -k_B{\mathrm{Tr}}_{(L)}\left[{\rho }_L\ln {\rho }_L\right]-k_B{\mathrm{Tr}}_{(L)}\left[{\rho }_L{\mathrm{Tr}}_{(C,R)}\left[{\Delta }_{\mu }{Q}^{(+)}-\Delta \beta E^{(+)}\right]\ln {\rho }_L\right].\end{array}$$

$$S_L^{\mathrm{NE}}\approx S_L^{\mathrm{eq}}+\mathrm{Tr}\left[\frac{\Delta {\mathcal{S}}^{\mathrm{NE}}}{k_b}{\mathcal{S}}_L^{\mathrm{eq}}\right],$$

$$S_C^{\mathrm{NE}}(\Delta \mu ,\Delta T)=-k_B\int \frac{\mathrm{d}\omega }{2\pi }{A}_C(\omega )\left[f_C^{\mathrm{NE}}(\omega )\ln f_C^{\mathrm{NE}}+(1-f_C^{\mathrm{NE}}(\omega ))\ln (1-f_C^{\mathrm{NE}})\right],$$

Full text

I Introduction

The understanding of irreversible phenomena is a long-standing problem in statistical mechanics. Explanations of the fundamental laws of phenomenological nonequilibrium (NE) thermodynamics have been given and applied to quantum open systems for several decades Spohn (1978); Alicki (1979). More recent discussions on the origin of thermodynamical laws at the nanoscale can be found in, for example, Allahverdyan et al. (2004). Originally the weak coupling limit of a finite-size central region interacting with thermal and/or particle baths was first considered Spohn (1978); Alicki (1979); Davies (1978); Spohn (1978); Kosloff:2013 . Methods for dealing with the strong coupling limit have been recently developed  Campisi:2009 ; Deffner:2011 ; Ajisaka:2012 ; Ludovico:2014 ; Esposito:2015 ; Topp:2015 ; Bruch:2016 ; Solano:2016 . One can study the NE thermodynamical properties and the entropy production in such systems when an external driving force is applied to the central region. The time dependence of the external force can be arbitrary or periodic Ludovico:2014 ; Esposito:2015 ; Bruch:2016 ; Ludovico:2016 . The long-time limit behaviour of the NE thermodynamics presents also very interesting properties Topp:2015 ; Bruch:2016 ; Solano:2016 .

Indeed, after some time much longer than some typical relaxation times of the finite system, a steady state can be obtained. Such a state arises from the balance between irreversible processes (fluxes of particle and/or energy) and the driving forces induced by the reservoirs. The NE steady state presents some analogy to its equilibrium counterpart in the sense that an equilibrium state represents a stationary state of a closed system, while the NE steady state is the time-invariant state of an open system. The fact that the NE steady state can be seen as a pseudo-equilibrium state is central to the construction of the corresponding generalised Gibbs ensembles McLennan (1959); Zubarev (1974, 1994); Zubarev et al. (1996, 1997); Morozov and G.Röpke (1998) and to the calculation of the entropy, heat or work production under NE conditions.

In the present paper, we construct such a generalised Gibbs expression for the NE density matrix and apply it to the calculation of the entropy production in the system under the presence of both a temperature difference and a chemical potential difference between the two reservoirs. Our approach has no restriction for the nature of the coupling (strong or weak) to the reservoirs, nor for the presence (or absence) of interaction between particle in the central region.

In Section II, we provide different, but fully compatible, expressions for the NE density matrix and show that the new terms in the NE density matrix (new in the sense that they do not appear in the equilibrium grand-canonical density matrix) are associated with the entropy production under the NE conditions. The entropy production rate is shown to be related to the fluxes of particle and heat across the system (Section III). We provide in Section III.2 some numerical calculations of the entropy production rate for a model system of a single electron resonance coupled to two Fermi reservoirs. A comparison to earlier results Esposito:2015 is also given. In Section IV, we consider the NE entropy production in the entire system obtained from the Gibbs - von Neumann expression based on the NE density matrix. Explicit derivations are provided in some limiting cases and it is shown that NE entropy is produced not only in the central region but also in the reservoirs. For the single resonance model, we also calculate the NE Gibbs - von Neumann entropy in the central region and present the corresponding results in Section IV.2. Our approach corroborate and extend earlier existing results. Furthermore it opens a new route to the calculations of the full NE response functions of the system, such as the NE charge susceptibility Ness and Dash (2012) or the NE specific heat of the central region.

II Non-equilibrium steady state

II.1 System and initial conditions

We consider a finite-size central region $C$, connected two electrodes (left $L$ and right $R$) acting as thermal and particle reservoirs. These electrodes are described within the thermodynamics limit. Initially they are at their own equilibrium, characterized by two temperatures $T_{L}$ and $T_{R}$, and by two chemical potentials $\mu_{L}$ and $\mu_{R}$. Furthermore, we ignore the interaction between particles in the electrodes, although the central region $C$ may contain such kind of interaction. We are interested in steady state regime, and therefore we take the initial state of the system to be in the far remote past. The system is then characterised by an Hamiltonian $H_{0}$. After all parts of the system are “connected” together and after some time elapses, the full system is considered to reach a NE steady state. The system is then characterised by the total Hamiltonian $H=H_{0}+W$.

The questions related to the possibility of reaching a NE steady-state have been addressed in Ruelle (2000); Tasaki and Matsui (2003); Fröhlich et al. (2003); Tasaki and Takahashi (2006); Maes and Netočný (2010); Tasaki et al. (2011); Moldoveanu et al. (2011); Cornean:2011 ; Cornean:2014 . It is also been argued that a system will always reach a steady-state if it is a (or if it is connected to another) system in the thermodynamic limit regardless the presence (or absence) of adiabatic switching of the interactions Ojima (1989); Cornean et al. (2008).

In the present paper, we consider that the full system is described by the Hamiltonian $H=H_{0}+W$ where $H_{0}$ is the non-interacting Hamiltonian $H_{0}=H_{L}+H_{C}+H_{R}$ built from the three independent regions $L,C,R$. The interaction $W$ is decomposed into several parts $W=V_{LC}+V_{CL}+V_{RC}+V_{CR}+V_{C}^{\rm int}$ where the interaction between particles in region $C$ is given by $V_{C}^{\rm int}$ and the coupling between the $C$ region and the $\alpha=L,R$ reservoirs is given by $V_{\alpha C}+V_{C\alpha}$. Without specifying explicitly the form of $H_{0}$, there exist different important commutation relations, i.e.

[TABLE]

with $\alpha,\beta=L,C,R$ and $N_{\beta}$ is the occupation number operator of the different regions $\alpha,\beta=L,C,R$.

Initially, all regions $L,C,R$ are isolated and characterised by their respective density matrix $\rho_{\alpha}$ with $\alpha=L,C,R$. The macroscopic $L$ and $R$ regions are represented by a density matrix $\rho_{L,R}$ expressed in the grand canonical ensemble, with temperature $T_{\alpha}=1/k\beta_{\alpha}$ and chemical potential $\mu_{\alpha}$ ($\alpha=L,R$):

[TABLE]

with $Z_{\alpha}={\rm Tr}_{(\alpha)}[e^{-\beta_{\alpha}(H_{\alpha}-\mu_{\alpha}N_{\alpha})}]$ and ${\rm Tr}_{(\alpha)}$ implies a summation only over the states of the region $\alpha$. The initial density matrix of the central region is assumed to take any arbitrary form $\rho_{C}$ as this region is not in the thermodynamic limit. Furthermore, considering $\rho_{C}$ to be given by a canonical or a grand canonical ensemble would imply the presence of the third reservoir, which is not ideal in the present case. Therefore, we define $\rho_{C}$ from a microcanonical ensemble where

[TABLE]

with the eigenstates $H_{C}|C_{n}\rangle=\epsilon_{n}|C_{n}\rangle$. The $\delta_{\Delta}$ function is the “regularized” delta function defined by $\delta_{\Delta}(\epsilon_{n}-E_{C})=1$ for $E_{C}\leq\epsilon_{n}\leq E_{C}+\Delta$ and 0 otherwise, and $Z_{C}={\rm Tr}_{(C)}[\delta_{\Delta}(H_{C}-E_{C})]$.

The total density matrix $\rho_{0}$ as the non-interacting state defined by $H_{0}$ is given by the direct product $\rho_{0}=\rho_{L}\otimes\rho_{C}\otimes\rho_{R}$.

II.2 The NE density matrix $\rho^{\rm NE}$

In Ref. Ness (2014), we used some concepts developed for asymptotic steady-state operators Doyon and Andrei (2006); Fujii (2007); Gelin and Kosov (2009); Hershfield (1993); Hyldgaard (2012); Bernard and Doyon (2013) and we have shown that the average of any arbitrary operator $A$ in the NE asymptotic steady state is given by

[TABLE]

where $\rho^{\rm NE}$ is the NE steady state density matrix. One should also note that the trace in Eq. (4) runs over all the states of the three $L,C,R$ regions. The NE density matrix is defined from $\rho^{\rm NE}=\Omega^{(+)}\rho_{0}\Omega^{(+)-1}$ where the Moeller operator Gell-Mann and Goldberger (1953); Akhiezer and Peletminskii (1981); Bohm (1993); Baute et al. (2001), characterising the asymptotic steady state, is given by: $\Omega^{(+)}={\rm lim}_{\tau\rightarrow-\infty}\ e^{iH\tau}e^{-iH_{0}\tau}$. Such an operator presents a central property, the intertwining relation Gell-Mann and Goldberger (1953); Akhiezer and Peletminskii (1981); Bohm (1993); Baute et al. (2001); Bernard and Doyon (2013), $\Omega^{(+)}H_{0}=H\Omega^{(+)}$, or equivalently $H_{0}^{(+)}=\Omega^{(+)}H_{0}\Omega^{(+)-1}=H$.

By defining any asymptotic operator as $X_{\alpha}^{(+)}=\Omega^{(+)}X_{\alpha}\Omega^{(+)-1}$ , it can be shown from Eq. (1) that, when $X_{\alpha}=H_{\alpha}$ or $X_{\alpha}=N_{\alpha}$, we have the following relation:

[TABLE]

Hence, any linear combination $Y^{a}=\sum_{\alpha}a_{\alpha}X_{\alpha}^{(+)}$ also commutes with $H$: $[Y^{a},H]=0$. The quantity $Y^{a}$ will be called a conserved quantity in the following. Furthermore, for $Y^{b}=\sum_{\beta}b_{\beta}X_{\beta}^{(+)}$, it can be shown that $[Y^{a},Y^{b}]=0$ when $X_{\alpha,\beta}=H_{\alpha,\beta}$ or $N_{\alpha,\beta}$. This follows from Eq. (1) and consequently from $[X_{\alpha}^{(+)},X_{\beta}^{(+)}]=0$.

We have now all the ingredients to study different expressions of the NE density matrix $\rho^{\rm NE}$ in the steady state Fujii (2007); Bernard and Doyon (2013). The latter can be recast as

[TABLE]

where, in the last equality, we used $[X_{\alpha}^{(+)},X_{\beta}^{(+)}]=0$ and $Z=Z_{L}Z_{R}$. Finally, one should also note that $H_{L}^{(+)}+H_{R}^{(+)}=\Omega^{(+)}(H_{L}+H_{R})\Omega^{(+)-1}=\Omega^{(+)}(H_{0}-H_{C})\Omega^{(+)-1}=H-H_{C}^{(+)}$.

II.3 Three equivalent expressions for $\rho^{\rm NE}$

Upon regrouping the different terms in the exponential of Eq. (6), one obtains different, but equivalent, expression for the density matrix. The first expression is a generalisation of the density matrix derived by Hershfield Hershfield (1993); Hyldgaard (2012), the third expression is the so-called McLennan-Zubarev NE statistical operator, while the second expression is an intermediate between the two.

First, we have generalised in Ness (2014) the results of Hershfield Hershfield (1993); Hyldgaard (2012) to the presence of both a temperature and a chemical potential differences ($\beta_{L}\neq\beta_{R}$, $\mu_{L}\neq\mu_{R}$) between the reservoirs. The NE density matrix is then recast as follows: follows

[TABLE]

where we have used the definitions and commutators given in Sec. II. Note that the generalised Gibbs form of the NE density matrix in Eq. (7) is given with an effective temperature $T_{\rm eff}$ defined from $\bar{\beta}=\frac{1}{2}(\beta_{L}+\beta_{R})$. This temperature is different from the temperature of the reservoirs $T_{L,R}$ since $T_{\rm eff}=1/k_{B}\bar{\beta}=2T_{L}T_{R}/(T_{L}+T_{R})$.

The conserved quantities $Y^{Q}$ and $Y^{E}$ are related to the charge and energy currents respectively via:

[TABLE]

Second, following Bernard and Doyon (2013), one can re-expressed the density matrix in a slightly different form (closer to the grand-canonical ensemble) involving quantities with a more explicit physical meaning. Indeed, by writing $E^{(+)}=\Omega^{(+)}E\Omega^{(+)-1}$ with $E=\frac{1}{2}(H_{L}-H_{R})$ and

[TABLE]

with $Q=\frac{1}{2}(N_{L}-N_{R})$, the NE density matrix takes the following form (with $X_{L+R}=X_{L}+X_{R}$):

[TABLE]

where $\bar{\mu}=(\beta_{L}\mu_{L}+\beta_{R}\mu_{R})/(\beta_{L}+\beta_{R})$, $\Delta_{\mu}=\beta_{L}\mu_{L}-\beta_{R}\mu_{R}$ Note:ablinearcombination . Note that the passage from the first to the second line requires the use of an intertwining relation for $N$ note:N+ .

Furthermore, the initial density matrix $\rho_{C}$ could be given by any other form different from Eq. (3) as such a choice is completely arbitrary. Indeed, in the steady state, the initial correlation vanishes Velicky:2010 and the final stationary properties should not dependent on the initial conditions taken for the statistics of the central region. Hence, for convenience and to simplify the notation, we chose $\rho_{C}$ such that $\rho_{C}\ e^{+\bar{\beta}(...)}=\mathbbm{1}$ in Eq. (6) and Eq. (10) Note:comment_on_rhoC .

It is also important to note that, in Eq. (10), the quantities $Q^{(+)}$ and $E^{(+)}$ are conserved quantities and are directly related to the charge and energy currents. Indeed, in the Heisenberg representation, the energy current operator is given by $j_{E}(t)=\partial_{t}E(t)=\frac{i}{\hbar}[H,E(t)]$ and the charge current operator $j_{Q}(t)$ is given by $j_{Q}(t)=e\partial_{t}Q(t)=\frac{ie}{\hbar}[H,Q(t)]$. This results permits us to connect the expressions Eq. (6) and Eq. (10) for the NE density matrix to the third formulation, i.e. the McLennan - Zubarev NE statistical operator.

Such a statistical operator is given by Zubarev (1994); Tasaki and Takahashi (2006); Ness (2013):

[TABLE]

The quantity $J_{S}(s)$ is called the non-systematic energy flows Tasaki and Takahashi (2006) and is related to the entropy production rate of the system Fröhlich et al. (2003); Tasaki and Takahashi (2006); Ness (2013). It is given by

[TABLE]

where all operators are given in the Heisenberg representation, $A(u)=e^{iHu/\hbar}Ae^{-iHu/\hbar}$. In the literature, it is also customary to call $J_{S,\alpha}$ the heat current with results from the energy flux $J^{E}_{\alpha}(t)=\frac{d}{dt}H_{\alpha}(t)$ measured with respect to the so-called convective term $\mu_{\alpha}J^{Q}_{\alpha}(t)=\mu_{\alpha}\frac{d}{dt}N_{\alpha}(s)$ Sierra:2015 . Because $J_{S}(u)$ is the sum of heat flows divided by the subsystem temperature, it is the entropy production rate of the whole system Fröhlich et al. (2003); Tasaki and Takahashi (2006). The time integration of $J_{S}(u)$ in Eq. (11) provides the asymptotic steady state value of the energy and charge fluxes $J^{E}_{\alpha}$ and $J^{Q}_{\alpha}$ respectively. Hence the quantity $\int{\rm d}uJ_{S}(u)$ is the entropy production in the NE steady state.

Recently we have shown Ness (2014) the full equivalence between the McLennan - Zubarev NE statistical operator $\rho^{\rm NESO}$ and the other expressions Eq. (7) and Eq. (10) for $\rho^{\rm NE}$. The equivalence is based on the so-called Peletminskii lemma Peletminskii:1972 which states that the time integral of an operator given in the Heisenberg representation (for example $\int{\rm d}uJ_{S}(u)$) can be obtained from an infinite series expansion of the time integral of the related quantities expressed in the interaction representation (see Appendix B in Ness (2014)).

Hence the such an equivalence implies that the quantities $Y^{Q,E}$ in Eq. (7) and the quantities $Q^{(+)},E^{(+)}$ in Eq. (10) can be calculated from the same formal iterative scheme:

[TABLE]

where we have used the notation $Y\equiv Y^{Q},Y^{E},Q^{(+)}$ or $E^{(+)}$ and the interaction representation, $A_{I}(t)=e^{iH_{0}t/\hbar}Ae^{-iH_{0}t/\hbar}$, for all quantities. The first values ($n=0$) of the series are $Y_{0,I}=a^{Q}_{L}N_{L}+a^{Q}_{R}N_{R}$ when $Y\equiv Y^{Q}$ or $Q^{(+)}$ and $Y_{0,I}=a^{E}(H_{L}-H_{R})$ when $Y\equiv Y^{E}$ or $E^{(+)}$. The different constants $a^{Q}_{\alpha}$ and $a^{E}$ are given by $a^{Q}_{\alpha}={2\beta_{\alpha}\mu_{\alpha}}/(\beta_{L}+\beta_{R})$ for $Y=Y^{Q}$. For $Y=Q^{(+)}$ we have $a^{Q}_{L}=\frac{1}{2}=-a^{Q}_{R}$. For the energy flux, we have $a^{E}={\beta_{L}-\beta_{R}}/(\beta_{L}+\beta_{R})$ for $Y=Y^{E}$ or $a^{E}=\frac{1}{2}$ for $Y=E^{(+)}$.

III Entropy production

The equations Eq. (7) and Eq. (10) correspond to the most general expressions of the steady-state NE density matrix in the presence of both heat and charge currents (for a two-reservoir device). We now use them to calculate the entropy production in the system under general NE conditions.

III.1 Entropy production rate

As mentioned in the previous section, the different quantities $Y=Y^{Q},Y^{E},Q^{(+)}$ or $E^{(+)}$ are related to the entropy production (rate) in the system. We can then define the NE entropy production $\Delta S^{\rm NE}$ in the steady state from Eq. (7) and Eq. (10) in the following way:

[TABLE]

where $\Delta\dot{S}^{\rm NE}$ is the NE entropy production rate. Hence, from the definition of $Q^{(+)}$ and $E^{(+)}$, the NE entropy production rate is directly related to the asymptotic steady state NE current of charge $\langle j_{Q}\rangle$ and energy $\langle j_{E}\rangle$:

[TABLE]

Such a result has also been used in early work Note:comments_on_entropy .

A few remarks are now in order. At equilibrium when $T_{L}=T_{R}$ and $\mu_{L}=\mu_{R}$, there are obviously no current flows and no extra entropy is produced (apart from the equilibrium entropy arising from the thermal fluctuations in the $L$ and $R$ reservoirs). When $T_{L}=T_{R}$, there is no energy flow and there is a charge current when $\mu_{L}\neq\mu_{R}$. By convention, the NE averaged charge current $\langle j_{Q}\rangle$ is taken positive (negative) when flowing from $L(R)$ to $R(L)$, i.e. when $\mu_{L}>\mu_{R}$ ($\mu_{L}<\mu_{R}$). Hence the contribution $(\mu_{L}-\mu_{R})\langle j_{Q}\rangle$ to the entropy production is always positive. Similarly, the NE averaged energy current $\langle j_{E}\rangle$ is taken positive (negative) when flowing from $L(R)$ to $R(L)$, i.e. when $T_{L}>T_{R}$ ($T_{L}<T_{R}$). Hence the contribution $-(\beta_{L}-\beta_{R})\langle j_{E}\rangle$ to the entropy production is also always positive. However, the sign of the total contribution from $\langle j_{Q}\rangle$ and $\langle j_{E}\rangle$ to the entropy production (under general NE conditions) is not obvious without further investigation. A general argument for the positiveness of the NE entropy production rate was given in Topp:2015 and a few numerical examples were given in Esposito:2015 .

Using the same model system, i.e. the non-interacting single level coupled to two reservoirs, we provide in the next Section results for the entropy production rate for a wide range of parameters.

III.2 An example

In the absence of interaction, the Hamiltonian for the central region $C$ is simply given by $H^{0}_{C}=\varepsilon_{0}d^{\dagger}d$ where $d^{\dagger}$ ($d$) creates (annihilates) an electron in the level $\varepsilon_{0}$. The non-interacting reservoirs are also described by a quadratic Hamiltonian $\alpha=L,R$ with $H_{\alpha}=\sum_{i}\varepsilon_{\alpha}c^{\dagger}_{\alpha i}c_{\alpha i}+t_{\alpha}(c^{\dagger}_{\alpha i}c_{\alpha i+1}+c.c.)$ where ${\alpha i}$ is an appropriate composite index to label the free electrons on the site $i$ of the $\alpha$ reservoirs. The coupling between the central region and the electrodes is given via some hopping matrix elements $v_{\alpha}$, and we have $\sum_{\alpha}(V_{C\alpha}+V_{\alpha C})=\sum_{\alpha}v_{\alpha}\left(c^{\dagger}_{\alpha 0}d+d^{\dagger}c_{\alpha 0}\right)$. We recall that, by definition, we have $W=H^{0}_{C}+\sum_{\alpha}(V_{C\alpha}+V_{\alpha C})$. The only non vanishing anti-commutators are $\{d,d^{\dagger}\}=1$ and $\{c_{\alpha i},c^{\dagger}_{\beta j}\}=\delta_{ij}\delta_{\alpha\beta}$.

The charge and energy currents can be calculated from the NE average expression in Eq.(4) Ness (2014), from asymptotic steady state scattering techniques Han (2006, 2007); Han and Heary (2007); Han (2010a, b); Topp:2015 or from a NE Green’s functions (NEGF) approach Esposito:2015 ; Bruch:2016 . The full equivalence between the asymptotic steady state scattering and the NEGF techniques has been shown in Han et al. (2012).

In Ness (2013), we have stressed that calculating the NE averages with the NE density matrix and the series expansion of the operators $Y$ in Eq. (13), is equivalent to the NEGF approach in the steady-state regime. The Green’s functions are correlation functions whose thermodynamical averages are formally identical to those given in Eq.(4). Both perturbation series used in the NEGF approach and in the derivations of the equations for the $Y$ operators start from the same nonequilibrium series expansion. They are just two different ways of summing that series. For a non-interacting problem for which the series can be resumed exactly, the NEGF and the NE density matrix with the $Y$ operators approach provide the same result Schiller:1995 ; Schiller:1998 . For an interacting system, one must resort to approximations to partially resum the series, and therefore the two approaches are similar only when the same approximations are used. For the purpose of the present section, we then use the NEGF approach as the calculations are more straight forward in the non-interacting case. We also note that the NEGF formalism permits us to include local interaction in the central region in a compact and self-consistent scheme, as we have done in Dash et al. (2010); Ness:2010 ; Dash et al. (2011); Ness and Dash (2011); Dash et al. (2012); Ness and Dash (2012); Ness (2014); Ness and Dash (2014).

For the non-interacting system, the charge and energy currents are related to the transmission coefficient $T(\omega)$ of the junction via the moments $M_{n}$:

[TABLE]

where $f_{\alpha}(\omega)$ is the equilibrium Fermi distribution of the reservoir $\alpha$. The charge current is $\langle j_{Q}\rangle=eM_{0}$ and the energy current is $\langle j_{E}\rangle=M_{1}$. The transmission is obtained from $T(\omega)=G^{r}(\omega)\Gamma_{L}(\omega)G^{a}(\omega)\Gamma_{R}(\omega)$ where the NEGF $G^{r,a}$ are given by $G^{r}(\omega)=[\omega-\epsilon-\Sigma^{r}_{L+R}(\omega)]^{-1}=(G^{a}(\omega))^{*}$, with $\Sigma^{r}_{L+R}=\Sigma^{r}_{L}+\Sigma^{r}_{R}$ being the reservoirs self-energy. Furthermore we have $\Gamma_{\alpha}(\omega)=\Sigma^{a}_{\alpha}(\omega)-\Sigma^{r}_{\alpha}(\omega)$ and the reservoir $\alpha$ self-energy is defined by $\Sigma^{r}_{\alpha}(\omega)=v_{\alpha}^{2}e^{-ik_{\alpha}(\omega)}/t_{\alpha}$ with the dispersion relation $\omega=\varepsilon_{\alpha}-2t_{\alpha}\cos(k_{\alpha}(\omega))$.

Figure 1 shows the NE entropy production rate $\Delta\dot{S}^{\rm NE}$ calculated for different transport regimes. The main conclusion is that $\Delta\dot{S}^{\rm NE}$ is always a positive quantity, as expected. Such a behaviour is obtained for a system with a single chemical potential (see panel (a) in Fig. 1). It is also obtained when there are both a chemical potential and temperature differences between the reservoirs, regardless of the respective direction of the charge and energy currents, see panel (b) for currents flowing in the same direction and panel (b) for currents flowing in opposite directions.

In panel (d) Fig. 1, we have tried to reproduce the results shown in Fig. 3(b) of Esposito:2015 . The results are qualitatively reproduced apart from the behaviour of the over whole amplitude of the entropy production rate. Indeed, in our model, the transmission coefficient $T(\omega)$ has roughly a Lorentzian lineshape, with a maximum amplitude of unity (whatever the values of the parameters are) and a width which scales approximately as $\sum_{\alpha}v_{\alpha}^{2}/t_{\alpha}$ versus the coupling parameters between the central region $C$ and the reservoirs. Hence the width of $T(\omega)$ increases with the strength of the coupling to the reservoirs, and therefore the currents will always have a larger values when increasing the strength of this coupling. Consequently, the entropy production rate $\Delta\dot{S}^{\rm NE}$ defined by Eq.(15) is always larger for larger values of the coupling to the reservoirs.

In Figure 2, we show how the NE entropy production rate $\Delta\dot{S}^{\rm NE}$ depends on the NE conditions, i.e. on the chemical potential difference $\Delta\mu$ (see left panels (a) and (c) in Fig.2) or on temperature differences $\Delta T$ between the reservoirs (see right panels (b) and (d) in Fig.2). First it is important to note that, for all the parameters used, the NE entropy production rate $\Delta\dot{S}^{\rm NE}$ is always a positive quantity (as expected). Furthermore, $\Delta\dot{S}^{\rm NE}$ increases when the NE conditions are more important, i.e. when $\Delta\mu$ or $\Delta T$ increases. In other words, the more the system is out of equilibrium, the larger the entropy production becomes.

One should note that the dependence of $\Delta\dot{S}^{\rm NE}$ on $\Delta\mu$ shows some form of linearity when $\Delta\mu\geq\Delta\mu^{*}$. This is simply due to the fact that the currents saturate: $\langle j_{Q,E}(\Delta\mu)\rangle=I^{\rm sat}_{Q,E}$ when $\Delta\mu\geq\Delta\mu^{*}$. Indeed, in the saturation region, increasing $\Delta\mu$ does not change the value of the moments $M_{n}$ as the transmission $T(\omega)$ is zero in energy range where $f_{L}(\omega)-f_{R}(\omega)$ is modified by an increase of $\Delta\mu$. In this regime, one can easily see that the dependence of $\Delta\dot{S}^{\rm NE}$ on $\Delta\mu$ is simply linear with a slope given by $(\beta_{L}+\beta_{R})I^{\rm sat}_{Q}/2$. Furthermore, the slope is maximal when $\beta_{L}=\beta_{R}$ and smaller for any $\beta_{L}\neq\beta_{R}$ as clearly exemplified by the results shown in panels (a) and (c) in Figure 2. Such a saturation regime does not exist for increase $\Delta T$ differences (at fixed $\Delta\mu$) as shown in the panels (b) and (d) in Figure 2.

IV Nonequilibrium Gibbs - von Neumann entropies

In the previous Section, we have shown how the NE entropy production rate is related to the charge and energy currents. We have also shown that the NE steady state can be considered as a pseudo equilibrium state with a corresponding (time-independent) density matrix which is given in the form of a generalised Gibbs ensemble. It is therefore be very interesting to be able to define a NE entropy Penrose:1979 ; Maroney:2007 from the NE density matrix by using the equivalence between pseudo equilibrium states and equilibrium states. In other words: when we build a NE entropy from the equilibrium expression Penrose:1979 ; Maroney:2007 $S=k_{B}\langle\ln\mathcal{W}\rangle=-k_{B}{\rm Tr}[\rho\ln\rho]$, which density matrix should be used?

IV.1 Which density matrix?

The first natural choice would be to take the NE density matrix $\rho^{\rm NE}$ derived in the previous Section. However such a choice does not bring any information about entropy production under the NE conditions. Indeed, if we consider the asymptotic operator $A^{(+)}$ being obtained from a unitary transformation ($\Omega^{(+)}\Omega^{(+)-1}=\mathbbm{1}$), we can show that for any function $f(A)$:

[TABLE]

By taking $A=\rho_{0}$ and $f(.)=\ln(.)$, one easily see that ${\rm Tr}[\rho^{\rm NE}\ln\rho^{\rm NE}]={\rm Tr}[\rho_{0}\ln\rho_{0}]$. The quantity ${\rm Tr}[\rho_{0}\ln\rho_{0}]=\sum_{\alpha={L,C,R}}{\rm Tr}_{(\alpha)}[\rho_{\alpha}\ln\rho_{\alpha}]$ defines the entropy of the three separated $L,C,R$ regions. It does not contains any information about the charge and energy currents following through the entire system under general NE conditions.

Another possibility would be to take the NE average of the density matrix of the coupled system at equilibrium, i.e. $S^{\rm NE}=-k_{B}{\rm Tr}[\rho^{\rm NE}\ln\rho^{\rm eq}]$ where $\rho^{\rm eq}=Z^{-1}{\rm exp}\left\{-\beta^{\rm eq}(H-\mu^{\rm eq}N)\right\}$. However, from the intertwining relation, we have $\Omega^{(+)-1}H\Omega^{(+)}=H_{0}$ and $\Omega^{(+)-1}N\Omega^{(+)}=N$ note:N+ , and we obtain $S^{\rm NE}=k_{B}{\rm Tr}[\rho^{\rm NE}\ln\rho^{\rm eq}]=-k_{B}{\rm Tr}[\rho_{0}\ln\rho_{0}^{\rm eq}]$ with $\rho^{\rm eq}\propto{\rm exp}\left\{-\beta^{\rm eq}(H_{0}-\mu^{\rm eq}N)\right\}$. Such an entropy contains some information about the NE conditions, considering that $\beta_{L,R}$ and $\mu_{L,R}$ in $\rho_{L,R}$ are different from the equilibrium $\beta^{\rm eq}$ and $\mu^{\rm eq}$. However this entropy is defined from the non-interacting Hamiltonian only, and it lacks the presence of the operator $W$ which is the generator of the different charge and energy currents. Hence such an entropy does not contain any information about the fluxes which are the responses to the applied forces $\Delta\mu$ and $\Delta T$.

One has to go back to the definition of the NE steady state averages given in Eq. (4). The asymptotic time-dependence, in such average, has been passed on to the NE density matrix which we use to calculate the average of quantum operators. Hence, it follows that one should define the entropy from the NE average of the nominal density matrix $\rho_{0}$, i.e. $S^{\rm NE}=-k_{B}{\rm Tr}[\rho^{\rm NE}\ln\rho_{0}]$.

As the density matrix $\rho_{0}$ is the direct product of the individual density matrices of each separate $L,C,R$ regions, it is easy to show that

[TABLE]

$S_{\alpha}^{\rm NE}$ is the contribution of the region $\alpha$ and $\rho^{\rm NE}_{{\rm red},\alpha}$ is the corresponding reduced density matrix obtained from $\rho^{\rm NE}_{{\rm red},\alpha}={\rm Tr}_{(\beta,\beta^{\prime})}[\rho^{\rm NE}]$ with $\beta,\beta^{\prime}=L,C,R$ and $\beta\neq\beta^{\prime}\neq\alpha$. For example, the NE reduced density matrix in the central region $\rho^{\rm NE}_{{\rm red},C}$ is obtained from $\rho^{\rm NE}_{{\rm red},\alpha}={\rm Tr}_{(L,R)}[\rho^{\rm NE}]$.

The corresponding entropy $S_{C}^{\rm NE}$ has been the object of recent studies Ludovico:2014 ; Esposito:2015 ; Topp:2015 ; Bruch:2016 ; Solano:2016 but it is clearly only a part of the entire entropy production in the system. For example, the contributions $S_{L}^{\rm NE}$ and $S_{R}^{\rm NE}$ are different from their (isolated) equilibrium counter parts $S_{\alpha}^{\rm eq}=-k_{B}{\rm Tr}_{(\alpha)}[\rho_{\alpha}\ln\rho_{\alpha}]$ ($\alpha=L,R$) since $\rho^{\rm NE}_{{\rm red},\alpha}\neq\rho_{\alpha}$.

We now further comment on this point. For that, we consider small deviations from the equilibrium where $\mu_{L,R}=\mu\pm\frac{1}{2}\Delta\mu$ and $\beta_{L,R}=\beta\pm\frac{1}{2}\Delta\beta$ with $\Delta\mu\ll 1$ (hence $\Delta_{\mu}\ll 1$) and $\Delta\beta\ll 1$. Hence

[TABLE]

where we kept only the lowest order terms in $\Delta_{\mu}$ and $\Delta\beta$. Furthermore, if we assume a lowest order expansion of the density matrices $e^{-\beta(H_{\alpha}^{(+)}-\mu N_{\alpha}^{(+)})}\approx e^{-\beta(H_{\alpha}-\mu N_{\alpha})}$ and $\rho_{C}(H_{C}^{(+)})\approx\rho_{C}(H_{C})$, one gets:

[TABLE]

and therefore

[TABLE]

The first term in the above equation is simply the entropy $S^{\rm eq}_{L}$ of the isolated $L$ region with the associated grand canonical density matrix given by Eq. (2). The second term can be re-arranged as follows: $-k_{B}{\rm Tr}_{(L)}[\rho_{L}{\rm Tr}_{(C,R)}[\Delta_{\mu}Q^{(+)}-\Delta\beta E^{(+)}]\ln\rho_{L}]=-k_{B}{\rm Tr}[\rho_{L}(\Delta_{\mu}Q^{(+)}-\Delta\beta E^{(+)})\ln\rho_{L}]={\rm Tr}[(\Delta_{\mu}Q^{(+)}-\Delta\beta E^{(+)})(-k_{B}\rho_{L}\ln\rho_{L})]$. Finally, we have

[TABLE]

where $\Delta\mathcal{S}^{\rm NE}$ is the operator defining the entropy production in Eq. (14), and $\mathcal{S}^{\rm eq}_{L}=-k_{B}\rho_{L}\ln\rho_{L}$ with $S^{\rm eq}_{L}={\rm Tr}_{(L)}[\mathcal{S}^{\rm eq}_{L}]$. Similar expressions can be found for $S_{C}^{\rm NE}$ and $S_{R}^{\rm NE}$.

The result show that, under general NE conditions, NE entropy is produced in the central region and in the reservoirs as well. Such an entropy is always related to the charge and energy currents flowing at the interfaces between the central region and the $L$ and $R$ regions.

The full calculation of the entropy from Eq. (18) is a non trivial task, especially for arbitrary interaction $V_{C}^{\rm int}$ in the central region. This can however achieved by either determining the asymptotic scattering states $|L_{k}^{(+)}\rangle=\Omega^{(+)}|L_{k}\rangle$ for the $L$ region (and for the states $|R_{k}\rangle$ and $|C_{n}\rangle$ for the $R$ and $C$ regions respectively). Following Han (2006, 2007); Han and Heary (2007); Han (2010a, b); Topp:2015 , the scattering states of the $L$ region, for the model described in Section III.2, are given by: $|L_{k}^{(+)}\rangle=|L_{k}\rangle+v_{L}G^{r}(\epsilon_{Lk})|C_{n}\rangle+\sum_{\alpha=L,R;k^{\prime}}v_{L}v_{\alpha}G^{r}(\epsilon_{Lk})/(\epsilon_{Lk}-\epsilon_{\alpha k^{\prime}}+i0^{+})|\alpha_{k^{\prime}}\rangle$. For the non-interacting case in Section III.2, the calculations of the entropy can also be easily performed using the NEGF formalism which we consider in the next Section.

IV.2 An example for the entropy of the central region

We now consider numerical calculations for the Gibbs - von Neumann entropy using the single level model described in Section III.2. We have shown that the NE steady state can be considered as a pseudo equilibrium state with a corresponding generalised Gibbs ensemble given by $\rho^{\rm NE}$. Following the same principles of equilibrium statistical mechanics, one can defined from the generalised Gibbs ensemble a local NE distribution functions Ness (2014) in the $L,C,R$ regions. From these NE distributions functions, one can also defined the corresponding Gibbs - von Neumann entropies. For example, the NE entropy $S_{C}^{\rm NE}$ in the central region $C$ can be defined as follows Note:comments_on_entropy :

[TABLE]

where $f_{C}^{\rm NE}(\omega)$ is the NE distribution function of the central region and $A_{C}(\omega)$ is the corresponding spectral function defined from the NEGF as $A_{C}(\omega)=-{\rm Im}G^{r}(\omega)/\pi$. It should be noted that the entropy $S_{C}^{\rm NE}$ is only a part of the total entropy $S^{\rm NE}$ in Eq.(18), which is produced in the entire system under the general NE conditions.

For non-interacting system considered in Sec. III.2, the NE distribution function in the central region is just an weighted averaged of the equilibrium Fermi distributions of the reservoirs $f_{C}^{\rm NE}(\omega)=[\Gamma_{L}(\omega)f_{L}(\omega)+\Gamma_{R}(\omega)f_{R}(\omega)]/\Gamma_{L+R}(\omega)$.

Figure 3 shows the dependence of the entropy $S_{C}^{\rm NE}$ calculated for different transport regime. Once more, we can see that $S_{C}^{\rm NE}$ is always a positive quantity. The positiveness of $S_{C}^{\rm NE}$ is obtained when the system has a single chemical potential (see panel (a) in Fig. 3) as well as when there are both a chemical potential and temperature differences between the reservoirs (see panels (b) and (c) in Fig. 3). In panel Fig. 3.(d), we show the behaviour of the entropy for the same parameter used in Fig. 1.(d) and we recover the same qualitative behaviour as shown in Fig. 3(b) of Esposito:2015 . The amplitude of the entropy is larger in the weak coupling limit in comparison to the strong coupling limit to the reservoirs. One should however note that in Fig. 3, $S_{C}^{\rm NE}$ has the dimension of an entropy, i.e. [energy]/[temperature], while in Fig. 3(b) of Esposito:2015 and Fig. 1, we are dealing with an entropy production rate, i.e. a quantity with dimension [energy]/[temperature $\times$ time].

In Figure 4, we show the dependence of the entropy $S_{C}^{\rm NE}$ on the NE conditions, i.e. on the chemical potential difference $\Delta\mu$, as shown in the left panels (a) and (c), and on temperature differences $\Delta T$ between the reservoirs, as shown in the right panels (b) and (d). One can see that, for the range of parameters we used, the NE entropy production $S_{C}^{\rm NE}$ is once more a positive quantity (as expected). Furthermore, the entropy $S_{C}^{\rm NE}$ increases with the NE conditions, i.e. it increases for increasing values of $\Delta\mu$ and/or $\Delta T$. We also observe a saturation regime in $S_{C}^{\rm NE}$ with increasing $\Delta\mu$. In the saturation regime, an increase of $\Delta\mu$ changes the features of the NE distribution function $f_{C}^{\rm NE}(\omega)$ in an energy range where the spectral function $A_{C}(\omega)$ has no weight, i.e. where $A_{C}(\omega)=0$. Therefore the energy integral in Eq. (23) does not change with increasing $\Delta\mu$ and the entropy $S_{C}^{\rm NE}$ saturates.

Finally, one should note that the calculation of the Gibbs - von Neumann NE entropy for the central region as well as the calculation of the entropy production rate Eq. (15) can also be performed when interactions are present in the central region. Our expressions for the entropies are generally applicable to the cases with and without interactions in the central region. The latter is directly related to the charge and energy currents which can be calculated for different kinds of interaction in the central region. For example, in Dash et al. (2010); Ness:2010 ; Dash et al. (2011, 2012); Ness and Dash (2012), we have studied the effect of electron-vibration interaction on the electron current. For Gibbs - von Neumann NE entropy, one can also defined a NE distribution function $f_{C}^{\rm NE}$ which contains all the effects of the interactions as shown in Ness (2014). The interactions will affect the entropy production, however we expect that, with the so-called conservative approximations for the interaction, the positiveness of the entropy will be conserved. In the presence of interaction with extra degrees of the freedom (vibration or other boson modes), the contribution of their respective entropy production will need to be taken into accout. However such an in-depth study is out of the scope of the present paper.

V Discussion and conclusion

We have studied the steady state NE thermodynamical properties of an open quantum system connected to two reservoirs $\alpha$, the latter are acting as equilibrium (particle and heat) baths with their respective temperature $T_{\alpha}$ and chemical potential $\mu_{\alpha}$. We have shown that the steady state of the entire system can be seen as a pseudo equilibrium state. The corresponding NE density matrix is expressed in the form of a generalised Gibbs ensemble $\rho^{\rm NE}=e^{-\bar{\beta}(H-\bar{\mu}N+\sum_{a=Q,E}\lambda_{a}\mathcal{J}_{a})}/Z$.

The NE density matrix is time independent and built from the so-called conserved quantities: the total Hamiltonian $H$ and the total number of electrons $N$ and the $\mathcal{J}_{Q,E}$ quantities which are related to the fluxes of charge and energy flowing in between the central region $C$ and the reservoirs. We have given different forms for the NE density matrix and shown their mutual equivalence. The extra terms entering the definition of $\rho^{\rm NE}$ which do not exist in the equilibrium grand canonical representation have been clearly identified and have been shown to be related to the entropy production in the entire system. From their expression, the entropy production rate is given in terms of the charge and energy currents.

We have calculated such an entropy production rate for a model system consisting in a single electron resonance coupled to two Fermi reservoirs. Numerical results performed for different transport regimes have shown that the entropy production rate is always a positive quantity.

Furthermore, based upon the pseudo equilibrium properties of the steady state, we have also calculated a Gibbs - von Neumann entropy for the entire system. Our results show that the NE conditions create extra entropy in the central region as well as in the reservoirs. The former can be derived from the equilibrium expression of the entropy by using the appropriate NE distribution function in the central region.

Our numerical results for the entropy production and production rate corroborate and expand earlier studies Ludovico:2014 ; Esposito:2015 ; Topp:2015 ; Bruch:2016 ; Solano:2016 . These results also open a new route for determining the NE thermodynamical properties of quantum open systems under general conditions. For example the corresponding NE specific heat or charge susceptibility Ness and Dash (2012) can be directly obtained from the derivative of the entropy versus the applied temperature or chemical potential biases.

Acknowledgements.

The author thanks Benjamin Doyon and Lev Kantorovich for fruitful discussions. The UK EPSRC is acknowledged for financial support under Grant No. EP/J019259/1.