Melting a Hubbard dimer: benchmarks of ‘ALDA’ for quantum thermodynamics

The competition between evolution time, interaction strength, and temperature challenges our understanding of many-body quantum systems out-of-equilibrium. Here, we consider a benchmark system, the Hubbard dimer, which allows us to explore all the relevant regimes and calculate exactly the related average quantum work. At difference with previous studies, we focus on the effect of increasing temperature, and show how this can turn the competition between many-body interactions and driving field into synergy. We then turn to use recently proposed protocols inspired by density functional theory to explore if these effects could be reproduced by using simple approximations. We find that, up to and including intermediate temperatures, a method which borrows from ground-state adiabatic local density approximation improves dramatically the estimate for the average quantum work, including, in the adiabatic regime, when correlations are strong. However at high temperature and at least when based on the pseudo-LDA, this method fails to capture the counterintuitive qualitative dependence of the quantum work with interaction strength, albeit getting the quantitative estimates relatively close to the exact results.


Introduction
In the last decades, we have witnessed remarkable progress in density functional theory (DFT) with the development of various tools to study many-body quantum systems [1][2][3][4][5]. In particular, time dependent DFT (TD-DFT) [6] allowed to go beyond ground-states and opened a new path [7][8][9] in the ab-initio formulation for electronic transport [10][11][12], electronic excitations [9], and calculations of thermodynamical properties of solid state systems [13,14]. Currently, with the urge to describe accurately realistic devices for quantum technologies, theoretical and computational physicists have been devoting great efforts to improve the treatment of out-of-equilibrium systems [12,15].
As in the classical world, thermodynamics will impose limits in the fabrication and operation of devices for quantum technologies [16][17][18][19][20]. The descriptions and the laws as formulated within the conventional thermodynamics are no longer valid at the scales where these technologies are being developed [21]. At this level energy fluctuations become important, and quantities such as work, heat, and entropy production are treated as stochastic variables [22]. The quantum thermodynamics research has been fueled recently by experiments, on small and noninteracting systems [23][24][25][26][27][28]. However, the experimental study of the thermodynamics of quantum many-body systems remains a challenging task, and even theoretical studies could require an enormous computational power. The investigation of the thermodynamics of the emergent collective phenomena in quantum many-body systems is, without a doubt, a fascinating subject. In this context, some discussions about quantum thermodynamic properties in out-of-equilibrium systems have been reported (see e.g. [29][30][31][32][33][34][35]). In particular, a method to calculate quantum thermodynamic properties of interacting systems subject to driving fields was presented recently in reference [14], where it was applied to the calculation of the average quantum work in a Hubbard dimer driven by a time-dependent external potential. Such method uses an approximated framework based on tools from DFT, where the many-body problem dynamics is mapped onto the non-interacting dynamics of Kohn-Sham Hamiltonians. As it is shown in [14], DFT can offer a reliable approach for estimating thermodynamical properties in out-of-equilibrium systems: this method provides good accuracy when a suitable approximation for the so-called exchange-correlation potential is used.
Using the approach in [14], the present paper aims to examine in more detail the use of different levels of approximation for the proposed DFT-protocol. In particular, we focus on the effect of increasing temperatures, and discuss the competition between temperature, coupling regime and evolution time on the average quantum work W extracted from a two-qubit system subject to a non-linear dynamics. We will compare the exact W extracted from the interacting system with its estimate from the noninteracting picture, and from approximations based on the local density approximation (LDA). We will examine the improvement in the calculations when the dynamical effects are included into the exchange-correlation functional by means of an adiabatic LDA-inspired method (ALDA-i) and explore the limits of a ground-state formalism as the system temperature rises.

