DFT-inspired methods for quantum thermodynamics

In the framework of quantum thermodynamics, we propose a method to quantitatively describe thermodynamic quantities for out-of-equilibrium interacting many-body systems. The method is articulated in various approximation protocols which allow to achieve increasing levels of accuracy, it is relatively simple to implement even for medium and large number of interactive particles, and uses tools and concepts from density functional theory. We test the method on the driven Hubbard dimer at half filling, and compare exact and approximate results. We show that the proposed method reproduces the average quantum work to high accuracy: for a very large region of parameter space (which cuts across all dynamical regimes) estimates are within 10% of the exact results.

a crucial quantity in DFT. Despite, being a very popular and successful approach to quantum many-body physics, to our knowledge, DFT or related methods have not been applied, so far, to quantum thermodynamics. Here we will consider tools from a particular flavour of DFT, lattice-DFT 50 , but our method could be straightforwardly extended to different DFT flavours. Lattice-DFT applies to model Hamiltonians such as the Hubbard or Heisenberg Hamiltonians, which are particularly important for quantum information processing and that describe well the type of spin systems of interest to the quantum thermodynamic community. As a test-bed example we will study the out-of-equilibrium dynamics of the Hubbard dimer, which provides a variety of interesting regimes and behaviours to benchmark our method, since it can be solved exactly. We will compare exact and approximate results and show that our method provides very accurate estimates in most regimes of interest opening attractive further possibilities for applications.

The Hubbard Model
The Hubbard model was conceived in 1963 separately by Gutzwiller 51 , Kanamori 52 , and Hubbard 53 to describe the interaction of electrons in solids. This model gives a microscopic understanding of the transition between the Mott-insulator and conductor systems 54 , and allows for tunnelling of particles between adjacent sites of a lattice and interactions between particles occupying the same site. Here we consider the time-dependent one-dimensional Hubbard Hamiltonian 55  , , J is the hopping parameter, Δ i (t) is the time-dependent amplitude of the spin-independent external potential at site i, and U is the on-site particle-particle interaction parameter. Analytical or numerically exact solutions for static many-body problems are rare. For the Hubbard model, there exist an exact solution to the homogeneous one-dimensional case, however numerically exact solutions to the non-uniform case becomes quickly problematic as the number of particles and sites increases. Solutions to time-dependent many-body problems are even more difficult to achieve, and the Hubbard model is no exception. However quantum work due to an external driving of a many-body dynamics is one of those problems in which including time dependence and the effects of many-body interactions is essential to capture -even qualitatively -the system behaviour. In the following, we propose a relatively simple approach to this problem which allows to include the key features stemming from the many-body interactions and dynamics, and yet maintaining the simplicity of simulating a non-interacting dynamics. We will show that, for a wide range of parameters, this simple approach allows to reproduce quantitatively the exact results.

