Thermal density functional theory approach to quantum thermodynamics
Antonio Palamara Francesco Plastina Antonello Sindona Irene D'Amico




















corresponding author: irene.damico@york.ac.uk


§ ABSTRACT 
Understanding the thermodynamic properties of many-body quantum systems and their emergence from microscopic laws is a topic of great significance due to its profound fundamental implications and extensive practical applications. 
Recent advances in experimental techniques for controlling and preparing these systems have increased interest in this area, as they have the potential to drive the development of quantum technologies. 
In this study, we present a density-functional theory approach to extract detailed information about the statistics of work and the irreversible entropy associated with quantum quenches at finite temperature. 
Specifically, we demonstrate that these quantities can be expressed as functionals of thermal and out-of-equilibrium densities, which may serve as fundamental variables for understanding finite-temperature many-body processes. We, then, apply our method to the case of the inhomogeneous Hubbard model, showing that our density functional theory based approach can be usefully employed to unveil the distinctive roles of interaction and external potential on the thermodynamic properties of such a system.


Thermal density functional theory approach to quantum thermodynamics
    Irene D'Amico
    September 9, 2024
====================================================================











§ INTRODUCTION
 
Density-functional theory (DFT) <cit.> and its time-dependent (TD) extension <cit.> are powerful and well-established methods for studying the electronic properties of interacting many-body systems at zero temperature, with DFT providing comprehensive access to ground state properties and TDDFT extending this capability to include the prediction of excited states.
Thermal density-functional theory (ThDFT), introduced by Mermin <cit.>, extends the Hohenberg-Kohn (HK) framework of DFT <cit.> to address the electronic properties of many-body systems under conditions where accounting for finite temperature effects is indispensable <cit.>.

Renewed attention in ThDFT has been driven by the study of thermal properties of out-of-equilibrium interacting quantum systems, which represents a major focus in quantum thermodynamics (QT).
This field is rapidly evolving thanks to advances in preparing and coherently controlling quantum systems at the microscopic scale, enabling experimental verification of fundamental properties such as fluctuation theorems <cit.>. 
These developments also hold potential for new quantum technologies based on complex quantum systems (see e.g. <cit.>).
A major goal in QT is to understand the role of purely quantum features, such as coherence and correlations, in thermodynamic processes.
This notably includes work processes, where work is extracted from or performed on a quantum system, and the generation of irreversible entropy <cit.>. 



In this realm, the emphasis is rapidly shifting towards coupled many-body systems; see, e.g,  <cit.>. Indeed, recent studies have demonstrated that particle interactions can enhance the efficiency of quantum heat engines <cit.>.
On the theoretical side, addressing finite-temperature quantum many-body systems poses significant challenges, often requiring approximations to manage their complexity. 


In this context, and drawing inspiration from previous works <cit.>, our objective is to establish a robust theoretical framework for applying DFT to the study of non-equilibrium thermodynamics in quenched interacting many-body systems. Our formalism focuses on the canonical ensemble, of particular relevance for QT thermal machines.
Specifically, we demonstrate that the thermal and out-of-equilibrium densities form the basis of an ab initio framework for deriving thermodynamic properties of quantum systems that experience a sudden quench.
A key advantage is that the out-of-equilibrium thermodynamics of interacting many-body systems can be effectively investigated using the Kohn-Sham (KS) approach to DFT <cit.>. 
The strength of this approach lies in its ability to evaluate, in principle exactly, the thermal densities by mapping the original interacting many-body system onto a fictitious non-interacting one.
We thus establish a general framework for the DFT approach to QT, specifically for the canonical ensemble, and validate it through the analysis of the quenched inhomogeneous Hubbard model. 
Accordingly, we first compare results from exact diagonalization for small systems with those obtained using our finite temperature KS  mapping, and then extend our analysis to larger systems. 
This extension allows us to clarify how interactions influence thermal densities and, consequently, the work performed.

The reminder of the paper is organized as follows. 
In Sec. <ref>, we recall the Mermin-Hohenberg-Kohn (MHK) theorem <cit.>, adapting it to the context of lattice Hamiltonians in closed quantum systems. 
In Sec. <ref>, we review the fundamental concepts of out-of-equilibrium QT, showing that for work protocols of infinitesimal duration, the probability distributions of work and irreversible entropy production are completely determined by the finite-temperature equilibrium densities of the pre-quench and post-quench Hamiltonians. 
In Sec. <ref>, we recall the Mermin-Kohn-Sham (MKS) equations <cit.>, employing the finite-temperature KS mapping to calculate the thermodynamic quantities of interest. 
In Sec. <ref>, we present a method for calculating the thermal densities of a system of indistinguishable particles within the canonical ensemble using the KS mapping.
In Sec. <ref>, we apply our theoretical framework to short Hubbard chains, solving the problem both exactly and via the finite-temperature KS mapping to validate the robustness and accuracy of our DFT-based approach. 
We then extend our DFT-based analysis to the Hubbard model with a larger number of sites. Finally, we draw our conclusions in Sec. <ref>.