The driven Hubbard dimer
The two-site Hubbard Hamiltonian is widely used for benchmarking density-functional approaches for strongly correlated systems [14,36,37]. It also provides a simplified model to study two qubit systems in all coupling and dynamical regimes. The size of the problem allows for exact solution as well as, potentially, for experimental verifications. The Hubbard dimer has in fact been simulated using different types of systems, such as quantum dots [38][39][40] and cold atoms [41], and could be simulated using small molecules driven by NMR techniques.
An interacting Hubbard dimer driven by a nonhomogeneous potentialV (t) is described by the HamiltonianĤ (t) =K +Û +V (t). (1) represents the kinetic energy with hopping parameter J; c † iσ (c iσ ) is the creation (annihilation) operator acting on sites 1 and 2 with spin σ =↑, ↓; represents the electrostatic Coulomb interaction of strength U , and number operatorn iσ =ĉ † iσĉ iσ , with i = 1, 2; and, finally, is the external, non-homogeneous, time-dependent potential. Fig. 1. Non-linear sinusoidal potential from equation (5). During the dynamics, the driving potential promotes a flux of electrons from site 1 to site 2. The site-occupations n1 and n2 at t = τ will depend not only on the parameters A0 and Aτ of equation (5), but also on the coupling U , the temperature kBT and the hopping parameter J.

Non-linear dynamics and quantum work
In our analysis, we consider a two-site Hubbard model driven by a non-linear potential with sinusoidal form where j = 1, 2 labels each site. As indicated below, ω 4τ corresponds to a period of 4τ .
In this paper, we consider the parameters A 0 = J and A τ = 7J and explore coupling regimes from fully noninteracting (U/J = 0) to the strongly coupled (U/J ≈ 10). The potential is plotted in Figure 1.
We will focus on the extraction of quantum work W during a dynamical process at finite temperature which starts at t = 0 and ends at t = τ = π/(2ω 4τ ). At t = 0, the system is described by the equilibrium thermal state where Z(t) = Tr e −Ĥ(t)/k B T is the instantaneous partition function, k B is the Boltzman constant and T is the temperature. Apart from the driving field, the system is considered isolated. At the end of the process, ρ(t = τ ) is not, in general, described by an equilibrium state. The extracted quantum work W is defined in terms of the work probability distribution P (w), a thermodynamical quantity yielding information about the energetic transitions and non-equilibrium fluctuations taking place while the quantum system evolves through a driven dynamics. We define W according to [42] as where P (w) is given by Here, p n (t = 0) denotes the population of the nth eigenstate |ψ n (t = 0) of the initial Hamiltonian; p m(t=τ )|n(t=0) defines the probability to find the final system ρ(t = τ ) in the mth eigenstate |ψ m (t = τ ) of the final Hamiltonian H(t = τ ), when starting in |ψ n (t = 0) ; and ∆ n,m is the energy difference between eigenenergies E n (t = 0) and E m (t = τ ).
Let us examine energy, temporal and thermal scales of the problem. The competition between J, U and the strengths of the external potentials A 0 and A τ dictates the first. In the weakly coupled regime U J, one might expect a non-interacting description to be adequate for representing the system. The changes in the instantaneous densities n j (t) are ultimately determined by a delicate balance with the external time-dependent potential. For increasing values of U , the external potentials, which transform the second site in an electron sink, enter in competition with the Coulomb repulsion.
The parameter τ controls the speed at which the particles are driven. For τ × J 1, the dynamics corresponds to a sudden quench: the evolution is so fast that the system does not have enough time to adapt. At the opposite limit, τ × J 1, the evolution can be considered adiabatic. The dynamic regimes are illustrated in Figure 2, where we show the instantaneous densities n 1 (t) (red tones) and n 2 (t) (blue tones) for τ × J = 0.5 (dark shades) to τ × J = 10 (light shades), U = 2J, and k B T = 2J. In the sudden quench regime the system dynamics is so out of synchronicity with the applied potential that the final value of the site occupations will strongly fluctuate with τ (see inset). As the adiabatic regime sets in, the average value reached by the site occupation probabilities becomes much better defined, though a finer oscillation persists, due to the competing effects on electronic transport of external potential and Coulomb repulsion. The amplitude of this oscillation decays with increasing τ .
The temperature introduces an additional dimension to the system parameter space. Figure 3 depicts the populations p n (n = 0, 1, 2, 3 for the dimer at half-filling) of the initial thermal state ρ(t = 0) as a function of the temperature T (in units of J/k B ) for U/J = 0 (noninteracting) and U/J = 10 (strongly coupled). At very low temperatures, the initial thermal state corresponds to the pure ground-state ρ(t = 0) T →0 = |ψ 0 ψ 0 |. Increasing temperature allows the system to populate different levels E n (t), giving rise to a larger number of potential transitions during the dynamics. At very high-temperatures, all states are equally contributing to the initial thermal state and a non-interacting picture arises.

