Nonequilibrium Thermodynamics and Steady State Density Matrix for Quantum Open Systems

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.


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 1,2 . More recent discussions on the origin of thermodynamical laws at the nanoscale can be found in, for example, 3 . Originally the weak coupling limit of a finite-size central region interacting with thermal and/or particle baths was first considered 1,2,4-6 . Methods for dealing with the strong coupling limit have been recently developed [7][8][9][10][11][12][13][14] . 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 10,11,13,15 . The long-time limit behaviour of the NE thermodynamics presents also very interesting properties [12][13][14] .
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 [16][17][18][19][20][21] 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 B 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 11 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 B. 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 22 or the NE specific heat of the central region.

A. 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 µ L and µ 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 [23][24][25][26][27][28][29][30][31] . 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 32,33 .
In the present paper, we consider that the full system is described by the Hamiltonian the three independent regions L, C, R. The interaction W is decomposed into several parts C where the interaction between particles in region C is given by V int C and the coupling between the C region and the α = L, R reservoirs is given by V αC + V Cα . Without specifying explicitly the form of H 0 , there exist different important commutation relations, i.e.
with α, β = L, C, R and N β is the occupation number operator of the different regions Initially, all regions L, C, R are isolated and characterised by their respective density matrix ρ α with α = L, C, R. The macroscopic L and R regions are represented by a density matrix ρ L,R expressed in the grand canonical ensemble, with temperature T α = 1/kβ α and chemical potential µ α (α = L, R): with Z α = Tr (α) [e −βα(Hα−µαNα) ] and Tr (α) implies a summation only over the states of the region α. The initial density matrix of the central region is assumed to take any arbitrary form ρ C as this region is not in the thermodynamic limit. Furthermore, considering ρ 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 ρ C from a microcanonical ensemble where with the eigenstates H C |C n = n |C n . The δ ∆ function is the "regularized" delta function defined by δ ∆ ( n −E C ) = 1 for E C ≤ n ≤ E C +∆ and 0 otherwise, and The total density matrix ρ 0 as the non-interacting state defined by H 0 is given by the B. The NE density matrix ρ NE In Ref. 34 , we used some concepts developed for asymptotic steady-state operators [35][36][37][38][39][40] and we have shown that the average of any arbitrary operator A in the NE asymptotic steady state is given by where ρ 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 ρ NE = Ω (+) ρ 0 Ω (+)−1 where the Moeller operator [41][42][43][44] , characterising the asymptotic steady state, is given by: Ω (+) = lim τ →−∞ e iHτ e −iH 0 τ . Such an operator presents a central property, the intertwining relation [40][41][42][43][44] , Ω (+) H 0 = HΩ (+) , or equivalently H By defining any asymptotic operator as X