§ THERMAL {A}-FUNCTIONAL THEORY

In its original formulation, DFT is founded on the HK theorems <cit.>, which were later generalized to finite temperatures by Mermin <cit.>.
These theorems establish a one-to-one correspondence between the external potential v(𝐫) (or v(𝐫) - μ, where μ is the chemical potential), the quantum (equilibrium) state of the system, and the ground state (or thermal) electron density n(𝐫) (or n^β(𝐫)).
More recently, it has been demonstrated that the HK theorems at zero temperature are still valid for quantum systems whose Hamiltonian is defined on a lattice <cit.>, with some limitations <cit.>.
Similarly, it can be shown that this extension also holds true for finite-temperature closed quantum systems. 

To demonstrate this, let us consider a closed quantum system governed by a Hamiltonian ℋ̂ of the following form:

    ℋ̂[{λ_i}]= ℋ̂_0+ 𝒱̂_ext[{λ_i}]
    =ℋ̂_0 + ∑^ℒ_i=1λ_i Â_i,

where ℋ̂ is quite general and can describe spin chains, fermionic, or bosonic systems. 
ℋ̂_0 is the `universal' HK Hamiltonian, and 𝒱̂_ext is the external potential, controlled by a set of ℒ parameters {λ_i} that multiply the local operators {Â_i}.

Let us denote the thermal expectation value of Â_i as 

    a^β_i := {ρ̂^th_β[{λ_i}]Â_i},

where

ρ̂^th_β[{λ_i}]= 
e^-βℋ̂[{λ_i}]/𝒵[{λ_i}]

represents the Gibbs state and 𝒵[{λ_i}] = {e^-βℋ̂[{λ_i}]} expresses the canonical partition function.
By definition, each a^β_i is a function of the parameters {λ_i}.
It can be shown that exactly one parameter set {λ_i} corresponds to a given mean value set {a^β_i} (see appendix <ref>).

Since {a^β_i} uniquely determines {λ_i}, which in turn determines ρ̂^th_β, the thermal state is also a unique functional of {a^β_i}.
This relationship highlights the direct connection between observable mean values and the underlying thermal state, emphasizing the role of {a^β_i} as the fundamental descriptors of the system:

    {λ_i}{a^β_i}ρ̂^th_β≡ρ̂^th_β[{a^β_i}].


By the one-to-one relations established in Eq. (<ref>), the thermal HK theorem, as demonstrated by Mermin <cit.>, remains valid for closed quantum systems defined by the  Hamiltonian (<ref>). 
Consequently, the free energy 

    ℱ[ρ̂]:={ρ̂(ℋ̂+lnρ̂/β)},

minimized by the equilibrium Gibbs state, can be expressed as a unique function of {λ_i} or, equivalently, as a unique functional of {a^β_i}:

    ℱ[{a^β_i}] = Ω[{a^β_i}]+ ∑^ℒ_i=1λ_i a^β_i.

Here, Ω[{a^β_i}] represents the generally unknown `universal'  functional <cit.>, associated with the {λ_i}-independent Hamiltonian ℋ̂_0.
Similarly to the zero-temperature case <cit.>, Eq. (<ref>) allows us to interpret the HK theorem as an expression of duality in the sense of a Legendre transform. 
This duality relates the thermal expectation values {a^β_i}, which play a role analogous to the thermal electron densities n^β(𝐫), to the set of work parameters {λ_i}, which fulfill a role analogous to the external potentials v(𝐫).

A special case of the Hamiltonian (<ref>) corresponds to the following single-parameter scenario:

    ℋ̂[λ]= ℋ̂_0+ λ∑^ℒ_i=1Â_i,

where the HK theorem remains valid. 
However, to ensure a unique mapping between λ and 

    ⟨Â|≡⟩{ρ̂^th_β[{λ_i}] ∑_i=1^ℒÂ_i},

it is necessary to require that [ℋ̂_0, ∑_i=1^ℒÂ_i] ≠ 0 <cit.>.

Another important tool is offered by the Hellmann-Feynman (HF) theorem <cit.>,  
which establishes a relationship between the first derivative of the free energy with respect to the i-th external parameter and the i-th thermal density:

    ∂ℱ/∂λ_i = {ρ̂^th_β[{λ_i}] ∂ℋ̂/∂λ_i} 
    ={ρ̂^th_β[{λ_i}]Â_i}=a^β_i.

This equation will be especially useful in the subsequent sections, particularly in the context of thermal sudden quenches, where thermal expectation values serve as fundamental variables for deriving the thermodynamic quantities of interest.




