{'fraction_non_alphanumeric': 0.03790877796901893, 'fraction_numerical': 0.008885542168674699, 'mean_word_length': 4.829800576947197, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Numerical Approach to Non-Equilibrium Quantum Thermodynamics: Non-Perturbative Treatment of the Driven Resonant Level Model based on the Driven Liouville von-Neumann Formalism

Non-equilibrium thermodynamics of the driven resonant level model is studied using numerical simulations based on the driven Liouville von-Neumann formalism. The approach is first validated against recently obtained analytical results for quasi-static level shifts and the corresponding first order corrections. The numerical approach is then used to study far-from-equilibrium thermodynamic properties of the system under finite level shift rates. The proposed methodology allows the study of unexplored non-equilibrium thermodynamic regimes in open quantum systems.


Introduction
The quantum technology revolution promises unprecedented advances in our computation and technological capabilities (see Ref. [1] and references therein). As machines are scaled-down to the quantum regime, it is of prime importance to understand how quantum effects are manifested in their operation. Within this effort, an extension of the thermodynamic description to nanoscale systems out of equilibrium, where dynamic quantum effects dominate, constitutes an important challenge [2,3].
A generic problem of this kind is the thermodynamic analysis of a small quantum system that is coupled to multiple reservoirs, which are out of equilibrium with respect to each other, possibly with an additional external time-dependent force that performs work on the system [2][3][4][5][6][7][8][9]. As in standard macroscopic thermodynamics, such an analysis requires the partitioning of energetic variations in the system into heat and useful work components. However, the small system size makes the characterization of its thermodynamic properties uncertain because energy parameters associated with the system are of the same order as those characterizing the system-baths interactions. Furthermore, quantum mechanics implies that energy transfer and relaxation are associated with broadening of energy levels, so even energy observables of the system alone are not well characterized.
Models constructed to consider these issues are often studied in the weak system-bath coupling limit, where the thermodynamic functions of the system are well defined and it is possible to work with a standard quantum master equation (see Ref. [2] and references therein). Studies that consider the implications of strong system-bath coupling have recently emerged [4][5][6][10][11][12][13][14][15]. In the latter, the use of weak system-reservoir(s) coupling is replaced by the assumption that the dynamics imposed on the system is slow. "Slow" here implies that the timescale on which the system Hamiltonian is changed by external force(s) (as work is done or extracted from the system) is long relative to the characteristic timescales on which the system exchanges energy and particles with its environment. We also note that these treatments are usually limited to non-interacting systems (e.g., free-electron or Harmonic bath models), although extensions of such considerations to include electron-electron interactions in the bath have recently been published [15,16].
While such studies provide fundamental insights into dynamics and thermodynamics of small systems in these specific limits, it is obviously of interest to explore less restrictive conditions, e.g. strong coupling under arbitrary external driving. To this end, one may resort to numerical simulations [15,[17][18][19][20][21][22][23][24][25][26][27][28]. A convenient numerical framework, relevant to problems where a small electronic system is coupled to one or more free electron reservoirs, each in its own equilibrium but not necessarily in equilibrium with each other, is the Driven Liouville von Neumann (DLvN) method [29][30][31][32]35]. In this approach, the Liouville von-Neumann (LvN) equation of motion (EOM) for an extended system comprising the system of interest plus finite lead models is solved with open boundary conditions imposed on the leads. These boundary conditions are enforced by augmenting the LvN equations with non-unitary source/sink terms constructed to drive each lead towards its own equilibrium state. When the driving is done by imposing different equilibrium states on different leads, the DLvN method was shown to reproduce the results of explicit time-dependent non-equilibrium Green's function treatments [33], while avoiding violations of Pauli's exclusion principle and preserving density matrix positivitytwo issues that were encountered in previous related implementations [25,34]. This was formally supported by a derivation of the DLvN EOM as an approximation resulting from the non-equilibrium Green's function formalism [33,35] and its recasting in Lindblad form for non-interacting systems [33,37].
The purpose of the present work is to examine the applicability of the DLvN methodology for nonequilibrium thermodynamics simulations. For simplicity and clarity, we limit ourselves to the driven resonant level model considered in Refs. [4,5,7,10,13,38], represented by a time-dependent Hamiltonian comprising of a single electronic level (termed as "dot") coupled to a single equilibrium reservoir and driven by an external force that performs work on the system or extracts work from it. Some technical considerations relevant to the implementation of the DLvN method for evaluating thermodynamic properties of this model system were recently discussed [39]. Here, we use the method to calculate the dynamic evolution of this system numerically and investigate the thermodynamic implications of the driving without being limited to the weak coupling and/or slow driving limits considered in previous studies.

Methodology The Resonant Level Model System
As mentioned above, for our model system we choose the resonant level model previously considered in Refs. [4,5,7,10,13,38]. The model consists of a single spin-less electronic level that represents a quantum dot, coupled to a spin-less free-electron reservoir representing a metallic lead. The dot level is driven (e.g. by an external gate voltage or an electromagnetic field) such that its energy (and consequently its occupation) varies with time at a finite rate (see Figure 1). The Hamiltonian 4 describing this system is: where the dot and lead Hamiltonians, and , and the coupling between them are given by: and = ∑ ( † + ℎ. . ).
In Eqs.
(2)-(4) † ( ) and † ( ) are the creation (annihilation) operators for an electron in the dot and lead levels of energies and , respectively. For simplicity, the lead level energies and the dot/lead coupling matrix elements, , are kept real and constant during the dynamics. Furthermore, the lead is assumed to be at (or sufficiently close to) thermal equilibrium, characterized by electronic temperature and chemical potential . Without loss of generality, the latter is set to be at the energy origin such that = 0, as discussed below. Before considering far-from-equilibrium scenarios, we first demonstrate the performance of our numerical calculations at near-equilibrium conditions, where analytical results can be obtained from standard equilibrium quantum statistical mechanics [4]. As will become obvious below, the DLvN numerical scheme is not limited to some approximations invoked in analytical derivations, such as the wide band limit (WBL) and slow driving rates with respect to typical relaxation times in the system.
On the other hand, our numerical scheme can only address finite systems, which here implies a freeelectron lead modeled by a finite number of levels spanned within a finite bandwidth. Still, to facilitate comparison with the analytical results, we construct our model to be a good representation of the WB 5 approximation by representing the lead by a single energy band of width and a finite number, , of equally spaced (spacing of Δ = / and denisty of states = Δ −1 ) single electron levels and an energy independent lead-dot coupling, . The spectral function of the dot state is then given, to a good approximation, by a Lorentzian centered around : provided that ℏ ≫ Δ , ≫ ℏ and is far (relative to the Lorentzian width, ℏ ) from the band edges. For brevity of notation we use below ( ) to denote ( ; , ). When these conditions are not fully satisfied, proper corrections can be applied as described in Ref. [39]. It should be emphasized that deviations from the WB limit do not invalidate the thermodynamic calculations described below, only their comparison with the analytical WBL results.

The Driven Liouville-von-Neuman Approach
In the realm of the non-interacting spinless resonant level model described above, all observables can be obtained from the 1-electron density matrix ( ) = 〈 ( ) † ( )〉, where † ( ) and ( ) are the Heisenberg representations of the electron creation and annihilation operators of state i and the average is taken over the initial many-electron wavefunction. For the present treatment may be conveniently expressed in the basis of eigenstates of the isolated dot and lead sections: where is the dot population, is the density matrix block of the lead, and = † are the dot/lead coherences vectors. Correspondingly, the system Hamiltonian matrix representation is represented in the same basis as: where is a diagonal block of lead level energies, , that, as mentioned above, are taken to be uniformly distributed over the bandwidth, and is a column vector containing the identical coupling matrix elements that are related to a given dot-level broadening value of through Fermi's golden rule of Eq. (5).

6
As discussed above, a numerical evaluation of the dynamics described by the Hamiltonian (1)- (4) necessarily needs to employ a finite lead model, comprising a finite number of states spanned within a finite energy band, . The time evolution of the single-particle density matrix in such a closed system is given by the standard single-particle LvN EOM: = − ℏ [ , ]. In an out-of-equilibrium simulation, this dynamics results in Poincaré recurrences with a period determined by the energy spacing in the lead, Δ = / [16,17,22,26,40,41]. Therefore, the closed system model can mimic the behavior of its open counterpart only for short periods limited by the typical reflection time of the electronic wavepacket from the finite system boundaries [39,42]. To overcome this problem we adopt the recently developed DLvN approach for simulating time-dependent electronic transport in open quantum systems [26,[30][31][32][33][34]. Within this approach, open system boundary conditions are imposed by augmenting the LvN evolution with additional rate processes (see Eq. (9)) that guide the lead section towards its equilibrium form. The latter is represented by a diagonal matrix with populations given by the Fermi-Dirac distribution of the corresponding bath, to which the lead is implicitly coupled: where is the chemical potential of the implicit bath, = ( ) −1 , is Boltzmann's constant, and is the bath's electronic temperature. For brevity of notation we use below ( ) to denote ( ; , ).
The DLvN EOM for the resonant level model considered herein is of the following form: where The lead driving rate, Γ, represents the timescale on which thermal relaxation takes place in the lead due to its coupling to the implicit bath, which is generally assumed to be fast relative to all other processes of interest. While a physically motivated value for Γ can be extracted from the electronic properties of the explicit bath [34], in the present study we set ℏΓ, the lead levels broadening due to their coupling to the implicit bath, to be uniform and of the order of the lead level spacing, ℏΓ~Δ .
We have verified that the results presented below are fairly insensitive to the specific choice of lead driving rate in this range (see section 1 of the Supporting Information and relevant discussion in Refs. [37,43]), to the bandwidth of the finite lead model (see section 2 of the Supporting Information), and to the density of lead states (see section 3 of the Supporting Information).

Thermodynamic Functions and Fluxes
We now turn to describe our approach for evaluating thermodynamic functions and fluxes using the numerical scheme detailed above. We distinguish between quasi-static or sudden jump processes, where thermodynamic functions can be obtained from equilibrium calculations, and finite-rate processes, where kinetic equations are required.

The Quasi-Static Limit
In a quasi-static process, the driving parameter(s) (here ) is changed slowly (reversibly) relative to system relaxation processes so that the system can be assumed to remain at equilibrium, characterized by the equilibrium density matrix ( ) of the entire dot-lead system, throughout the process. The equilibrium value of a system observable that corresponds to an operator is given by: and its time-dependence (used in section 4 below for the evaluation of thermodynamic quantities in the quasi-static limit) is derived from the variation of , such that: The superscripts in Eqs. (11) and (12) denote the order in ̇ of the calculated quantity. Specifically, the WBL expressions for -dependent contributions to the particle number, energy, and entropy are given by (See Ref. [4]): (0), The corresponding quasistatic fluxes of particles, energy, work, heat, and entropy are given by ( [4]): In the numerical finite-bandwidth model calculation, the equilibrium density matrix, ( ) can be obtained using two different approaches. In one, for any instantaneous the Hamiltonian (7) is diagonalized and the equilibrium density matrix of the dot-lead system in the eigenstate representation, {| ⟩}, [ eq ( )] , ′ = , ′ ( ), is used to calculate the equilibrium expectation value of any singleelectron operator. For the thermodynamic functions of interest we get: To obtain -dependent expressions equivalent to Eqs. (13)-(15) these functions need to be projected onto the dot section as follows: The work done in the quasistatic process is obtained as an integral over the population: and the corresponding heat can be obtained from the first law: Alternatively, an approximate equilibrium density matrix may be calculated as the solution of a Sylvester equation [26,39] obtained by setting ̇= in Eq. (9) for any given value of . This can then be used to evaluate all needed expectation values of any single-electron observable 〈 〉 = Tr( ). The approximate nature of the latter stems from the fact that we impose thermal equilibrium only on the lead states and not on the eigenstates of the full dot-lead system. This would not matter if the lead was infinite, but small differences are expected when using finite lead models. Comparing the two solutions provides a measure of the suitability of the chosen finite lead model as an approximation for a wide-band lead [39] and can guide our effort to balance between the required accuracy and the computational cost. Within this approach, the occupation of the dot, which can be substituted (0) in expressions (26) and (27), is calculated using: In either case, the corresponding quasistatic fluxes are obtained as products of ̇ and the numerical derivatives with respect to d  of the above functions.

The Sudden Jump Limit
On the other extreme limit, the system starts at equilibrium with = 1 and at a given time, marked as = 0, the dot energy makes a sudden jump to = 2 without changing its population, after which the system is let to relax to the new equilibrium associated with 2 in place. Obviously, the work, which is performed only during the sudden jump, is given by: In the subsequent relaxation to equilibrium, the dot-lead system exchanges heat and particles with the external bath while no further work is done. The change in thermodynamic state functions is given by Δ = (0) ( 2 ) − (0) ( 1 ), ( = , , ) and the heat exchanged with the bath during the relaxation processes can be written in terms of the eigenstates and eigenenergies, (| ⟩, ) and (| ⟩, ), of the initial ( = 1 ) and final ( = 2 ) Hamiltonians of the dot-lead system, respectively. This leads to: ̇= Tr(̇) (32) ̇= Tr(̇( − )), where is a unit matrix of appropriate dimensions. Clearly, more heat is generated by the thermalization of electrons following the sudden jump than in the quasistatic process, hence > .

Finite Dot-Level Driving Rate
When the dot energy is shifted at a finite rate, the system goes out of equilibrium and irreversibility (1), (1), and These expressions are useful to describe processes involving vanishingly slow (with respect to the typical lead relaxation rate) dot level shifts. Since obtaining higher order analytical expressions rapidly become intractable, for any finite dot driving rate the particle and energy fluxes need to be evaluated from the kinetic equations. In what follows we discuss the numerical evaluation of these functions via the DLvN EOM (Eq. (9)).

Particle fluxes
Particle (electron) fluxes flowing between different system segments can be readily evaluated using the DLvN EOM. Given the 1-electron density matrix, , the instantaneous total number of electrons in the entire (dot-lead) system is given by ( ) = Tr[ ( )], which may vary with time due to the electron exchange with the implicit external bath. can be formally divided into contributions from the two system sections as follows: Since the dot is coupled only to the lead section, the particle flux flowing between the dot and the lead sections can be calculated using Eq. (9) as: where the sum runs over all lead state indices, . Similarly, for the lead section we obtain: As in Eq. (40), the first term on the right hand side of Eq. (41) is the particle flux flowing between the lead and the dot sections, with the opposite sign resulting from the directionality of the current. The second term can be identified as the outgoing particle flux from the lead to the implicit bath, to which it is coupled. Therefore, the particle influx into the bath is given by: Eqs. (40)- (42) can also be used to identify the energy resolved particle fluxes [38] that were used in Ref. [4] to evaluate heat currents flowing between the dot and the lead. However, due to technical considerations discussed in Ref. [38] we do not follow this route in the present paper.

Energy and energy fluxes
Next, consider the energy fluxes. We first consider the expectation value of the total electronic energy of the entire (dot-lead) system, and express it as a sum of the dot ( ) and lead ( ) components: where we have used the cyclic property of the full trace operator and defined Note that, since partial traces are used in the definitions of and , symmetrization is being employed to assure that the obtained values are real. We also note in passing that the contributions to these expressions from the coupling term of Eq. (4), , are equal (see Appendix B) such that This demonstrates that the contribution of the coupling to the total energy is split evenly between the dot and lead subsystems as was assumed in several recent papers [4,10].
The time variation of the dot's section contribution to the total electronic energy can now be evaluated as: The first term on the right-hand-side of Eq. (45) represents energy variations due to the flow of particles in-and out-of the dot, as obtained by the time-derivative of the single-particle density matrix.
The second term, which depends on the time-derivative of the Hamiltonian operator, corresponds in the present model to energy variations due to the external time-dependent perturbation acting on the dot section and is related to the work performed on the system, as discussed below. The former can thus be identified with the energy flux entering (or leaving if negative) the dot: where we have used Eq. (9) for ̇ and defined (see Eq. (4)): Similarly, the energy flux entering the lead section is given by: where we have used the fact that and (whose elements are given by Eq. (10)) are diagonal matrices. We may now define: where the second equality in Eq. (49) stems from the cyclic property the full trace operation: (49) can be obtained from the definition of (Eq. (7)) and (Eq. (6)) assuming that is a real coupling vector, and the last equality in the definition of results from the structure of and (Eq. (47)). With these definitions we can rewrite the energy fluxes as: and Since the lead section is coupled to both the dot (explicitly) and the bath (implicitly), energy conservation requires that the total energy flux into the lead equals the sum of energy fluxes out of the bath (− (E) ) and the dot sections: . Therefore, the energy flux into the implicit bath is given by: The structure of Eqs.

Heat fluxes
Turning to the consideration of heat currents, we encounter the usual conceptual problem associated with the fact that the distinction between energy and heat fluxes involves some form of relaxation. In the quasistatic limit discussed above, no ambiguity arises since this relaxation is assumed to be instantaneous on the driving timescale. When the driving is done at a finite rate, however, the definition of the heat flux becomes ambiguous unless some physical choices are made, as explained below.
Formally, the evolution of the system's (dot-lead) energy, = Tr( ( ) ( )), when the system undergoes a non-equilibrium process may be cast in the form: The first term on the right can be identified as the power (work per unit time) done on the system by the agent that makes the Hamiltonian time-dependent: In our model ̇( ) =̇( )| ⟩⟨ |, so that Energy conservation implies that the second term in (55)  For an open system (Г ≠ ), the second term on the RHS of (55) is determined by the processes that take place at the interface between the lead and the external bath, and expresses the effect of energy and particles exchange between the dot-lead system and this bath. In such a grand-canonical ensemble the total energy variation of the dot-lead system satisfies ̇=̇+̇+̇, where we recall that is the chemical potential of the bath. Using Eqs. (7), (9) and (55)-(56) we may write: or where is the unit matrix in the single electron d-L space and we have used again the fact that the trace of the commutator vanishes.
While formally exact, expressions (58)-(59) may miss the essence of the intended calculation. Since was introduced as a mathematical tool that depends on choices made for the finite lead size and the corresponding density of lead states, the heat flux calculated from equation (59) will depend on these unphysical choices. Specifically, we note that in the limit where the lead becomes infinite and its density of states becomes continuous, the dynamics of interest can be described by taking the limit → . In this limit, the calculated heat flux into the implicit bath vanishes (see Eq. (59)), while the lead itself is assumed to remain in thermal equilibrium. Indeed, this is the way the analytical calculation is done, by imposing a thermally equilibrated wide-band lead of uniform density-of-states and a coupling between the dot state and the eigenstates of the free lead, while disregarding any external bath that may be required to affect the lead thermalization. While such a picture of a constantly equilibrated lead may be valid when the dot level shift rate is not too fast, it cannot be implemented for the calculation of heat and particle exchange at finite shift rates within our finite-size model that addresses the lead dynamics explicitly.
An alternative approach for evaluating the heat flux can be devised based on the following two assumptions: (a) The process of interest, namely relaxation by transport of energy and particles between the dot and the lead, is irreversible; and (b) any change in the state of the lead resulting from this transport is insignificant in the sense that it does not affect the relaxation dynamics at the dot-lead interface. The former can be satisfied by taking > − , whereas the latter is obtained by choosing a sufficiently large lead model. Under these conditions, the actual relaxation of the lead to its thermal equilibrium by exchange of particles and energy with the external (implicit) bath can occur on a timescale much slower than the timescale of interest, namely the characteristic relaxation time at the dot-lead interface, without affecting the dynamics of interest. In these circumstances, all energy and particle fluxes at the dot-lead interface, being irreversible, are eventually expressed as heat and particle currents into the external bath. It is therefore convenient to define a book-keeping procedure associating ̇ and ̇ with the fluxes exchanged between dot and the lead (from lead to dot when positive), although their actual realization as fluxes in the external bath might take place on a different, 16 possibly much longer, timescale. Within this picture, the heat current (negative when out of the dot) is defined as: where is given by Eq. (50)  Eq. (60), on the other hand, is a meaningful representation of a quantity that would eventually be expressed as heat exchanged with the external thermal bath, however, it does not represent the actual heat exchanged in real-time. In this regard, it should be kept in mind that augmenting the finite dotlead system by coupling it to an external thermal bath was primarily done in order to impose irreversibility on the numerical procedure. It has no other dynamical consequences provided that the lead is taken large enough, namely fluxes at the dot-lead interface should be identical to what they would be if this lead was infinite once the calculation is converged with regards to lead size. Still, we find that the external bath is required for the conceptual definition of the heat flux, although, once the bookkeeping argument is adopted, the heat flux itself is calculated independently of this bath.
Another pertinent numerical point should be mentioned. In what follows we compare analytic results obtained at the WBL with numerical results calculated with a finite-band model. The latter are valid at any driving rate, ̇, whereas the former were obtained only for small ̇ with respect to the relaxation processes that thermalize the lead. While the numerical procedure described above is sufficiently accurate for many applications [29][30][31][32]35], the small shift rate corrections to the quasi-equilibrium behavior can be of the order of the small difference between the finite-and infinite-band results.
Furthermore, these differences between the finite-and infinite-band results converge slowly with the lead's bandwidth, in particular for the energy calculation (see Eq. (14)). We have discussed these numerical issues in details in Ref. [39], and have found that the following modification to the straightforward numerical calculation provides good agreement with the WBL analytical results. Let be, respectively, the differences between the equilibrium values and between the quasi-static change rates of a system property = , , calculated analytically in the wide-band limit from Eqs. (13) and (14) and numerically for the corresponding finite bandwidth model. Note that the latter cannot be simply represented by truncating the integrals in Eqs. (13) and (14), because the spectral function It should be noted that these flux expressions represent the difference between the quasistatic currents and the instantaneous currents obtained under a finite driving rate. However, care should be exercised when interpreting integrals of these corrections between two values of the externally driven parameter (here ). The reason is that these instantaneous currents do not contain information on the residual relaxation that takes place after the final value of the driving parameter is reached (note that in the sudden limit discussed above this relaxation dominates the non-equilibrium process). At the same time, it is obvious that the integral over the power ̇ accounts for the full work associated with the driven process.

Entropy and entropy fluxes
Finally, we consider entropy generation and entropy fluxes associated with the dot driving process.
Two points should be made at the outset: (i) in the quasistatic limit, the heat and entropy currents are related by ̇( 1) ( ) =̇( 1) ( )⁄ and are given in the WBL by Eq. (19); (ii) since the equilibrium entropy is a state-function, any driving protocol that carries the system between two equilibrium states corresponding to 1 and 2 induces the same entropy change ( 1 , 2 ) = (0) ( 1 , 2 )⁄ , determined by the corresponding quasistatic (reversible) process. Any excess heat generation resulting from the finite driving rate will be eventually expressed in the external bath, after the dot-lead system has reaches its final equilibrium. The total thermodynamic entropy production, expressed as entropy increase in the external bath which is taken to be at equilibrium throughout (i.e. the implicit bath of the numerical model), is: Here, ( 1, 2 |̇( )) is the heat produced under a given dot level shift protocol, ̇( ), that drives the dot between 1 and 2 , (e.g. ̇( ) = ( 2 − 1 ) ( ) for a sudden jump at = 0), and the negative sign indicates that this heat was transferred to the external bath. The first law implies that this excess heat is determined by the excess work namely (setting = 0): In departure from Eqs. (64) and (65) that explore thermodynamic changes associated with driving protocols that move a system between two equilibrium states, several recent studies, e.g. [4,7,10,44] have considered the instantaneous rates of variations in the particle number, work, heat, as well as entropy, to lowest non-trivial order in the dot driving speed. While such quantities can be calculated, (Eqs. (36)-(38)), their interpretation needs a careful examination as demonstrated by our discussion of the heat flux above. For the model under discussion there is no ambiguity concerning the power (work done or produced per unit-time) and the energy of the dot-lead system as well as its time-derivative are also well defined. With regards to the heat, we saw that we need to distinguished between the rate of heat escaping to the external (implicit) bath in real-time, Eq. (59) and that part of the energy current between the dot and the lead that will eventually be realized as heat in the external bath, Eq. (60).
We can use the bookkeeping argument discussed above for the entropy as well by identifying the excess entropy production as / , being the excess heat transferred to the external bath calculated from equation (60). Alternatively, we can follow Refs. [4,7,10,42] in calculating the effect of driving on the time-dependent entropy of the dot-lead system defined as: where is the many electron density matrix. The following interesting result was obtained for the lowest order excess quantity [4,7,10,44]: which serves to identify Ẇ (2) as the rate of entropy production. Indeed, Eq. (67) expresses the fact that this entropy production (positive if work is done on the system) is expressed as a sum of a term − Q̇ ( 2) (positive when heat exits the system) associated with the excess heat given to the external bath plus the excess change in the system (dot-lead) entropy, ̇( 2) . With this interpretation Q̇( 2) has to be the heat calculated from Eq. (59), which, as explained above, is model dependent and not a physically meaningful quantity. Furthermore, as was shown in Ref. [15], the result (67) cannot be extended to higher orders in the driving speed, so our finite rate calculation is not expected to rigorously show it.
In what follows, we examine our ability to numerically evaluate quantities analogous to those To calculate the corresponding entropy of the dot-lead system we evaluate Eq. (66) for the many-body density matrix , which needs to be expressed it in terms of the single-electron density matrix ( ) of Eq. (9). In Appendix A we show that the resulting expression for non-interacting systems is given by: where { ( )} is the set of eigenvalues of ( ). At equilibrium → ( ( ); , ), where { ( )} are the eigenvalues of the Hamiltonian of the dot-lead system at a given . In this limit Eq. (70) reduces to its equilibrium counterpart. Under a finite dot-level shift rate the excess entropy is given by: where ( ( )) is calculated via Eq. (70) by diagonalizing the dot-lead Hamiltonian for any given value of .
To conclude this section, we note that other entropy expressions, associated with the inner (dot) system only, can be considered. One can adopt for the dot entropy the definition of Esposito et al. [45]: where = Tr and B is the total system without the dot. However, properties derived for the evolution of this function assume an uncorrelated initial state, (0) = (0) , [45]. Alternatively, one can follow Ingarden [46], in focusing on the quantum mechanical operator entropy, , associated with a Hermitian operator ̂. quantifies the amount of missing information about the observable in a given system state . Following Ben Naim [47] (see also Lent [48]), for such an operator written in its eigenbasis representation, ̂= ∑ | ⟩ ⟨ |, the operator entropy is defined in terms of the probability = |⟨ | ⟩| 2 as = − ∑ 2 . For the operator ̂=̂ †̂, the corresponding operator entropy is the binary entropy associated with the dot, where = is the probability that the dot is occupied.

The Quasi-Static and Sudden Jump Limits
In this section, we present and discuss numerical results calculated for the quasistatic, sudden jump, and finite driving rate non-equilibrium processes. These results were obtained from the model of Eqs.
We start by considering the extreme limits of quasi-static driving and a sudden jump. As discussed above, the calculation of these equilibrium functions can be done by diagonalizing the Hamiltonian and occupying the eigen-levels according to the Fermi-Dirac distribution, or by calculating the 1electron equilibrium density matrix for the dot-lead system by solving Eq. (9) for ̇= 0, which may be recast in a Sylvester-type equation [25].    These results generally show good agreement between the different calculations, implying that our finite numerical model provides a reasonable representation of the wide-band limit, and that the Sylvester solution provides a decent approximation to the exact equilibrium density matrix. Note that small deviations do exist and are more appreciable in the heat and entropy differences calculations due to their small numerical values. These discrepancies are expected to become smaller as the size of the 24 numerical basis increases, i.e. when a larger lead model is used (see SI). We note that Figs. 2-4 show changes in thermodynamic functions assigned to the dot subsystem. We could alternatively show the changes in the corresponding quantities associated with the full dot-lead system. For non-interacting electrons and in the wide band limit the results for these changes should be identical. However, small differences are obtained for the present finite system model (see section 4 of the Supporting Information) because the system's spectral function slightly depends in this case on the dot level energy.

Finite Dot-Level Driving Rate
Next, consider shifting the dot's energy at a finite rate. Some preliminary notes are in place: (i) As mentioned above, recent analytical studies [4][5][6][7][10][11][12][13][14][15] have evaluated the lowest order (in the driving rate) corrections to thermodynamic functions and fluxes. To approach this limit in the numerical simulations we need to use relatively slow shift rates, for which deviations from the quasistatic limit are small, implying relatively large numerical errors. Such numerical errors become smaller at higher driving rates, where the numerical approach is obviously most useful. Still, to facilitate comparison with analytical results we have chosen to display, for any thermodynamic function , the normalized excess flux:

25
Focusing first on the dot's population, energy, and the work performed to move its level, we find that for slow dot-level shift velocities the corresponding excess thermodynamic fluxes show good agreement with the analytical results. Small differences can be attributed to contributions beyond first order and to the difference between the dynamics of our relatively small model system with a discrete lead level spectrum and the analytical resonant level model in the wide band limit. More important are the deviations associated with increasing dot level shift rates, that reflect the fact that relaxation to equilibrium lags behind the evolution of ( ). This demonstrates the capability of our numerical scheme to explore dynamical regimes that are difficult to access using analytical treatments.   For the heat flux, calculated using Eqs. (60) and (69) (Fig. 8a), we also find qualitative agreement with the analytical results at a slow driving rate and a similar lag behavior of the excess heat flux with increasing dot energy shift rate. Here, as well, the deviations between the numerical and analytical results can be associated with higher order contributions in the dynamical results and to the finite discrete nature of the numerical model. Another source of deviation is the above-mentioned approximate nature of the relaxation imposed on the system, driving the zero-order lead levels rather than the exact system eigenstates, to equilibrium. Naturally, the resulting differences between the numerical and analytical results are expressed more strongly in the heat calculation due to its overall smaller absolute magnitude obtained under the studied conditions. As discussed above, the results obtained by using Eq. (59) instead of (60) for the heat flux (Fig. 8b) represent the actual heat exchanged instantaneously between the dot-lead system and the implicit bath in the finite numerical model. These results do not correspond to the heat exchanged between the dot subsystem and the wide-band lead adopted in the analytical treatment and they depend on the choice of the driving rate Γ, which in turn is related to the discrete spectrum of the specific finite lead model [34].

Conclusions
The results presented above demonstrate the suitability of the Driven Liouville von-Neumann methodology for the study of non-equilibrium thermodynamic properties of open quantum systems, even in the regime of strong coupling between the subsystem of interest and its environment.
Specifically, we have focused on the resonant level model subjected to a time-dependent driving of the dot energy, but the same numerical approach can obviously be applied to a variety of other and more realistic models. Unlike recent analytical treatments of this problems, which rely on expansions in powers of the driving rate and are therefore limited to slow driving scenarios, our numerical approach applies also to systems subjected to high driving rates and, in fact, becomes more efficient and accurate in this regime. Furthermore, the numerical model is not limited to the wide band limit often used to simplify analytical treatments and can be used with arbitrary lead band structure. Here, however, we have deliberately attempted to stay close to the wide band and slow driving limits, at which the available analytical results are reliable, in order to facilitate comparison.
Our numerical results for the dot occupation, energy, and entropy as well as for the work and heat fluxes show excellent agreement with the analytical theory in the quasi-static limit. Furthermore, the DLvN predictions (with the exception of entropy) correspond well to the analytical results also for the case of finite slow dot level driving -the most challenging regime for our numerical simulations. We therefore conclude that the DLvN approach can be used as a complementary tool to analytical calculations, providing valuable information in dynamic and thermodynamic regimes that are difficult to explore by current analytical treatments.
One important aspect of applying thermodynamics to non-equilibrium situations clearly comes to light in the numerical analysis. In contrast to the energy, particle number, and work, the definitions of heat and entropy in non-equilibrium situations is subject to the usual uncertainty encountered whenever the bath or part of it (here the lead) is simulated explicitly as part of the system. In the present calculation, while the heat exchanged with the external bath and the entropy of the dot-lead system are readily evaluated, their evolution is determined by the DLvN driving rate (Γ) that is used as a tool to stabilize the numerical solution and depends on the choice of lead model size. Since this choice is somewhat arbitrary and depends on the convergence of the results, the physical meaning of the Γ dependent heat and entropy expressions become model dependent. More generally, the numerical calculation sheds light on the ambiguity associated with the fact that particle and energy fluxes exchanged by a system of interest and the realization of these fluxes as thermodynamic variations in an equilibrated external bath occur on different timescales. Alternative local estimates that focus on changes that take place at the dot-lead interface were explored. For example, the reasonable assumption that all the dot outgoing energy eventually translates into heat at the equilibrated bath was used to define, as a bookkeeping tool, an instantaneous heat flux out of the dot. For the entropy, local quantities such as the system entropy projected on the dot and the dot information entropy appear to be useful. It should be 30 emphasized that these aspects of our numerical results reflect not any difficulty in evaluating the time evolution of the driven system but only the ambiguity in interpreting these results in a thermodynamic framework.
Finally, the present numerical work has focused on a non-interacting particles (here fermions) model.
Generalization of the numerical procedure to better descriptions of the system, e.g. interacting electrons, will be explored in future studies.

Appendix A: Entropy in Terms of the Single-Particle Density Matrix.
In terms of the its many-electron density matrix operator, , the entropy of a system is given by: For a system of non-interacting electrons, Wick's theorem implies that can be written as an exponential of a free Fermion operator [42]: where is a scalar, = ( ) −1 , is Boltzmann's constant, is the temperature, and the manyelectron operator ̂ is given by: We start by finding an explicit expression for the pre-factor . To this end, we consider the unitary transformation matrix that diagonalizes . is constructed from the column eigenvectors of the matrix , which we denote by { }. Given we can transform the Fermion operators according to: where ̃ † denotes a row vector of creation operators, † , given in the general representation of the matrix , ̃ † denotes a row vector of creation operators, † , given in its eigenbasis representation and and are column vectors of the corresponding annihilation operators. The elements of the matrix itself can be written in terms of its eigenvalues, , and eigenvectors as follows: where is the ℎ element of the column vector .
We can write Eq. (A.2) in the diagonal representation as follows: where is the number of states and we have used the fact that all the operators in the exponents are commutable to write the corresponding Kronecker product.
The pre-factor can now be found from the condition Tr( ) = 1 by using the fact that the trace of a Kronecker product of two matrices is the multiplication of the traces of the individual matrices and that for spinless Fermionic systems Tr ( − † ) = 1 + − , where the two terms on the right-handside stand for the unoccupied and occupied ℎ eigenstate, respectively. The resulting expression is: . (A.8) Next, to find a relation between and we show that the two matrices share the same eigenbasis. The The expectation value of the total electronic energy of the entire (dot-lead) system is given by: We symmetrize this expression to ensure that all partial traces remain real-valued, and write is as a sum of the dot ( ) and lead ( ) components: where we have used the cyclic property of the full trace operator and defined: We can divide the Hamiltonian into its diagonal ( 0 ) and off-diagonal ( ) components: where the former represents the isolated dot and lead eigenstates, and the latter their mutual coupling, whose matrix elements are assumed to be real for simplicity.
We may now write as:

DLvN Driving Rate ( ) Sensitivity Check
The results presented in the main text were obtained using a finite lead model of

Finite Bandwidth Lead Model Sensitivity Check
The results presented in the main text were obtained using a finite lead model size of 100 levels spanning a bandwidth 10 ℏ . To verify that our results are sufficiently insensitive to this bandwidth choice we repeated the quasi-static dot occupation, entropy, and heat calculations using a lead Furthermore, the minor changes associated with doubling the bandwidth of the finite lead model make the agreement between the numerical calculation and the WBL expression even better. Therefore, we may conclude that our numerical results using the finite lead model are converged to the WBL case.
We note again that our numerical calculations are not limited to the wide-band approximation and the choice to approach this limit here is made deliberately for the purpose of comparison with the analytical expressions.

Lead Density of States Sensitivity Check
The results presented in the main text were obtained using a finite lead model of

Evaluation of the Lead Contribution to Quasi-Static Observable Variations
In the main text, when comparing the results of our numerical calculations to the WBL analytical expressions, we have focused on the dot's contribution to the variations in particle number and energy.
Indeed, in the analytical treatment the lead (representing the entire bath) is assumed to be constantly at equilibrium regardless of the dot dynamics. Hence, all variations occur at the dot itself and its interface with the lead. On the contrary, in our numerical treatment the coupled dot-lead dynamics is considered explicitly. Since the lead is driven towards equilibrium at a finite rate, its dynamical state in any given instance only approximates the desired equilibrium.
To test how much this influences the comparison with the WBL analytical results we have repeated the