it can be shown from
Eq. (1) that, when X α = H α or X α = N α , we have the following relation: Hence, any linear combination Y a = α a α X (+) α also commutes with H: [Y a , H] = 0.
The quantity Y a will be called a conserved quantity in the following. Furthermore, for β , it can be shown that [Y a , Y b ] = 0 when X α,β = H α,β or N α,β . This follows from Eq. (1) and consequently from [X We have now all the ingredients to study different expressions of the NE density matrix ρ NE in the steady state 35,40 . The latter can be recast as where, in the last equality, we used [X  First, we have generalised in 34 the results of Hershfield 38,39 to the presence of both a temperature and a chemical potential differences (β L = β R , µ L = µ R ) between the reservoirs.
The NE density matrix is then recast as follows: follows 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 The conserved quantities Y Q and Y E are related to the charge and energy currents respectively via: Second, following 40 , 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 (+) = Ω (+) EΩ (+)−1 with E = 1 2 (H L − H R ) and with Q = 1 2 (N L − N R ), the NE density matrix takes the following form (with X L+R = X L + X R ): Note that the passage from the first to the second line requires the use of an intertwining relation for N 46 .
Furthermore, the initial density matrix ρ 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 47 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 ρ C such that ρ C e +β(...) = 1 in Eq. (6) and Eq. (10) 48 .
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 and the charge current operator j Q (t) is given by Such a statistical operator is given by 18,26,49 : The quantity J S (s) is called the non-systematic energy flows 26 and is related to the entropy production rate of the system 25,26,49 . It is given by where all operators are given in the Heisenberg representation, A(u) = e iHu/ Ae −iHu/ . In the literature, it is also customary to call J S,α the heat current with results from the energy is the sum of heat flows divided by the subsystem temperature, it is the entropy production rate of the whole system 25,26 . The time integration of J S (u) in Eq. (11) provides the asymptotic steady state value of the energy and charge fluxes J E α and J Q α respectively. Hence the quantity duJ S (u) is the entropy production in the NE steady state.
Recently we have shown 34  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: where we have used the notation Y ≡ Y Q , Y E , Q (+) or E (+) and the interaction representation, A I (t) = e iH 0 t/ Ae −iH 0 t/ , for all quantities. The first values (n = 0) of the series are

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 tworeservoir device). We now use them to calculate the entropy production in the system under general NE conditions.

A. Entropy production rate
As mentioned in the previous section, the different quantities are related to the entropy production (rate) in the system. We can then define the NE entropy production ∆S NE in the steady state from Eq. (7) and Eq. (10) in the following way: where ∆Ṡ 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 j Q and energy j E : Such a result has also been used in early work 52 .
A few remarks are now in order. At equilibrium when T L = T R and µ L = µ 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 there is no energy flow and there is a charge current when µ L = µ R . By convention, the NE averaged charge current j Q is taken positive (negative) when flowing from L(R) to 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.

B. An example
In the absence of interaction, the Hamiltonian for the central region C is simply given by creates (annihilates) an electron in the level ε 0 . The non-interacting reservoirs are also described by a quadratic Hamiltonian where αi is an appropriate composite index to label the free electrons on the site i of the α reservoirs. The coupling between the central region and the electrodes is given via some hopping matrix elements v α , and we have The only non vanishing anti-commutators are {d, d † } = 1 and {c αi , c † βj } = δ ij δ αβ . The charge and energy currents can be calculated from the NE average expression in Eq.(4) 34 , from asymptotic steady state scattering techniques 12,53-57 or from a NE Green's functions (NEGF) approach 11,13 . The full equivalence between the asymptotic steady state scattering and the NEGF techniques has been shown in 58 .
In 49 , 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 59,60 . 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 22,[61][62][63][64][65][66][67] .
For the non-interacting system, the charge and energy currents are related to the transmission coefficient T (ω) of the junction via the moments M n : where f α (ω) is the equilibrium Fermi distribution of the reservoir α. The charge current is j Q = eM 0 and the energy current is j E = M 1 . The transmission is obtained from T (ω) = G r (ω)Γ L (ω)G a (ω)Γ R (ω) where the NEGF G r,a are given by G r (ω) = [ω − − Σ r L+R (ω)] −1 = (G a (ω)) * , with Σ r L+R = Σ r L + Σ r R being the reservoirs self-energy. Furthermore we have Γ α (ω) = Σ a α (ω) − Σ r α (ω) and the reservoir α self-energy is defined by Σ r α (ω) = v 2 α e −ikα(ω) /t α with the dispersion relation ω = ε α − 2t α cos(k α (ω)). Figure 1 shows the NE entropy production rate ∆Ṡ NE calculated for different transport regimes. The main conclusion is that ∆Ṡ 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 11 . The In Figure 2, we show how the NE entropy production rate ∆Ṡ NE depends on the NE conditions, i.e. on the chemical potential difference ∆µ (see left panels (a) and (c) in Fig.2) or on temperature differences ∆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 ∆Ṡ NE is always a positive quantity (as expected). Furthermore, ∆Ṡ NE increases when the NE conditions are more important, i.e. when ∆µ or ∆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 ∆Ṡ NE on ∆µ shows some form of linearity when ∆µ ≥ ∆µ * . This is simply due to the fact that the currents saturate: j Q,E (∆µ) = I sat is modified by an increase of ∆µ. In this regime, one can easily see that the dependence of ∆Ṡ NE on ∆µ is simply linear with a slope given by (β L + β R )I sat Q /2. Furthermore, the slope is maximal when β L = β R and smaller for any β L = β 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 ∆T differences (at fixed ∆µ) 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 68,69 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 68,69 S = k B ln W = −k B Tr[ρ ln ρ], which density matrix should be used? ). For the dependence on ∆µ, we take µ L = µ eq + ∆µ/2 and µ R = µ eq − ∆µ/2 with µ eq = 0.2. For the dependence on ∆T , we take and T R = T 0 = 0.1,