§ THERMAL {A}-FUNCTIONAL THEORY APPROACH TO QUANTUM THERMODYNAMICS
 
ThDFT can be effectively employed to extract information about the out-of-equilibrium thermodynamics of a closed quantum system by leveraging its ability to handle thermal and quantum fluctuations.
To establish the formalism, we first review some key concepts in QT.
Our focus is on a closed quantum system that has been driven out of equilibrium by a unitary quantum process.  
Specifically, we consider a generic closed quantum system characterized by the Hamiltonian (<ref>), which depends on a set of time-dependent work parameters {λ^t_i}. 
The system is initially in equilibrium with a bath at inverse temperature β. 
In this configuration, at t=0, the set of work parameters {λ^0_i} defines the thermal state, 
represented by the Gibbs density operator ρ̂^th_β[{λ^0_i}]. 
After the system is decoupled from the bath, it undergoes unitary dynamics, from t=0 to t=τ,  governed by the time evolution operator 𝒰̂(τ,0). 
This evolution is driven by a protocol that changes the set of work parameters from {λ^0_i}, with corresponding mean values {a^β 0_i},  to {λ^f_i}, with corresponding mean values {a^β f_i}, over the finite time interval τ.



 §.§ Probability distributions of work and irreversible entropy production

The work performed or extracted during the protocol is not an observable and cannot be represented by a Hermitian operator. Rather, work is a stochastic variable characterized by a probability distribution, which is determined by performing two projective measurements at the initial and final times, respectively <cit.>. 
These measurements involve the instantaneous eigenbasis of the system Hamiltonian

denoted {ϵ_n(t), |n(t)⟩}.
The probability distribution of work (PDW) is defined as follows:

    P_τ(w)=∑_nmp_n(0)p_n → m(τ)δ(w-w_mn),

where w_mn = ϵ_m(τ) - ϵ_n(0) is the work performed in a single realization, p_n(0) ≡⟨n(0)|ρ̂^th_β|n(0)⟩ is the probability of finding the system in the n-th eigenstate at time t=0, and p_m → n(τ) is the transition probability, between the n-th and m-th eigenstates, due to the protocol. 
Fluctuations arising from the protocol and the measurements are encoded by p_n(0)p_m → n(τ) and constrained by the Jarzynski's equality <cit.> 

    ⟨e^-β w|=⟩ e^-βΔℱ,

where  Δℱ= ℱ[{λ^f_i}]-ℱ[{λ^0_i}] is the free energy difference between the two equilibrium configurations, corresponding to the initial and final Hamiltonian.

By the Jensen's inequality, Eq. (<ref>) implies that ⟨w|≥⟩Δℱ, which reflects the second law of thermodynamics. 
This leads to the definition of the average irreversible work <cit.>

    ⟨w_irr|:⟩= ⟨w|-⟩Δℱ,

which is directly related to the average  irreversible entropy production <cit.>

    ⟨𝒮_irr|:⟩=β⟨w_irr|.⟩

Both ⟨w_irr|$⟩ and⟨𝒮_irr|$⟩ give a measure of the irreversibility introduced by performing the unitary transformation  ρ̂(τ)=𝒰̂(τ,0)ρ̂^th_β[{λ^0_i}]𝒰̂^†(τ,0).
Strictly speaking, due to the unitary nature of the time evolution, no von Neumann entropy is generated during this process, with the entropy of the system remaining constant:

    𝒮(ρ̂(τ))=-{ρ̂(τ)logρ̂(τ)}≡𝒮(ρ̂^th_β[{λ^0_i}]).

Equation (<ref>) is referred to as a measure of irreversibility because, when the system is returned to the bath after the protocol, it relaxes from the out-of-equilibrium state to the thermal state ρ̂^th_β[{λ^f_i}]. 
This relaxation is a non-unitary process, and the entropy produced during this process is precisely ⟨𝒮_irr|$⟩.


We emphasize that, like the average work,⟨𝒮_irr|$⟩ is also the first moment of a probability distribution obtained within the two-point measurement framework.
Indeed, it is possible to define a stochastic variable,  associated with the production of irreversible entropy, as follows:  

    s_mn:=β(ϵ_m(τ)-ϵ_n(0))-βΔℱ.

Then, the probability distribution for entropy production (PDE), analogous to the PDW, takes the form

    P_τ(s)=∑_nmp_n(0)p_n → m(τ)δ(s-s_mn).