Zero-order DFT protocol
In this section, we review the methods proposed in [14] to calculate the thermodynamic properties of an out-ofequilibrium system.
The "zero-order DFT protocol" [14] approximates to zero-order the interacting Hamiltonian withĤ KS (t). H KS (t) describes a formally non-interacting system which reduces to the standard non-interacting system when many-body interactions are switched off. For finite manybody interactions, the single-particle potential in H KS (t) accounts for many-body correlations in such a way that the electronic density of the many-body system is, in principle, exactly reproduced. However, other properties of the many-body system are not guaranteed to be reproduced or even well-approximated, and especially so in the strongly correlated regime. One of the aims of this paper is indeed to explore how close the quantum work of the KS system is to the many-body quantum work in the various regimes. The advantage of the zeroorder DFT protocol over the perturbation theory is that H KS (t) is formally non-interacting, but contains some of the interaction effects. This method has proven useful in estimating both entanglement [43] and quantum thermodynamic quantities [14]. In respect to the latter, an additional advantage of this method is that it formally avoids having many-body operators acting during the quantum evolution, which simplifies both the numerical simulations as well as potential experiments: these, in appropriate regimes, could be designed based on this approximation [14].
Within the "zero-order DFT protocol", the initial thermal state is approximated as Afterwards the system evolves according to the Kohn-Sham HamiltonianĤ KS (t), and at t = τ an approximation for the final state ρ(t = τ ) can be extracted. By means of this protocol, we can obtain an estimate for the average quantum work discussed in Section 3.
For the driven dimer, the interacting systemĤ(t) is mapped into a non-interacting system according tô where the last term is the Kohn-Sham potential, defined asV where n(t) = (n 1 (t), n 2 (t)), and V H j = U 2 n j is the Hartree potential accounting for the classical electrostatic repulsion.
The first term on the right-hand side of equation (12) is an effective single-particle potential which accounts for exchange and correlation effects, the so-called exchangecorrelation potential [44]. In general, the exact functional form of the exchange-correlation potential is unknown, requiring approximations. In this paper we will consider local density-type of approximations, as described below.

Zero-order DFT protocol with Pseudo-LDA (p-LDA)
We consider time-independent V xc j and V H j . They are calculated at time t = 0, and V xc j through the pseudo-local density approximation [45]. They are given by Although expression (13) comply with the main idea of an ab-initio approach to approximate the exchangecorrelation potential, it is known to be a rough approximation for Hubbard-like Hamiltonians. As discussed in reference [45], LDA approximations based on the manybody solution of the homogeneous Hubbard Hamiltonian, such as the Bethe Ansatz Local Density Approximation (BALDA) [46], offers a more appropriate approach.
For the purposes of the present paper, considering p-LDA instead of BALDA suffices for a qualitative analysis of out-of-equilibrium systems: our aim here is to test the improvement in density-functional calculations with the inclusion of temporal dependences in the exchangecorrelation functional. Therefore, our comparison of static and adiabatic p-LDA illustrates the situation in which one does not have at hand the best approximation for V xc .