DFT-inspired estimate of quantum thermodynamic quantities
The Kohn-Sham system. The KS system is a fictitious, non-interacting, quantum system defined as having the same ground state particle density as the original interacting physical system 49 . The two systems have then the same number of particles. For each physical system, and a given many-body interaction, the KS system is uniquely defined and, in the limit in which the physical system is non interacting, the Hohenberg-Kohn theorem 48 ensures that the KS and the physical system coincide. Exceptions are generally pathological cases: for Lattice-DFT, see 56 , where the extension to lattice-DFT of the one-to-one correspondence between external potential and ground-state particle density is demonstrated and related exceptions are discussed. Given an interacting system of Hamiltonian = + +ˆˆĤ K V V ee , where K is the kinetic energy operator, V is the one-body external potential, and V ee is the electron-electron repulsion, the corresponding KS system is described by the non-interacting Hamiltonian Here V H is the Hartree potential corresponding to the classical electrostatic interaction and V xc is the exchange-correlation potential, the functional derivative of the exchange-correlation energy E xc , which contains additional contributions from the many-body interactions as well as the many-body contributions to the kinetic energy. Usually, the exchange-correlation potential is defined as a sum of exchange and correlation contributions, c , and, due to its unknown functional dependence on the ground-state particle density, approximations have to be used for calculating it.
Work in a quantum system. In a non-equilibrium driven quantum system work is defined as the mean value of a work probability distribution ∫ = W WP W dW ( ) that takes into account energy-level transitions (system histories) that can occur in a quantum dynamics 57 . In this scenario an external agent is performing (extracting) work on (from) a quantum system and the concept of work takes into account both the intrinsic nondeterministic nature of quantum mechanics and the effects of non-equilibrium fluctuations. The work probability distribution P(W) contains all the information about the possible transitions in the Hamiltonian energy spectrum produced by an external potential (field) that drives the system out-of-equilibrium between the initial time t = 0 and the final time t = τ of a protocol. This distribution is defined as n m n m n m n , , , where p n is the probability to find the system in the eigenstate |n〉 (with eigenenergy ε n 0 ) of the initial Hamiltonian = H t ( 0), p m|n is the transition probability to evolve the system to the eigenstate |m〉 (with eigenenergy ε τ m ) of the final Hamiltonian τ = H t ( ) given the initial state |n〉, and ε ε ε = − τ m n m n , 0 is the energy difference associated with the transition probability p m|n . We will consider the work done on a system which starts in the thermal equilib- Hamiltonian Ĥ t ( ) at time t. Then the system evolves up to time τ according to some driven protocol described by the time evolution operator generated by Ĥ t ( ) at the constant inverse temperature β. The final state of the system will not be, in general, an equilibrium state. This describes the non-equilibrium dynamics that we are interested in and it is the typical scenario explored in quantum thermodynamic protocols where fluctuations theorems can be applied.

DFT-inspired methods.
We consider the Hubbard system Eq. (1) and propose methods to accurately, quantitatively estimate the average quantum work 〈W〉 produced by its driven dynamics. The key idea will be to use the KS Hamiltonian Ĥ KS as a zeroth-order Hamiltonian in a perturbation expansion scheme which converges to the exact many-body Schrödinger equation [58][59][60] . Using this expansion at its zeroth-order, gives a simple method of including interactions within a formally non-interacting scheme through the DFT exchange-correlation (V XC ) and Hartree (V H ) terms. For the static case a similar scheme has been proven very effective to largely improve results over a wide range of parameters for the description of entanglement with respect to standard perturbation schemes 60 .
The formulation of DFT for the Hubbard model that we employ is the site-occupation functional theory (SOFT) 61,62 , where the traditional density in real space , . Therefore, using the SOFT we can write the KS Hamiltonian for the Hubbard model as For obtaining the mean value 〈W〉, we suggest three related protocols (as illustrated in Fig. 1).