A complementary approach to investigate the statistical properties of work processes and irreversible entropy production is based on the Fourier transforms of the corresponding probability distributions, namely, P_τ(w) from Eq. (<ref>) and P_τ(s) from Eq. (<ref>). 
It is therefore convenient to rely on the characteristic function of work <cit.>

    χ_ν(w,τ) :=   ∫ dw e^iν wP_τ(w)   
    
     =   {
    e^iνℋ̂[{λ^f_i}]𝒰̂(τ,0)
    e^-iνℋ̂[{λ^0_i}]
                 ×ρ̂^th_β([{λ^0_i}])
    𝒰̂^†(τ,0)},

with associated moments 

    ⟨w^n(τ)|=⟩ (-i)^n∂^n_νχ_ν(w,τ)|_ν=0.
 
It is further instructive to introduce the characteristic function for irreversible entropy production:

    χ_μ(s,τ) :=   ∫ ds e^iμ s P_τ(s) 
    
    =    e^-i βμΔℱ{
    e^iβμℋ̂[{λ^f_i}]𝒰̂(τ,0) 
           ×
    e^-iβμℋ̂[{λ^0_i}]ρ̂^th_β[{λ^0_i}]
    𝒰̂^†(τ,0)
    },

with associated moments

    ⟨s^n(τ)|=⟩ (-i)^n∂^n_μχ_s(w,τ)|_μ=0.
 
The two quantities expressed in Eqs. (<ref>) and (<ref>) play a crucial role in the development of the thermal density functional framework for specific protocols, as detailed in Sec. <ref> (sudden quench) and appendix <ref> (finite-time protocols).



 §.§ The sudden quench protocol

A sudden quench involves an instantaneous shift of the work parameters, from {λ^0_i} to {λ^f_i}.

This variation occurs in an infinitesimally short timeframe, unlike the finite-time protocols covered in appendix <ref>. 



It can be demonstrated that all of the moments of the PDW in this scenario are functionals of the initial thermal densities.





This outcome is due to the fact that, in a sudden quench, the time  evolution operator approaches the identity operator, as the quench duration becomes infinitesimally small <cit.>:  lim_τ→ 0^+𝒰̂(τ,0)=1̂. 
Therefore, the characteristic function of work (<ref>) simplifies to an ensemble average over the initial Gibbs state:

    χ_ν(w,0^+) = {
    e^iνℋ̂[{λ^f_i}]
    e^-iνℋ̂[{λ^0_i}]ρ̂^th_β[{λ^0_i}]}.





For a Hamiltonian of the form (<ref>), any thermal average over the initial state can be expressed as a functional of the initial mean values {a^β 0_i}, as dictated by the generalized HK theorem outlined in Sec. <ref> and appendix <ref>.
Consequently, we have: χ_ν(w,0^+)=χ_ν(w,0^+)[{a^β 0_i}], which implies that all the moments of P_0^+(w) 
are functionals of {a^β 0_i}, with parametric dependence on both {λ^0_i} and {λ^f_i}.
More explicitly, by using Eq. (<ref>), these moments can be expressed as the following thermal averages:

    ⟨w^n|=⟩{(ℋ̂[{λ^f_i}]-ℋ̂[{λ^0_i}])^nρ̂^th_β[{λ^0_i}]},

which are functional of the initial thermal densities: ⟨w^n|=⟩⟨w^n|[⟩{a^β 0_i}].
Now, the linearity of the Hamiltonian ℋ[{λ_i}] in the work parameters {λ_i} leads to 

    ℋ̂[{λ^f_i}]-ℋ̂[{λ^0_i}]=∑_i(λ^f_i-λ^0_i)∂ℋ̂/∂λ^0_i.

Then, using Eqs. (<ref>) and (<ref>), 
the average work becomes:

    ⟨w|=⟩∑_i(λ^f_i-λ^0_i)a^β 0_i


Similarly, the characteristic function of the irreversible entropy production (<ref>),  for a sudden quench protocol, takes the simplified expression

    χ_μ(s,0^+)=   e^-i βμΔℱ{
    e^iβμℋ̂[{λ^f_i}]
               × e^-iβμℋ̂[{λ^0_i}]ρ̂^th_β[{λ^0_i}]
    }.

Again, by virtue of the thermal HK theorem, the final and initial free energies in Δℱ are functionals of {a^β f_i} and {a^β 0_i}, respectively.
On the other hand, as seen in Eq. (<ref>), the trace in Eq. (<ref>) is a functional of {a^β 0_i} only. 
Consequently, we can assert that P_0^+(s), or χ_μ(s,0^+), and the associated moments ⟨s^n|$⟩ are functionals of both{a^β f_i}and{a^β 0_i}.
In particular, the average irreversible entropy production takes the form:

    ⟨𝒮_irr|=⟩   β∑_i(λ^f_i-λ^0_i)a^β 0_i
           -β{ℱ[{a^β f_i}]-ℱ[{a^β 0_i}]}