Zero-order DFT protocol with adiabatic local density approximation
So far, we have considered zero-order HamiltoniansĤ 0 where particle-particle interactions were included at most through time-independent functionals of the initial siteoccupation. However, a more accurate representation of the driven system evolution should be expected by including time-dependent functionals. We take inspiration from the ground-state adiabatic LDA (ALDA) [6] and include a time-dependence by considering the same functional forms as for the static DFT but calculated at every time using the instantaneous thermal site-occupation. This time-dependence is local in time.
To implement this protocol numerically, we use a self-consistent cycle to obtain the time-dependent siteoccupation n j (t) at all times, and from there the corresponding functionals V H j [n j (t)] and V xc j [n j (t)] at all times. We iterate the protocol by running the dynamics several time, until convergency for the site occupation at all times is reached. Figure 4 depicts the flowchart corresponding to the iterative protocol. In details: we use as starting point the exact density at the initial time, i.e., n    j (t), t]. We evolve the system using this Hamiltonian and we obtain the state of the system ρ (1) (t). From this state we can calculate the next iteration for the site-occupation n (1) j (t) = Tr ρ (1) (t)n j . Using this, we restart the cycle. This cycle is repeated until the convergence criteria In the present paper we use the pseudo-LDA approximation as a base to implement this time-dependent approach.