Figure 1.
Illustration of the application of the protocols to the Hubbard dimer. This results in four possible approximations to obtain out-of-equilibrium dynamics and related thermodynamic quantities for the Hubbard dimer. These approximations can be used depending on the desired precision and the regime of interest. The worse-case scenario is to use a standard zero-order non-interacting Hamiltonian to describe the system and the dynamics. For more precision, and depending on the regime of interest, we can use either the Kohn-Sham Hamiltonian to describe system and dynamics, or use the standard zero-order non-interacting Hamiltonian for the dynamics but refine the approximation of the initial and final energy spectra (here done up to first order perturbation, FOP): the implementations of both options have similar degrees of difficulty. Finally, we can combine the Kohn-Sham Hamiltonian for the dynamics and the FOP (or higher order precision) for the initial and final spectra.
'Zero-order' approximation protocol. We write the interacting Hamiltonian as = + ∆ˆĤ H H 0 where Ĥ 0 is a (formally) non-interacting Hamiltonian and ∆ ≡ −ˆĤ H H 0 . We will consider the case in which =Ĥ H KS 0 and compare it with the case in which Ĥ 0 corresponds to the standard non-interacting approximation to Ĥ , see Results section. Then we approximate the initial thermal state of the system as the non-interacting one, 0 0 0 , and the time-dependent evolution is calculated according to the (formally) non-interacting Ĥ t ( ) 0 up to the final time τ. The quantum work will then be estimated based on this time evolution and on the spectra of the initial and final zero-order Hamiltonians. We note that in this protocol the time dependency in Ĥ t ( ) 0 is included only through the (explicitly) time-dependent term appearing in Ĥ . For =Ĥ H KS 0 , we expect this method to reduce the magnitude of the perturbation ∆Ĥ as many-body interactions are already partially accounted in the zero-order term Ĥ KS , see Eqs (2), (3) and (4). We note that indeed the zero order term Ĥ KS of this approximation already reproduces a very important property of the many-body system, namely the ground-state site occupation distribution. We therefore expect this to produce more accurate results with respect to the standard perturbation expansion.
Protocol with first-order correction to eigenenergies. Here we propose to use the same formally non-interacting dynamics of the previous protocol, but now to include the first order correction in the estimate of the eigenenergies associated to the initial and final Hamiltonians of the system, with the corresponding correction to the approximation of the initial thermal state. As we will show, a better estimate of the eigenenergies may be important to preserve agreement with exact results for certain regimes, and especially so when the zero-order and exact Hamiltonians present qualitatively different eigenstate degeneracies, as in the case study below. While, in this paper, we will consider first-order corrections only, it should be possible to further improve accuracy by including higher order corrections to the eigenenergies.

Protocol including time-evolution effects of many-body interactions. This third protocol applies to the case in which
Here we include an implicit time-dependency in V H and V xc . These quantities are functionals of the site occupation, so time dependence can be included via the time dependence of n i . This implies a time-dependence which is local in time, i.e. that has no memory and could then not describe accurately non-Markovian processes. Nevertheless, it allows to mimic, at least in part, the variation with time of the interaction effects due to the particles dynamics. As we will show, this significantly improves results in certain regimes. The time-dependent site occupation n i (t) will be obtained by solving the system self-consistently. This protocol takes inspiration from the adiabatic LDA approximation for time-dependent DFT [63][64][65] . It may be further enhanced by improving the approximation for the eigenenergies of the initial and final time Hamiltonians, as described in the previous subsection.
We will consider two approximations for the KS exchange-correlation potential (see Results and Methods sections). We will compare the results from these protocols to the exact results and to the results obtained for same-order standard perturbation theory.