We focus on sudden quenches of infinitesimal variation, where the work parameters{λ^0_i}change by an elementary amount to{λ^0_i + δλ_i}. 
In this context, we can derive an explicit functional form for the average irreversible entropy production.
In particular, we can expand Eq. (<ref>) in a Taylor series and apply the HF theorem, as expressed by Eq. (<ref>). 
This yields:

    ⟨𝒮_irr|=⟩ -β/2∑_i,jδλ_iδλ_j∂ a^β 0_i/∂λ^0_j.

We emphasize that Eqs. (<ref>) and (<ref>) demonstrate that the mean values of work and irreversible entropy production are explicit functionals of the initial thermal densities.  
As discussed in the following sections, this is particularly important for extracting information about many-body systems using the KS mapping, which enables the computation of thermal electron densities within a formally non-interacting framework.



 §.§ Fluctuation-dissipation relations in the sudden quench limit

We now recall that classical quasi-adiabatic processes follow the fluctuation-dissipation relation (FDR) 

    ⟨𝒮_irr|=⟩β^2/2σ^2_w,

whereσ^2_w=⟨w^2|-⟩⟨w|^⟩2represents the variance in the PDW <cit.>.
Recently, it has been demonstrated that for slow processes in open quantum systems, close to equilibrium, the FDR is given by Eq. (<ref>) minus a positive, purely quantum  term, which arises from the non-commutativity of the thermodynamic protocol <cit.>. 
Here, with the aim of obtaining an explicit functional form in terms of initial thermal densities for the second moment of the work probability distribution, we reobtain a similar generalized FDR  that holds in the infinitesimal sudden quench regime.
This should not be surprising, as an adiabatic process, i.e., one that is close to equilibrium throughout, can be considered as a sequence of a large number of sudden quenches, each followed by thermalization towards the equilibrium state <cit.>.


To this end, we focus on the second moment of the PDW. 
Then, we distinguish the case where the final and initial Hamiltonians share a common eigenbasis, and the case where they do not.
Additional details on the following derivations are provided in appendix <ref>.
In the specific scenario where[ℋ̂[{λ^f_i}],ℋ̂[{λ^0_i}]]=0, the second moment of the PDW is given by:

    ⟨w^2|_⟩c=∑_i,jδλ_iδλ_j a^β 0_ia^β 0_j-1/β∑_i,jδλ_iδλ_j∂ a^β 0_i /∂λ^0_j.

Notably, Eq. (<ref>) expresses an explicit functional of the initial mean values, independently of the amplitude of the sudden  quench.
Nonetheless, with an infinitesimal sudden quench, we can utilize the expressions for the average work, Eq. (<ref>), and the average irreversible entropy production, Eq. (<ref>), to rewrite Eq. (<ref>) as: 

    ⟨w^2|_⟩c=⟨w|^⟩2+2/β^2⟨𝒮_irr|,⟩

This relation validates the FDR, in its classical form, as given by Eq. (<ref>), 

to the leading order in{δλ_i}.

Turning to the instance where the initial and final Hamiltonians do not commute, and using Eq. (<ref>), the second moment of the PDW is still a functional of the initial equilibrium thermal densities.
Specifically, the latter can be split into the following two parts: 

    ⟨w^2|=⟩⟨w^2|_⟩c+ Θ_2[{a^β 0_i}],

whereΘ_2[{a^β 0_i}]arises directly from the incompatibility of the two Hamiltonians. 

A possible approximation method for this functional is provided in sec <ref>.
Operating again in the infinitesimal sudden quench limit, we can plug Eqs. (<ref>) and (<ref>) into Eq. (<ref>).
By doing so, we recover the generalized FDR,

    ⟨𝒮_irr|=⟩β^2/2σ^2_w -β^2/2Θ_2[{a^β 0_i}],

which takes into account both thermal fluctuations and quantum fluctuations due to[ℋ̂[{λ^f_i}],ℋ̂[{λ^0_i}]]≠0.



§ THERMAL KOHN-SHAM MAPPING FOR QUANTUM THERMODYNAMICS
  
The results discussed so far in previous sections are formally exact, at least in the limiting conditions of the protocols investigated for the evolution of the coupling parameters.
However, as a many-body system grows in complexity, the number of interactions and possible configurations needed to determine the exact thermal density becomes computationally infeasible.  
To address this challenge, the KS scheme <cit.> provides a powerful approach within the framework of DFT for developing efficient approximations. 
This method relies on defining a formally non-interacting many-body system, the KS system, which is designed to replicate the same particle density as the original interacting physical system.

For systems governed by the Hamiltonian (<ref>), and building on methods developed in earlier studies <cit.> at zero temperature, the KS approach can be applied as follows.
We assume the existence of a set of auxiliary systems, each described by the Hamiltonian 

    ℋ̂^ks=ℋ̂^ks_0+ ∑^ℒ_i=1λ^ks_i Â_i,