Temperature effects on average quantum work, exact results
In this section, we analyze how the exact extracted quantum work W exact is modified by increasing the temperature. For the Hubbard dimer and 0 ≤ U/J ≤ 10, the gap ∆E between ground and first excited states is in the range 0.39J ≤ ∆E ≤ 2.82J, depending on the value of U . To consider all regimes of interest, we then consider the three temperatures: T = 0.2J/k B , for which excited states are barely populated; T = 2J/k B , for which k B T is comparable to ∆E and the excited states start to be populated; and T = 20J/k B for which the thermal energy is the largest energy in the system and all states become comparably populated (see Fig. 3).
The results for W exact are shown in Figure 5, right column, where the contour lines for the quantum work are plotted against the Coulomb interaction strength U and the evolution time τ . The hopping parameter J is our unit of energy. It can be observed that for k B T ∆E work can be mainly extracted in the low U -large τ region (panels (b) and (d)), while when k B T >> ∆E the situation is dramatically different, and the region with large U becomes the most productive region (lower panel).
We explain this behavior as follows. In the sudden quench regime (small τ values, τ × J < 1), the system will not evolve significantly, ρ(t = 0) ≈ ρ(t = τ ), and the quantum work will be in the great part determined by the characteristics of the eigenstates and spectra of the initial and final Hamiltonians, which are the same at all T . This explains why the variation of work production in this regime is mostly unaffected by temperature. However the physics becomes very different as τ increases, as the system becomes sensitive to the driving field, and eventually enters the adiabatic regime.
For τ × J > 1 and k B T ∆E, the dominant contribution to the system thermal state comes from the ground state. The dynamic of this will be affected by two competing forces: on the one side the driving field which drives the state towards double occupation in site 2; on the other the Coulomb repulsion which tends to decrease double occupation and, for very large U , would eventually lead to the precursor of the Mott-insulator phase transition. For increasing U , the reduced efficacy of the applied field means that the field can extract increasingly less work from the system, as it is clearly observed in Figure 5, right column, panels (b) and (d). This is a noticeable effect of many-body interactions on the quantum work. We note that the sign of the many body effect (reduction instead of increase of field efficacy) is due to the particular dynamics chosen, for which Coulomb repulsion and driving field are in competition. We would expect instead many-body interactions to help the process in the reverse dynamics: Coulomb interactions should help pumping work into the system when the system is driven from an initial state that favors double occupation towards a final state with a more homogeneous site potential. An example more subtle than this reverse dynamics will be given below.
This picture for k B T ∆E is supported by the contour plots of the site occupation density n 1 (τ ) (Fig. 5, panels (a) and (c)). We see that for τ × J 4 and small U , n 1 (τ ) is depleted by the action of the field well-below half-filling. As U increases and the system is not yet adiabatic (compare to the red adiabatic line in the panels), the Coulomb repulsion strongly affects n 1 (τ ) by increasing its value. As τ increases further and the system settles into the adiabatic regime, the impact of the applied field increases, giving rise to the decaying ripples already observed in Figure 2.
Let us now focus on the bottom panels of Figure 5. Here k B T >> ∆E and k B T is also larger than the maximum difference in the site potentials induced by the applied field. Because of this, the effect of the dynamic regime becomes much less striking, as can be noted by the much reduced range of extracted work and site occupation variations (compare the scale of the contour lines between the bottom panels and the ones above them). In addition, for the parameters considered, all states (ground state and excited ones) give comparable contributions to the thermal state. Also for all values of U considered the inequality k B T ≥ 2U holds, so that the effect of interactions becomes less important. Yet, quantitatively subtle, but qualitatively strikingly different behaviors from the ones observed in the upper panels, occurs and due to the presence of many-body interactions.
Counterintuitively and contrary to lower temperatures, the effect of increasing many body interactions is now to help the applied field to deplete site 1, as shown by the corresponding density contours. In turn this enhances the field efficacy and hence the work that can be extracted from the system. Maximum work can then be extracted in the adiabatic regime for large U values (the red line within each figure indicates the transition region between the sudden quench and the adiabatic regime). In fact the adiabatic regime is the one in which the least entropy is produced, and hence the system is able to produce the largest work. This is confirmed by the results shown in panel (f) of Figure 5.
The unexpected behavior of site occupation and average work with increasing Coulomb interaction arises from the subtle interplay between the evolution of the character of the eigenstates driven by the applied field, and the substantial occupation of higher energy states in the initial thermal state. When the evolution is adiabatic, the final Hamiltonian eigenstates will inherit the occupation from the initial thermal state. Strong Coulomb repulsion at t = 0 implies that eigenstates with strongly asymmetrical occupation are pushed further up in energy so that, for the same temperature, at time t = 0, higher energy states become less populated for large U 's than for lower values of U .
At strong coupling, the applied field must drive the system through an anticrossing so that the final ground state at time τ may give rise to the strongly asymmetric site occupation expected by a system dominated by a step potential. So, in the adiabatic regime this asymmetric state will have an increased weight at t = τ with respect to what it had at t = 0 due to the higher population of the ground state at time t = 0. The effect of this overall process is to transfer population to the second site.
At lower U 's, for U values less or comparable to the step in the potential at t = 0, the initial eigenstates are closer in energy as they are not too influenced by the Coulomb repulsion. Because of this, for the same temperature, the occupation of the highest energy state will be larger than in the presence of strong coupling. After the time evolution, this state corresponds to a distinctively asymmetric site occupation, which favors site one. As the state maintains the high thermal population acquired at time t = 0, the effect of the overall process is to transfer population to the first site. Note that, as the highest excited state is involved, this process would be negligible at low and intermediate temperature.
The combination of this high-temperature opposite population transfer for weak and strong coupling leads to the behavior observed, where, at high temperatures, increasing Coulomb repulsion favors asymmetry in the t = τ site occupation.
In the rest of the paper we will focus on the T = 2J/k B and T = 20J/k B cases, and study how various approximations capture the corresponding qualitative and quantitative changes in the production of quantum work.
6 Accuracy of zero-order DFT-inspired protocols 6.1 Intermediate-temperature regime, k B T = 2J In Figure 6, we consider the effect of approximations when k B T = 2J. In the upper three panels we show how climbing the ladder of "zero-order" protocols -from the standard non-interacting, and through the DFT-inspired approximations -improves the estimate of the average quantum work. From the top to the penultimate panel, we plot the approximated extracted work (left column) and its relative error with respect to W exact (right column) calculated from: a completely non-interacting Hubbard dimer (U = 0); an interacting Hubbard dimer approximated by a Kohn-Sham Hamiltonian within the pseudo-LDA approximation; and finally an interacting Hubbard dimer approximated by a Kohn-Sham Hamiltonian within an ALDA-inspired approximation and based on the pseudo-LDA (see Sect. 4.2).
Let us first concentrate on the non-interacting results (upper panels of Fig. 6). As should be expected, in this approximation the work is independent from U , and simply increases with τ as the system becomes more adiabatic. When compared to panel (d) of Figure 5, the non-interacting quantum work is qualitatively dramatically different from the corresponding exact one, with the exception of very small U and especially in the adiabatic regime ( Fig. 6 panel (b)).
When turning to the zero-order p-LDA results W p−LDA , Fig. 6 panel (c), we observe that now the exact work is reproduced up to U ≈ 2J both qualitatively and quantitatively. However, at larger U values, as correlations become stronger, the match fails even qualitatively, compare panel (c) of Figure 6 with panel (d) of Figure 5.
While only a relatively minor improvement is noted between non-interacting and zero-order pseudo-LDA results, a remarkable improvement is observed when including time-dependent correlations through the ALDAinspired scheme. Notably the bulk of this improvement is in the adiabatic parameter region, where non-Markovian effects are not important, in line with the fact that the ALDA approximation does not include memory effects. We stress though that in this intermediate-temperature regime the system state is clearly a mixed thermal state, while the scheme is based on an ALDA designed for pure states at zero temperature, so our result are non-trivial as in principle ALDA could badly fail for thermal states.
When considering W ALDA (Fig. 6 panel (e)), we notice that it qualitative reproduces the main features of the behavior of the exact work, and at all U values. In addition the range of variation of W ALDA quantitatively closely match the one of W exact in most regions, at difference with W p−LDA (compare scales at the bottom of the related panels). As this is the case also at large U 's, we consider this as a confirmation that this zero-order DFT-inspired approximation fairly treats correlation up to intermediate temperatures even when strong, at least for the system at hand. This approximation reproduces the exact work even quantitatively in the whole region corresponding to high extractable work (W exact > 5), also a very good result. We observe that the contour lines of panels (e) and (f) in Figure 6 present oscillations and we will come back to this later in the paper.
In the bottom panels of Figure 6 we consider the ALDAinspired estimate for the site-one particle occupation at time τ (left) and show its relative difference with the corresponding exact result. We note that the site occupation has, on average, the best (worst) agreement in the same parameter region where the average quantum work is closely (farthest) reproduced. Figure 7 shows the same quantities as Figure 6 but for high temperature (T = 20J/k B ). By comparing the exact work (Fig. 5, panel (f)) with the approximate estimates, we immediately notice that all the approximations qualitatively reproduce the low τ regime (τ × J 2) at all U values (with the DFT-inspired approximations improving slightly over the non interacting one). However at the same time they all fail to reproduce the behavior, even qualitatively, at larger τ 's and U/J 1, as the regime becomes adiabatic and the interactions increase beyond perturbative. Interestingly enough, in this parameter region the corrections to the contour lines provided by   the correlations introduced through the zero-order DFTinspired approaches go in the opposite direction with respect to the exact behavior. The extent of extracted work is though comparable to the exact one, although its range is slightly smaller. The above picture is well summarized in the second column of Figure 7, where the percentage variation with respect to the exact work is plotted, and the region where the exact work is worst reproduced indeed corresponds to the adiabatic, highly interacting parameter region. We note that, because the value range of approximated and exact work are similar, the error never goes below 30%.