Results
We will focus on a Hubbard dimer with two electrons of opposite spin (half filling), which is known to display a rich physical behaviour [66][67][68] , including an analogue to the Mott metal-insulator transition 67,68 . Because of its non-trivial dynamics, this model is ideal as a test bed for assessing the accuracy of approximations in reproducing quantities related to quantum fluctuations and quantum thermodynamics. When the system is reduced to a dimer with half-filling, the Eq. (1) can be written in the subspace basis set We will calculate the average work along the dynamics driven by the linear time-dependent on-site potentials τ , with the parameterized initial and final values Δ 0 = 0.5J, Δ τ = 5J, at the parameterized temperature T = J/(0.4K B ). The combined values of U and τ will then determine the dynamical regime (sudden quench to intermediate to adiabatic, see discussion of Exact results). Due to the small Hilbert space associated to (5), the quantum dynamics generated by it can be exactly solved by directly integrating the time-dependent Schrödinger equation. The Hamiltonian (5) can describe various physical systems, including, but not limited to, two gate-defined quantum dots [69][70][71] , cold atoms 67 , etc.
A schematic description of how we will apply to the Hubbard dimer the protocols proposed is provided in Fig. 1. Exact results. Average work. In the panel (a) of Fig. 2 we display the exact work that can be extracted from an interacting Hubbard dimer at half-filling when operated according to the dynamics described by the Hamiltonian (5). By varying the interaction U and the time length of the dynamics τ, we can access very different regimes: from non-interacting to very strongly interacting; from sudden quench, to intermediate and to adiabatic dynamics.
Crossover between non-adiabatic to adiabatic dynamics. This system has three characteristic energy scales, U, 1/τ, and J. For the parameter region for which > ∼ U J, the crossover between non-adiabatic to adiabatic dynamics depends mostly on the interplay between U and 1/τ, and occurs when the two energy scales become comparable, that is for τ ∝ U 1/ , see behaviour of dashed-red curve in Fig. 2(a). For each given U, when τ  U 1/ the work does not depend on the time length of the dynamics, showing that, for that particular U the work has converged to its adiabatic value. For τ <  U 1/ and a given interparticle interaction, the work instead strongly depends on τ, increasing with τ up to the maximum allowed for that particular interparticle interaction. For <  U J, the J energy scale starts to influence the crossover between non-adiabatic to adiabatic regime, so that the simple relation between U and τ described above breaks down in Fig. 2(a).
The contour curves for equal average work are strikingly different between the adiabatic and non-adiabatic regime: in the U, τ plane, they can be well approximated by U = constant for the adiabatic regime, while they rapidly and almost linearly increase with τ for non-adiabatic dynamics. The behaviour of these contour curves mirrors the fact that in the non-adiabatic regime the final state of the system, and so the work, is strongly influenced by the details of the dynamics, and hence strongly dependent on the time-scale on which the time-dependent driven protocol occurs; however, by definition, in an adiabatic (energy-level transitionless) process the system remains at all times with the same energy-level occupation of the instantaneous Hamiltonian as in the initial thermal state, which means that the final state of the system -and hence the work -is completely defined, independently from the time taken by the system for going from the initial to the final state. We note that, as the Hamiltonian changes due to the driven protocol, the final state in the adiabatic regime is not an equilibrium state at the inverse temperature β.
Transition to 'Mott insulator'. The very strongly interacting parameter region > ∼ U J 5 corresponds to the two-particle equivalent of the Mott-insulator 67 , where double site occupation becomes energetically very costly. Hence for the Hubbard dimer, in this regime both double and zero site occupation are highly suppressed. As the Hilbert space available to the system dynamics reduces across the transition, we observe a corresponding decreasing of the average work that can be extracted from the system.
Entropy production and irreversibility. The entropy production Σ is defined in terms of the dissipated work in the out-of-equilibrium dynamics, 0 is the free energy variation in the protocol. It is related with an uncompensated heat, that is the energy that will be dissipated to the environment in order for the out-of-equilibrium system to get back to the thermal equilibrium. In this sense Σ quantifies the degree of irreversibility of the dynamical process at hand 22,72 . It is then instructive to look at it as the system undergoes different dynamic regimes. In general the degree of irreversibility will be related to the size of the Hilbert space available to the dynamics. For the system at hand both entering the adiabatic regime and the Mott metal-insulator transition contribute to the reduction of the available Hilbert space, and hence to the decrease of entropy production. This phenomenon shows well in the data plotted in the panel (b) of Fig. 2, where we observe a combined decrease in the entropy production as τ increases and the adiabatic regime is entered and as U increases and the system tends to 'freeze' towards the n 1 = n 2 = 1 Mott-insulator configuration.