where the one-body operatorℋ̂^ks_0replaces the complex many-body termℋ̂_0in Eq. (<ref>).
The KS Hamiltonian simplifies the problem by focusing on non-interacting particles in an effective potential associated to specific coupling parameters.
We further assume thatℋ̂^ksyields the same set of thermal densities as the original Hamiltonian:

    {ρ̂^th_β[{λ_i}]Â_i}={ρ̂^th_βks[{λ^ks_i}]Â_i}.

In this setting,ρ̂^th_βandρ̂^th_β, ksdenote the thermal density matrices of the original and KS systems, respectively, both parameterized by inverse temperatureβand coupling parameters{λ_i}and{λ^ks_i}.

The MHK theorem clearly holds for the KS Hamiltonian, though with some restrictions at absolute zero temperature <cit.>. 
Consequently, the coupling parametersλ^ks_iare functionals of the thermal averages{a^β_i}, i.e.,λ^ks_i = λ^ks_i[{a^β_i}]. 
At this point, the following MKS equations can be solved self-consistently for{a^β_i}:

    ℋ̂^ks_0+∑^ℒ_i=1(λ^h-xc_i[{a^β_i}]+λ_i)Â_i|ϕ^i_β⟩=ϵ^i_β|ϕ^i_β⟩
       a^β_i={ρ̂^th_βks[{λ^h-xc_i[{a^β_i}]+λ_i}]Â_i}.

Here, the effective parametersλ^h-xc_i[{a^β_i}]=λ^ks_i-λ_iplay the role of the Hartree (H) and exchange-correlation (XC) potentials in the usual KS mapping,  which account for the effect of the many-body interaction term inℋ̂_0. 
It is worth recalling that while the eigensystem of the non-interacting KS Hamiltonian exactly reproduces the thermal density, it generally do not correspond to the eigensystem of the interacting Hamiltonian <cit.>.

The approach outlined here is particularly useful when the one-body Hamiltonianℋ̂^ks_0has a simple form, such as in the case of a chain of interacting fermions, whereℋ̂^ks_0reduces to a kinetic energy operator.
In these scenarios, as is typically done within the KS framework, suitable approximations can be employed for the functionalsλ^h-xc_i[{a^β_i}].
In other terms, any thermodynamic quantity expressed as an explicit functional of the thermal densities can be evaluated through a KS mapping, with an accuracy dictated by the approximations made for the functionalsλ^h-xc_i[{a^β_i}].
Nonetheless, not all quantities in the MKS equations can be directly expressed as functionals of the densities. 
For example, the functional form ofΘ_2[{a^β 0_i}]in Eq. (<ref>) requires reasonable approximations to be determined.



 §.§ Local density approximation for Θ_2[{a^β 0_i}]

The local density approximation (LDA) is the simplest and most widely used approach for modeling XC effects in DFT. 
For instance, the LDA has been effectively employed  to develop functionals for calculating the entanglement in spatially inhomogeneous many-fermion systems <cit.>.
To construct an LDA scheme for an inhomogeneous system, it is necessary to have an analytical solution for the corresponding homogeneous problem, where all coupling parameters are equal, i.e.,λ_i = λ.  
In the homogeneous case, the functionalΘ_2[{a^β 0_i}]reduces toΘ^hom_2[a^β 0], where:

    a^β 0=1/ℒ{ρ̂^th_β[{λ_i}]∑^ℒ_i=1Â_i}.

Based on this, the following LDA scheme can be put forward:

    Θ^lda_2[{a^β 0_i}]=∑_iΘ^hom_2[a^β 0]|_a^β 0→ a^β 0_i.
  
A crucial aspect of this implementation is that Eq. (<ref>) approximates the fluctuations inΘ_2[{a^β,0_i}]arising from the incompatibility between the pre- and post-quench Hamiltonians, as discussed in Sec. <ref>. 
Therefore, it is essential that the homogeneous system satisfies the condition:[ℋ̂_0, ∑_i=1^N Â_i] ≠0. 
Otherwise,Θ^hom_2[a^β 0]would be identically zero.



 §.§ Approximation for Θ_2[{a^β 0_i}] via perturbation treatment of the KS system