High-temperature regime
Small differences between the results from the three approximations may be accounted for by noticing that the Hamiltonians remain different in the three cases, the more different the higher the value of U . We expect then small differences to appear for large U values, as indeed is the case. Panels (g) and (h) of Figure 7 show the ALDA-inspired estimate for the site-1 occupation n 1 (τ ) and its relative error with respect to the exact values, respectively. Panel (g) demonstrates that in this case, apart from the suddenquench regime, the estimates from the ALDA-inspired zero-order approach behave qualitatively opposite to the exact behavior (Fig. 5, panel (e)) when parameters are varied.
We attribute this opposite behavior to the fact that in respect to the exact case, the instantaneous eigenvalues for the ALDA-inspired case are not as sensitive nor qualitatively nor quantitatively to the value of U . The gaps between higher excited states are poorly reproduced: any formally non-interacting Hamiltonian leads, for the Hubbard dimer, to a degeneracy between first and second excited states, while the widening of the gap between second and third excited states with increasing U is underestimated by the ALDA-inspired estimates. This implies that the delicate interplay between thermal occupation and interaction strength at t = 0, at the origin of the exact system behavior, is not reproduced.
Interestingly, due to the small variation of the siteone density values, quantitatively the overall picture is much better, as the relative error ( Fig. 7 panel (h)) is actually very good in most parameter areas. Comparing this with the corresponding relative error for the average quantum work (panel (f)) confirms that there is a very good correlation between the quantitative accuracy of the ALDA-inspired estimates for local thermal density and average quantum work. We note that, independently from U and T , in the adiabatic regime, after an initial transient, the density oscillates around a mean value. This explains the oscillations observed in panel (e) of Figures 6 and 7. These oscillations qualitatively reproduce the ones observed in the exact density, albeit with a different period e.g. compare to the grey curve in panels (e) and (f) of Figure 8. This explains the oscillatory patterns observed in the panels (f) and (h) of Figures 6 and 7. As in the exact case, see Section 3, the oscillations are due to the interplay between the Coulomb repulsion and the attractive potential generated by the driving field, that induce transport in opposite directions.
Numerically, strong interactions (and the consequent "stiffness" of the system) slows down the convergency for the site occupation in the self-consistent cycle necessary to obtain the time-dependant KS potentials. This is most accentuated at low-intermediate temperatures as can be observed by looking at panels (c) and (e) of Figure 8, where a much larger number of iterations is necessary to obtain the same level of convergency.