Results from 'zero-order' approximations (standard non-interacting and Kohn-Sham based).
In this section we compare results from the protocol which uses a zero-order approximation, where the dynamics is propagated according to the Hamiltonian Ĥ t ( ) 0 and time dependency is included only through the actual driving term Δ i (t). As the zero-order Hamiltonian, Ĥ t ( ) 0 , is always formally the sum of single-particle Hamiltonians, the dynamics it generates will then correspond to the dynamics of non-interacting systems. Being formally non-interacting, this dynamics is easier to calculate numerically and would be easier to simulate and measure experimentally (in a quantum simulator) than the one originated by the many-body Hamiltonian Ĥ t ( ). We underline once more that the formally non-interacting systems in Ĥ t ( ) 0 represent physical systems (non-interacting particles) in the case of standard zero-order perturbation , while represent fictitious systems ('Kohn-Sham particles') when =Ĥ H KS 0 . In Fig. 3 we present the contour plots of the relative error, with respect to the exact average extractable work for all 'zero-order' approximations. In the panel (a) we show the results for standard perturbation theory, where the dynamics is propagated according to the non-interacting part of the many-body Hamiltonian, H NI (t). In the panels (b) and (c) of Fig. 3 we present our results for the case where the Hamiltonian is approximated to zero-order by Ĥ t ( ) KS . We approximate the exchange-correlation potential in two ways: using the pseudo local density approximation (PLDA) as proposed by Gunnarsson and Schönhammer 61 (Fig. 3, panel (b)) and using the recent parametrization to the exact correlation energy E c proposed by Carrascal et al. 68 (PAR) (Fig. 3, panel (c)). The two approximations PLDA and PAR are described in the 'Methods' section. It is known that the PLDA is not a particularly good approximation for the Hubbard model 50 , but we chose it as we wish to show that even this provides already a good improvement over standard zero-order perturbation theory (compare panel (a) and (b) of Fig. 3). This is true especially for small and intermediate values of τ: for example, for τ ≈ 0 (sudden quench), the use of PLDA instead of standard zero-order perturbation doubles the U-interval for which the relative error in the average work is below 10%, and makes it almost four times larger for τ~2/J (non-adiabatic/adiabatic cross-over region).
The parametrization to E c by Carrascal et al. 68 reproduces the exact correlation potential V c for the Hubbard dimer quite accurately; as such, Fig. 3, panel (c), is close to the best results we can expect for this system from the approximation with =Ĥ H KS 0 we are proposing. We see that now the average work is reproduced to high accuracy in most of the (U, τ) parameter region, even in parameter regions with strong many body interactions and/or corresponding to a dynamics very far from adiabaticity.
For very small τ (the sudden quench dynamics) we can reproduce the extractable work within 10% accuracy for interaction strengths larger than U = 4J. This is an enormous improvement over results obtained from the non-interacting zero-order dynamics (Fig. 3, panel (a)), where, for the same values of τ, the exact extractable work could be reproduced within 10% error only for < .  U J 0 5 . In the non-adiabatic/adiabatic cross-over region, τ ≈ J 2/ , we reproduce very well the exact average extractable work up to interactions ≈ U J 6 , while standard perturbation theory does poorly, accounting for interactions only up to ≈ U J 1 for τ ≈ 2 and up to ≈ U J 3 for τ ≈ 3. In this respect we wish to remark that even the 'lighter patch' occurring in the panel (c) in Fig. 3 within the region < <     1 5 35 still corresponds to very good accuracy, with a maximum relative error of 12% for the τ = 2 cut and of 14% for the U = 4 cut. Finally, in the adiabatic regime, results from Carrascal's parametrization still outperforms substantially standard perturbation, almost doubling the 10% relative error accuracy region, which now extends to interaction strengths of ≈ . U J 4 5 , against the limit of ≈ . U J 2 5 for standard perturbation.
We note that, at least for the Hubbard dimer, a 'zero-order'-type of approximation will always start to deteriorate as the system enters the Mott metal-insulator-type transition, and that this is independent of how well many-body interactions are accounted for in Ĥ 0 . In fact Ĥ 0 , and hence Ĥ KS , by definition, does not include a many-body interaction term formally written as = ∑ = ↑ ↓Û n n i L i i 1, , , : this leads to a spectrum where the singlet and triplet eigenstates, and , are always non-degenerate. However, in the Mott-insulator-type regime described by the actual many-body system of Hamiltonian Ĥ , only the two aforementioned states remain energetically accessible and, most importantly, they become degenerate. Why the first feature may be mimicked (e.g. this is done by the exact Ĥ KS to reproduce the exact ground state site-occupation profile) the intrinsic qualitative difference in degeneracy between the interacting and the formally non-interacting spectra determines the failure of any 'zero-order'-type of approximation in the Mott-insulator-type region, which is what we observe in Fig. 3.   Figure 3. Relative error of the mean extracted work for the 'zero-order' perturbation protocol. Contour plot of the relative error is shown as function of the dimensionless particle-particle interaction strength U/J and the dimensionless evolution time τ × J. (a) Standard zero-order perturbation theory. (b) KS-based zero-order perturbation theory using the PLDA approximation for the exchange-correlation potential. (c) KS-based zeroorder perturbation theory using an accurate parametrization for the exact E c .
We note that the large improvement provided by the 'zero-order' , KS-based approximations comes at no additional computational cost with respect to standard perturbation, as in both cases we are propagating formally non-interacting Hamiltonians.
Adding first order perturbation corrections to the initial and final energy spectra. At the end of the previous section we have discussed how, in regimes where the extent of extractable work is dominated by the spectrum and details of the system dynamics become less relevant, the KS-based 'zero-order' approximation protocol may be seriously limited, and especially so if there exist different degeneracy patterns between the exact and the 'zero-order' approximation spectra. For the Hubbard dimer this happens in the Mott-insulator-type parameter region. In this section we explore if a potential solution to this issue could be to lift this degeneracy by applying higher order perturbation to the initial (t = 0) and final (t = τ) spectra, while performing the system dynamics according to the 'zero-order' Hamiltonian. In this paper, we have consider first order perturbation (FOP) corrections, and results are provided in Fig. 4. FOP corrections to the energy spectra can be considered accurate only for relatively low interactions <  U J 1 . For these values of U we see indeed either an improvement in accuracy or, where results were already within the 10% of the exact ones, this accuracy is maintained.
For larger values of the interaction, we can give a qualitative explanation of the influence of modifying the energy spectra. Let us first consider the adiabatic regime (τ  U 1/ ): here results for the average work are dominated by the accuracy of the spectrum as the system -in a perfectly adiabatic case -would remain at any t in a thermal state characterized by the same occupation probabilities determined at time t = 0. So for < <   J U J 2 5 , as the spectra provided by the FOP are increasingly quantitatively worsening, we see that this correction reduces the accuracy of the results. However, for > ∼ U J 5 the system undergoes the analogue to the Mott metal-insulator transition, which, as discussed in the previous section, cannot be properly accounted for by 'zero-order' protocols because there is a qualitatively different degeneracy between the zero-order and the exact spectra. In this parameter region then the protocol using the FOP spectrum, quantitatively inaccurate but qualitatively correct, provides a substantial improvement over the 'zero-order' protocol, as can be seen in Fig. 4. In particular, FOP corrections are very important at very large particle-particle interaction strength U, ≈ U J 10 : here as long as U is accounted for in the eigenenergy splitting, the system freezes in the ground state and the FOP approximation is then enough to reproduce the exact result W = 0 shown in Fig. 2. In the same parameter region, without the FOP correction, the standard non-interacting zero-order approximation, completely independent from U, would predict maximum average work (W ≥ 3.3J for τ > . ∼ J 2 5/ ), while the KS-based zero-order protocols provide some improvement over this result but still predict a way-too-high average work (W ≥ 2.1J).
For > ∼ U J 1 , and for non-adiabatic and transition regime (τ <  U 1/ ), both spectrum and dynamics contribute to the average work. Here results from the FOP corrections seem to depend on how well the 'zero-order' protocol was already reproducing the exact dynamics. In particular, for the KS-based zero order protocol which uses the accurate parametrization of the exact E c (panel (c) of Fig. 4), the contribution of the quantitatively incorrect spectra from the FOP worsen the results.