The fictitious KS world is governed by a Hamiltonian that differs from the one describing the actual system under investigation.
However, if the exact form of the functionalλ^h-xc_i[{a^β_i}]is known, the KS framework can accurately reproduce the density of the original interacting system.  
Therefore, a well-founded idea is to treat the KS Hamiltonian as a zeroth-order approximation to the `true' Hamiltonian <cit.>. 
This is because the KS Hamiltonian encapsulates key aspects of the many-body properties inherent the real system.
In this perturbation-like approach, we can express the original Hamiltonian asℋ̂ = ℋ̂^ks + Δℋ̂, whereΔℋ̂ = ℋ̂ - ℋ̂^ks.
Accordingly, any observable property𝒬[{a^β 0_i}]can be expanded as

    𝒬[{a^β 0_i}] = 𝒬^ks[{a^β 0_i}] + Δ𝒬[{a^β 0_i}],

where𝒬^ks[{a^β 0_i}]represents the zeroth-order approximation of the quantity of interest. 
This method is particularly advantageous when the LDA scheme is not applicable. 
For example, if the underlying homogeneous system satisfies[ℋ̂_0,∑^N_i=1Â_i]= 0, the zeroth-order approximation can effectively address the incompatibility between the initial and final forms of the interacting Hamiltonian.
In such cases,Θ_2[{a^β 0_i}]can be approximated using the corresponding quantity calculated within the KS framework.



§ THE KOHN-SHAM SCHEME IN THE CANONICAL ENSEMBLE

In Sec. <ref>, we introduced a finite-temperature KS mapping, which enables the accurate evaluation of equilibrium thermal densities through iterative solutions of the self-consistent MKS equations, i.e., Eqs.(<ref>) and (<ref>). 


However, dealing with statistical systems that have a fixed number of indistinguishable particles, even within the non-interacting KS framework, presents considerable challenges <cit.>, making the solution of the MKS equations computationally demanding.
This complexity is one reason why finite-temperature DFT calculations typically employ the grand canonical ensemble, where the average number of particles is fixed by⟨N|=⟩∑_i f_μ(ϵ_i). 
In a many-fermion system, the Fermi-Dirac distributionf_μ(ϵ_i)=(1+e^β(ϵ_i-μ))^-1describes the occupation probabilities of the KS eigenstates, withϵ_irepresenting the KS eigenvalues andμthe chemical potential.  
This framework offers the advantage of a straightforward expression for thermal electron densities:n^β(𝐫)=∑_i f_μ(ϵ_i)|ψ^ks_i(𝐫)|^2, based on the KS wavefunctionsψ^ks_i(𝐫)<cit.>. 
Here, we present a method to compute thermal densities within the canonical ensemble using the KS mapping while keeping the calculations feasible.  

To elucidate our approach, we consider a system ofNinteracting fermions on a lattice, characterized by the Hamiltonian (<ref>), with the specific form

    ℋ̂=ℋ̂_0+∑^ℒ_i=1 V_i n̂_i.

The universal part of this Hamiltonian reads

    ℋ̂_0=𝒯̂+𝒲̂,

where

    𝒯̂=-J
    ∑_i; σ=↑,↓(ĉ^†_i,σĉ_i+1,σ + h.c.)

is the kinetic term, while𝒲̂accounts for the two-body interaction.
In Eq. (<ref>),ĉ^†_i,σandĉ_i,σdenote the creation and annihilation operators for a fermion with spinσ=↑,↓, andn̂_i=n̂_i,↑+ n̂_i,↓is the total number operator for thei-th site.

Given the form of the external potential, the thermal densities are naturally defined asn^β_i = {ρ̂^th_β[{V_i}]n̂_i}. 
According to the MHK theorem, the following correspondence holds:

    {V_i}{n^β_i}ρ̂^th_β≡ρ̂^th_β[{n^β_i}].

As remarked in Sec. <ref>, there exists a KS system having exactly the same set of thermal densities as the original interacting system. 
The corresponding KS Hamiltonian is

    ℋ̂^ks=𝒯̂+∑^ℒ_i=1(V^h-xc_i[{n^β_i}]+V_i)n̂_i.

At this point, we seek a more computationally tractable equation for the thermal densities than Eq. (<ref>). 
In particular, we build on previous research <cit.> that examined the canonical partition function forNnon-interacting fermions.
These studies enable us to express the canonical partition function for the KS fermions using the following recursive formula: 

    𝒵_N = 1/N∑^N_m=0(-1)^m-1𝒵_1(mβ)𝒵_N-m(β),

where𝒵_1(0)=1and𝒵_1(mβ)=∑_i e^-mβϵ_iform>1.

Given the structure of the KS Hamiltonian (<ref>), the partition function for a KS system withN_↑spin-up andN_↓spin-down fermions becomes:𝒵^ks_N=𝒵^ks_N_↑𝒵^ks_N_↓, 
with particle number conservation ensured byN = N_↑+ N_↓.
The corresponding equilibrium free energy is then:ℱ^ks=-1/βlog(𝒵^ks_N), 
from which the equilibrium thermal densities can be extracted using the HF theorem as

    n_i^β= ∂ℱ^ks/∂ V^ks_i.

The combination of Eq. (<ref>) and Eq. (<ref>) forms the self-consistent foundation of our finite-temperature KS approach, which, in principle, exactly reproduces the thermal densities of the original interacting system.



§ A NOTABLE EXAMPLE

This section is dedicated to validating our finite-temperature KS approach, as defined by Eqs. (<ref>) and (<ref>), within the context of the Hubbard model <cit.>.
Specifically, we examine a scenario where electrons are influenced by an inhomogeneous external potential dependent on the parameterv_0.
The system Hamiltonian is expressed as: 

    ℋ̂=ℋ̂_0+∑^ℒ_i=1 V_i(v_0) n̂_i,

whereV_i(v_0) = f_i v_0with the dimensionless set{f_i}that defines the spatial shape of the external potential along the chain. Hereℋ̂_0takes the form given in Eq. (<ref>), with the two-body interaction term given by:  

    𝒲̂=∑^ℒ_i=1 Un̂_i,↑n̂_i,↓.

As described in Sec. <ref> the system is  initially prepared in the Gibbs stateρ̂^th_β[{n^β_i}].  
Subsequently, it is decoupled from the thermal bath and undergoes an instantaneous quench in the work parameter, with amplitudesδV_i = f_i δv_0. 


We focus on the half-filled Hubbard model, with the total spin along thezdirection set to zero. 
Our analysis considers two external potentials. One that decreases linearly along the chain, 

    𝒱̂_ext = ∑^ℒ_i=1[v_0 - 2v_0(i-1)/ℒ-1] n̂_i,

and another one with a harmonic dependence,  

    𝒱̂_ext = ∑^ℒ_i=11/2v_0[i - ℒ+1/2]^2n̂_i.
 
We begin by analyzing the exact results for the Hubbard dimer, as presented in Sec. <ref>. We then compare these results with those obtained from our KS  mapping, as detailed in Sec. <ref> and further explored in Sec. <ref>.
Next, we extend our study to longer Hubbard chains probed by the linear potential defined in Eq. (<ref>), as outlined in Sec. <ref> and Sec. <ref>. 
In particular, in Sec. <ref>, we compare exact results for systems with up to8sites with corresponding ones from our KS mapping. 
Finally, in Secs. <ref> and <ref>, we investigate the impact of electron-electron interactions on work extraction in longer chains, considering both linear and harmonic potentials given in Eqs. (<ref>) and (<ref>).



 §.§ Exact results for the Hubbard dimer

In the two-particle subspace with total spin zero along thez-axis, the Hamiltonian (<ref>) characterizes a two-site Hubbard chain and is represented by the matrix

    ℋ̂=̇[  U+2V_1      -J       J       0;      -J V_1+V_2       0      -J;       J       0 V_1+V_2       J;       0      -J       J  U+2V_2 ],

in the basis{|↑↓, 0⟩, |↑, ↓⟩, |↓, ↑⟩, |0, ↑↓⟩}.
This straightforward, exactly solvable model exhibits a diverse range of physical phenomena <cit.>, including a precursor to the Mott metal-insulator transition and, influenced by the external potential, a precursor to the ionic insulator transition.
The two transitions are in competition, with the former favoring single-site occupation and the latter promoting double-site occupation.
The metal phase emerges in the narrow region whereU ∼2v_0, driven by the interplay of the interaction term and the external potential.

In Fig. <ref>(a), we examine the average extracted work,⟨w|_⟩ex = -⟨w|$⟩, as a function of U and v_0, following a sudden quench of amplitude δ v = 0.05J. 
For U > 2v_0, the system enters the Mott insulating phase, leading to a decrease in the extractable work as the interaction strength U increases.

In this phase, double occupancy of sites becomes energetically unfavorable, rendering work extraction through the external potential quench impractical. 

This behavior can be understood by examining the thermal densities in the limit of large U: for U ≫ 2v_0, we have n^β_1 ∼ n^β_2 ∼ 1, which, according to Eq.(<ref>), results in ⟨w|∼⟩0.

Conversely, for U < 2v_0, the system is in the ionic insulating phase, where the extractable work increases as U decreases, reaching its maximum value as U approaches zero.


This is because for U≪2v_0 we find  n^β_1∼ 2 and n^β_2∼ 0, leading to ⟨w|∼⟩-2δv_0 as per Eq. (<ref>).


In Fig. <ref>(b), we show the average  irreversible entropy production ⟨𝒮_irr|$⟩ for the same process.

As expected,⟨𝒮_irr|$⟩  exhibits a pronounced peak in the metallic region separating the Mott insulating phase from the ionic insulating phase. 



This behavior can be understood through the dependence of the average irreversible entropy production on the thermal densities, as described in Eq. (<ref>). 
The thermal densities are sensitive to small variations in the work parameter v_0 when U ∼ 2v_0, which is also reflected in the peak observed in the derivatives of the thermal density with respect to v_0, shown in Fig. <ref>(c).