A. Which density matrix?
The first natural choice would be to take the NE density matrix ρ NE derived in the 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 ρ 0 , i.e.
As the density matrix ρ 0 is the direct product of the individual density matrices of each separate L, C, R regions, it is easy to show that S NE α is the contribution of the region α and ρ NE red,α is the corresponding reduced density matrix obtained from ρ NE red,α = Tr (β,β ) [ρ NE ] with β, β = L, C, R and β = β = α. For example, the NE reduced density matrix in the central region ρ NE red,C is obtained from ρ NE red,α = Tr We now further comment on this point. For that, we consider small deviations from the equilibrium where µ L,R = µ ± 1 2 ∆µ and β L,R = β ± 1 2 ∆β with ∆µ 1 (hence ∆ µ 1) and ∆β 1. Hence where we kept only the lowest order terms in ∆ µ and ∆β. Furthermore, if we assume a lowest order expansion of the density matrices e −β(H (+) and therefore The first term in the above equation is simply the entropy S 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: where ∆S NE is the operator defining the entropy production in Eq. (14), and S eq L = −k B ρ L ln ρ L with S eq L = Tr (L) [S eq L ]. Similar expressions can be found for S NE C and S NE R . 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 int C in the central region. This can however achieved by either determining the asymptotic scattering states |L (+) k = Ω (+) |L k for the L region (and for the states |R k and |C n for the R and C regions respectively). Following 12,53-57 , the scattering states of the L region, for the model described in Section III B, are given by: For the noninteracting case in Section III B, the calculations of the entropy can also be easily performed using the NEGF formalism which we consider in the next Section.

B. 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 B. We have shown that the NE steady state can be considered as a pseudo equilibrium state with a corresponding generalised Gibbs ensemble given by ρ NE . Following the same principles of equilibrium statistical mechanics, one can defined from the generalised Gibbs ensemble a local NE distribution functions 66 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 NE C in the central region C can be defined as follows 52 : where f NE C (ω) is the NE distribution function of the central region and A C (ω) is the corresponding spectral function defined from the NEGF as A C (ω) = −ImG r (ω)/π. It should be noted that the entropy S NE C is only a part of the total entropy S NE in Eq. (18), which is produced in the entire system under the general NE conditions. For non-interacting system considered in Sec. III B, the NE distribution function in the central region is just an weighted averaged of the equilibrium Fermi distributions of the Figure 3 shows the dependence of the entropy S NE C calculated for different transport regime. Once more, we can see that S NE C is always a positive quantity. The positiveness of S NE C 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 11 . 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 NE C has the dimension of an entropy, i.e.
In Figure 4, we show the dependence of the entropy S NE C on the NE conditions, i.e. on the chemical potential difference ∆µ, as shown in the left panels (a) and (c), and on temperature differences ∆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 NE C is once more a positive quantity (as expected). Furthermore, the entropy S NE Neumann NE entropy, one can also defined a NE distribution function f NE C which contains all the effects of the interactions as shown in 66 . 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 α, the latter are acting as equilibrium (particle and heat) baths with their respective temperature T α and chemical potential µ α . 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 ρ NE = e −β(H−μN + a=Q,E λaJa) /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 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 ρ 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 [10][11][12][13][14]  could note that another splitting of the total Hamiltonian H can be used to reduce the complexity of the expressions for the NE density matrix. In an earlier work 34 , we considered splitting the Hamiltonian H into H 0 + W where H 0 is only H 0 = H L + H R , hence the initial density matrix is only the direct product ρ 0 = ρ L ⊗ ρ R , and the expected form of the NE density matrix is obtained. A difference however occurs in the construction of the Y operators given in Eq. (13). In the present work, the operator W does not include H C while it does in 34 . In the calculation of the Gibbs-like entropy in the central region, one deals with products of terms including ρ(H + C ) ln ρ C . The asymptotic operator ρ(H + C ) can be expanded in a series of ρ n C from the series expansion of Moeller operators. Hence leading to a series of terms in ρ n C ln ρ C . By considering that initially the central region (of finite size) is fully isolate, there cannot be any partial occupation of the electronic levels, and hence the terms ρ n C ln ρ C expressed in the basis set of the central region will lead to the evaluation of either 1 n ln 1 or 0 n ln 0. This obviously leads to a zero contribution to the entropy, and therefore the terms in ρ(H + C ) in the density matrix can be ignored. 52 Similar results for the entropy production rate or for the Gibbs-von Neumann entropy have been also derived or used in Refs. [10][11][12][13][14] . A critical analysis of the results in 10,13 has been given in 71,72 . In 11 , no expression for the entropy production rate was given while its expression for the Gibbs -von Neumann entropy differs significantly from Eq. (23) due to the different time-dependent conditions. In the present paper, we do not consider that the central system is driven by an external time-dependent driving force. In 12 only the assumed standard definition for heat flux is used and hence Eq. (15) follows automatically.
Finally, only weak coupling regime was considered in 14 .