Including implicit time-dependency in many-body interaction terms (without and with FOP).
So far we have considered zero-order Hamiltonians Ĥ 0 where particle-particle interactions were included at most through time-independent functionals of the initial site-occupation. However a more accurate representation of the driven system evolution should be expected by including time-dependent functionals. In this subsection we take inspiration from the adiabatic LDA and propose to include a time-dependence in these functionals by considering the same functional forms as for the static DFT but calculated at every time using the instantaneous site-occupation. The time-dependence considered is then local in time. To implement this protocol numerically, a self-consistent cycle to obtain the time-dependent site-occupation n i (t), and from there the V n t [ ( )] functionals, is necessary. Figure 4. Relative error of the mean extractable work for the first-order correction to the eigenenergies' protocol. Contour plot of the relative error is shown as function of the dimensionless particle-particle interaction strength U/J and the dimensionless evolution time τ × J. (a) Dynamics is generated by standard zero-order perturbation theory. (b) Dynamics is generated by KS-based zero-order perturbation theory using the PLDA approximation for the exchange-correlation potential. (c) Dynamics is generated by KS-based zeroorder perturbation theory using an accurate parametrization for the exact correlation energy E c .
Scientific RepoRts | 7: 4655 | DOI:10.1038/s41598-017-04478-y We illustrate this by applying the protocol to the PLDA exchange-correlation functional and focusing on the non-adiabatic and crossover regimes, τ ≤ ≤ J 0 4 / . We use as starting point the exact density at the initial time, i.e., . From this density we obtain the exchange-correlation energy . We evolve the system using this Hamiltonian and we obtain the state of the system ρˆt ( ) (1) . From this state we can calculate the next iteration for the site-occupation . Using this, we restart the same cycle calculating the E t ( ) and consequently a new Kohn-Sham . This cycle is repeated until the convergence criteria Results are shown in Fig. 5, panel (a), to be compared with the panel (b) of Fig. 3 for τ ≤ ≤ J 0 4 / . As the system exits the sudden-quench regime (τ > ∼ J 1/ ), and the site-occupation starts to respond to the dynamics, we notice a marked improvement over using time-independent functionals. For τ > 1.5/J we now achieve an accuracy of at least 10% up to interaction strength ≤ ≤ J U J 5 6 , while in Fig. 3 it was only up to ≤ U J 4 . The wavy-pattern in the contour lines for τ > . J 1 5/ reflects the system charge transfer dynamics between the two sites: the PLDA functional is unable to reproduce correctly the Mott-insulator transition so some charge transfer dynamics persists at large values of U.
When including first order corrections to the initial and final energy spectra (Fig. 5, panel (b)), we recover, and for analogous reasons, a behaviour similar to what observed in Fig. 4.