Conclusions
We have presented a study for the temperature dependence of the average quantum work that can be extracted from a driven Hubbard dimer at half-filling, for varying evolution time and Coulomb interaction strength parameters. The driving field was set up as to oppose, as the system evolved, the effect of inter-particle Coulomb repulsion, by building a potential step between the two lattice sites at least twice as large as the maximum Coulomb interaction considered in this work. This competing behavior has allowed to uncover subtle and counterintuitive interplay between temperature, external driving field, and many-body effects on the production of quantum work. On the one side, at low and intermediate temperatures and all other parameters the same, increasing Coulomb interaction strength indeed opposes the effect of the applied field, and, by doing so, reduces its efficacy and hence the work produced. On the other side we discovered that, at high temperatures and medium and large system evolution times, Coulomb interaction strength favors the action of the applied driving field, and, consequently increases the production of work. We explain this counterintuitive behavior by the subtle interplay between interaction-dependent changes in the instantaneous eigenstates of the time-dependent Hamiltonian combined with the variation of their thermal occupation with temperature. We also explore the possibility of approximating the exact work and the thermal site occupation by using zeroorder DFT-inspired approximations. We find that quantitatively the ALDA-inspired approximation, based on the pseudo-LDA exchange-correlation potential, behaves well in most of the parameter space and at all temperatures considered. However, qualitatively the picture is more complex, and in particular even the ALDA-inspired approach, at least when based on the pseudo-LDA, fails to reproduce the subtle interplay between Coulomb interaction and temperature which, in the exact case, leads to the depletion of site-one for increasing Coulomb interaction at high temperature.