Conclusion
We have proposed a new method which uses some tools and concepts from density functional theory to study the non-equilibrium thermodynamics of driven quantum many-body systems, and illustrated it by the calculation of the average extractable work in a driven protocol. The method has the advantage of considering appropriate formally non-interacting systems (Kohn-Sham systems) to approximate the system dynamics, circumventing the theoretical and experimental problems of dealing with actual many-body interactions. It is easily scalable to large systems and can be used at different levels of sophistication, with increasing accuracy.
We have tested it on the Hubbard dimer, a two-spin system with a rich dynamics which includes the precursor to a quantum phase transition (Mott metal-insulator transition), and which can be embodied by various physical systems, including coupled quantum dots and cold atom lattices. Our results show that the proposed method reproduces the average extractable work to high accuracy for a very large region of parameter space: for all dynamical regimes (from sudden quench, to the non-adiabatic to adiabatic crossover region, to the adiabatic regime) and up to quite strong particle-particle interactions ( <  U J 6 ) our results are within 10% of the exact results.
These very encouraging results, together with the simplicity of the method make for a breakthrough in the calculation of non-equilibrium thermodynamic quantities, as the quantum work, in a complex many-body system. Future developments include the possibility to combine the method with quantum simulation techniques and an experimental implementation of quantum simulators based on this method.

Methods
Approximations for the exchange-correlation energy. The exchange-correlation energy is a functional of the site occupation density, but its functional form is unknown and needs approximations. In this work we will consider and compare results from two different types of approximations to the exchange-correlation energy for the Hubbard Model.
The first is the pseudo-LDA expression 61