Functional theories of thermoelectric phenomena

# Functional theories of thermoelectric phenomena

F. G. Eich Max-Planck-Institute for the Structure and Dynamics of Matter, Luruper Chaussee 149, D-22761 Hamburg, Germany Department of Physics, University of Missouri-Columbia, Columbia, Missouri 65211, USA    M. Di Ventra University of California at San Diego, La Jolla, California 92093, USA    G. Vignale Department of Physics, University of Missouri-Columbia, Columbia, Missouri 65211, USA
July 9, 2019
###### Abstract

We review the progress that has been recently made in the application of time-dependent density functional theory to thermoelectric phenomena. As the field is very young, we emphasize open problems and fundamental issues. We begin by introducing the formal structure of thermal density functional theory, a density functional theory with two basic variables – the density and the energy density – and two conjugate fields – the ordinary scalar potential and Luttinger’s thermomechanical potential. The static version of this theory is contrasted with the familiar finite-temperature density functional theory, in which only the density is a variable. We then proceed to constructing the full time-dependent non equilibrium theory, including the practically important Kohn-Sham equations that go with it. The theory is shown to recover standard results of the Landauer theory for thermal transport in the steady state, while showing greater flexibility by allowing a description of fast thermal response, temperature oscillations and related phenomena. Several results are presented here for the first time, i.e., the proof of invertibility of the thermal response function in the linear regime, the full expression of the thermal currents in the presence of Luttinger’s thermomechanical potential, an explicit prescription for the evaluation of the Kohn-Sham potentials in the adiabatic local density approximation, a detailed discussion of the leading dissipative corrections to the adiabatic local density approximation and the thermal corrections to the resistivity that follow from it.

## I Introduction

Density functional methods occupy a central position in the landscape of modern computational materials theory Dreizler and Gross (1990); Giuliani and Vignale (2005); Marques et al. (2012); Ullrich (2012). Not only do they provide an indispensable tool for the calculation of equilibrium properties of materials, but they are also being widely used for the calculation of transport properties at all length scales Ullrich (2012); Di Ventra (2008). Key to the success of these methods is their ability to incorporate important many-body effects, arising from electron-electron interaction, in an intuitive and computationally affordable scheme. In spite of all this, little attention has been paid so far to the possibility of applying density functional methods directly to the study of thermoelectric phenomena. These are transport phenomena, such as the Seebeck and the Peltier effect, in which electronic degrees of freedom are involved in an essential manner, with important contributions from electron-phonon interactions Nolas et al. (2001); Dubi and Di Ventra (2011). Therefore, they seem natural candidates for a time-dependent non equilibrium density functional treatment. Beyond its immediate practical value, such a treatment would forge a link between density functional theory – a formally exact theory valid at all length scales – and classical hydrodynamics, a theory of locally conserved particle, momentum and energy energy density, which is valid at long wavelengths and long time scales Landau and Lifshitz (1959); Andreev et al. (2011). Such a link is not complete at present, due to the absence of thermal and thermoelectric couplings in density functional theory.

In this article, we review the progress that has recently been made in the application of non equilibrium density functional theory to thermoelectric phenomena Eich et al. (2014a). Since this line of research is a very recent development, we focus specifically on open problems and fundamental issues. We start by introducing the basic variable of the time-dependent thermal density functional theory, namely, the energy density, and discuss the problems associated with its non-unique definition. Next, following a classic paper by Luttinger Luttinger (1964), we introduce the mechanical field conjugate to the energy density, which we call “thermomechanical potential”. This thermomechanical potential serves as mechanical proxy for local temperature variations Shastry (2009). We discuss how the results of the standard Landauer theory Landauer (1957, 1989), for mesoscopic thermal transport can be reproduced in the framework of the thermomechanical potential approach, while the latter exhibits higher flexibility in the nonlinear regime Eich et al. (2014b). We conclude by introducing the generalization of the thermomechanical potential to the “thermal vector potential” that couples to thermal currents Tatara (2015).

We then proceed to develop the formalism of the new functional theory from a generalized stationary action principle and introduce the Kohn-Sham reference system and the formal construction of its one-particle effective potentials Eich et al. (2014a). This formalism is suitable for the calculation of electrical and thermal currents, as well as particle and energy densities. At this point we present a comparison between the thermal density functional theory and alternative approaches that can potentially be applied to the same class of problems, e.g., the stochastic density functional theory for open systems D’Agosta and Di Ventra (2007); Biele and D’Agosta (2012); D’Agosta and Di Ventra (2013), and the time-dependent thermal transport theory Biele et al. (2015), and discuss ways in which these different theories could complement each other.

The last part of the review is concerned with the practical issue of constructing an explicit functional for thermal density functional theory: we examine the adiabatic local density approximation and the leading dissipative corrections to it in the light of recent numerical calculations of the free energy density Brown et al. (2013a); Karasiev et al. (2014) and thermal transport coefficient of the homogeneous electron gas at finite temperatures. Lastly, we present a few applications of the theory which are relevant to experiments such as the existence of thermal dynamical corrections to the electrical resistance of conductors Sai et al. (2005); Koentopp et al. (2006); Vignale and Di Ventra (2009); Kurth and Stefanucci (2013), current- and thermal-current-induced quantum oscillations in the local temperature Bergfield et al. (2015); Eich et al. (2016), transient flow of charge and energy and dynamical temperature waves in mesoscopic systems Eich et al. (2016).

## Ii Luttinger’s thermomechanical potential

Hydrodynamics is one of the most versatile and successful methods ever introduced to describe the dynamics of strongly interacting many-particle systems. In its full fledged form the hydrodynamic theory consists of three equations Huang (1987); Puff and Gillis (1968): the continuity equation, which connects the time derivative of the particle density, , to the divergence of the particle current density (proportional to the momentum density); the Euler equation, which expresses the time-derivative of the particle current density in terms of the divergence of an internal stress tensor , plus external volume forces; and the heat transport equation, which expresses the time derivative of the energy density, in terms of the divergence of the energy current density . Essential to the closure of the scheme is the existence of approximate linear relations that allow us to express the stress tensor (alias momentum current) and the heat current, , in terms of the gradients of local velocity and temperature fields, which in turn are determined by the particle current density and energy density. There is a very good reason why hydrodynamics employs the densities of particles, momentum, and energy, and no more. In the limit of slow spatial variation these are quasi-conserved quantities and therefore their time evolution is very slow. Provided the external potentials are slowly varying on the microscopic scale, the evolution of the quasi conserved quantities is sufficient to characterize the dynamics of the system, especially so when the particle-particle interactions are strong and frequent, resulting in small deviations from local thermal equilibrium.

Electrons in solid state devices rarely satisfy the conditions for the applicability of hydrodynamics: the external potentials arising from nuclei, impurities, lattice vibrations, are usually not slowly varying on the microscopic scale. Further, the description of the particles is necessarily quantum mechanical. Nevertheless, a powerful many-body formalism known as time-dependent density functional theory (TDDFT) Runge and Gross (1984); Tokatly (2005a, b); Ullrich (2012) allows us to derive closed equations of motion for the particle density and the current density, similar to the hydrodynamic equations, but valid, in principle, at all length scales. The foundation of the theory is the Runge-Gross theorem Runge and Gross (1984), according to which the quantum stress tensor is uniquely determined by the density and the initial state of the system. In practice, however, the dependence of the stress tensor on the underlying densities is strongly nonlocal–both, in space and time–and needs to be a subject of drastic approximations. Further, the Fermi statistics of the electrons forces in many cases an orbital representation of the densities, the so-called Kohn-Sham representation Kohn and Sham (1965), which tends to obscure the hydrodynamic structure of the theory. In spite of these limitations, TDDFT and its companion time-dependent current-density functional theory (TDCDFT) Vignale (2004) have grown to be widely used approximation methods in the study of electronic dynamics, especially optical excitations and electric transport, at the nanoscale Marques et al. (2012).

The topic of this review is the progress that has recently been made in extending the TDDFT so that it becomes capable to deal with thermoelectric phenomena, such as heat transport, local heating in current-carrying systems, temperature-voltage conversion in nanoscale electronic devices. Our focus is on phenomena that depend essentially on electronic degrees of freedom. There is, of course, significant heat transport in electrically insulating materials, where the excitations responsible for the transport are phonons or magnons. As a first step we concentrate here on the problem of calculating heat currents in electric, non magnetic conductors, where heat is primarily carried by electrons in the vicinity of the Fermi surface. Parts of our system may be in contact with external reservoirs, which enforce a local temperature. Other parts may be allowed to “float” so as to adjust their temperature – assuming that such a notion still makes sense – to the value that is most compatible with the microscopic state of the system. In Sec. IV.2 we use these floating probe leads to compute the local temperature profile in a conducting nano wire.

### ii.1 Thermal energy density

A major difficulty in carrying out such a program is that quantities such as temperature and entropy have a statistical significance, but not a clear mechanical one. TDDFT, on the other hand is designed as a quantum mechanical theory, in which the quantum states evolve deterministically starting from some initial equilibrium ensemble, according to the equations of quantum mechanics (for a stochastic approach, in which the electrons are allowed to interact with a “thermal bath”, see section IV.1). The unperturbed Hamiltonian for such a system is

 ^H0=^T0+^V0+^U0 (1)

where

 ^T0=∫d3rℏ22m[∇^Φ†(r)]⋅[∇^Φ(r)] (2)

is the kinetic energy,

 ^V0=∫d3r v(r)Φ†(r)^Φ(r) (3)

is the external (static) potential energy, and

 ^U0=12∫d3r∫d3s^Φ†(r)^Φ†(s)w(|r−s|)^Φ†(s)^Φ(r) (4)

is the electron-electron interaction energy. The basic variables of standard TD(C)DFT are the particle density

 ^n(r)=^Φ†(r)^Φ(r) (5)

and the particle current density

 (6)

The corresponding expectation values are denoted by the same symbols without the hat, i.e., and . All these quantities have a clear mechanical significance and are unambiguous insofar as they are the sources of a measurable electromagnetic field.

To introduce thermal effects, the natural local variable would be the entropy density. According to equilibrium statistical mechanics the entropy is the expectation value

 S0=β⟨^H0−μ^N−Ω0⟩ (7)

where is the inverse equilibrium temperature, is the chemical potential, is the total number operator, and is the grand thermodynamic potential defined as

 e−βΩ0=Tr e−β(^H0−μ^N). (8)

The equilibrium expectation value of any operator is

 ⟨^A⟩=Tr[^Ae−β(^H0−μ^N−Ω0)] . (9)

Since is just a number, which commutes with the Hamiltonian, Eq. (7) strongly suggests that the operator

 ^Q0≡^H0−μ^N0 (10)

is the mechanical operator that corresponds to the entropy (times the temperature). In fact, this is not quite true, due to the presence of the term in Eq. (7). Nevertheless, the role played by this quantity in the theory of thermal transport is so central that we will from now on refer to as the “thermal Hamiltonian” and we introduce the associated thermal energy density operator , in the following manner

 ^q0(r) =^h0(r)−μ^n(r) , (11a) ^h0(r) =^t0(r)+^v0(r)+^u0(r) , (11b) ^t0(r) =ℏ22m[∇^Φ†(r)]⋅[∇^Φ(r)] , (11c) ^v(r) =v(r)^n(r) , (11d) ^u0(r) (11e)

It is easy to verify that and .

Before proceeding, we must observe at this point that the thermal energy density (or the energy density, for that matter) is not unique. For example, focusing on the kinetic energy density, the expression

 ^t1(r)=−ℏ24m(^Φ†(r)[∇2^Φ(r)]+[∇2^Φ†(r)]^Φ(r)) , (12)

is a priori as good as the expression

 ^t2(r)=ℏ22m[∇^Φ†(r)]⋅[∇^Φ(r)] . (13)

Both expressions integrate to the total kinetic energy , even though they are obviously different at each point. On the one hand, the form has the advantage that it can be derived from the non-relativistic limit of the time-component of the relativistic stress tensor for Dirac electrons, the latter being arguably a physical quantity, namely, the source of the gravitational field. On the other hand is strictly positive, which is intuitively expected of a kinetic energy. More importantly, the form will allow us to construct an energy current that has a simple scaling property Qin et al. (2011) upon introduction of the Luttinger potential, discussed below. In the following we will therefore adopt the form of the kinetic contribution to the thermal energy density.

### ii.2 Thermal potential

Following an original idea introduced by Luttinger, we will now consider the effect of an external field (for the time being, static) that couples to the thermal energy density. The perturbed thermal Hamiltonian is

 ^Q=∫d3r^q0(r)[1+ψ(r)]. (14)

From the form of the thermodynamic potential it should be evident that, as long as equilibrium conditions are maintained, a variation in the Luttinger potential corresponds to a local variation of the temperature:

 δψ(r)=δββ≃−δTT . (15)

In the next section we will show that, under such equilibrium conditions, the relation between and the local energy density is invertible so we can take as an indicator of the temperature that would correspond to a given local energy density if the system were in equilibrium.

However, the real purpose of Luttinger’s potential is to drive the system out of equilibrium, simulating the effect of a temperature gradient in driving a thermal current. To accomplish this we must allow to be time-dependent, i.e., we have and the Hamiltonian becomes time-dependent. It turns out that, under non equilibrium conditions, the time-dependent Luttinger potential acts as a mechanical proxy for a temperature gradient, i.e., . Notice the change in sign compared to the equilibrium expression (15). The identification of with under non-equilibrium conditions is simple, but surprisingly subtle. The point is that under equilibrium conditions we have both a non-uniform and a non-uniform temperature, which is associated with the nonuniform thermal energy density , and the thermal current is zero. At the same time, we can think of the equilibrium state as the net result of two non-equilibrium processes, which exactly cancel each other. More precisely, the thermal current driven by must be the opposite of the thermal current driven by , so that the sum of the two currents is zero (notice that this is essentially the argument leading to Einstein’s relation between diffusion constant and conductivity). But we have already seen that is the opposite of : it follows that the thermal conductivity, which relates the thermal current to , is exactly the same as the thermal conductivity that relates the thermal current to . In this precise sense acts as a mechanical proxy of  Shastry (2009). We will often refer to as to the “thermomechanical potential” in what follows.

In practical applications we consider an electronic device that is connected by leads and thermal contacts to a set of thermal reservoirs, which are separately in equilibrium at chemical potentials and temperatures , as shown in Fig. 1. In analogy with the theory of mesoscopic electric transport, we will require some of the reservoirs to have fixed temperature and inject a thermal current (more precisely defined below) into the system, while other reservoirs do not exchange any current and their temperatures are allowed to “float” to a local equilibrium value. The first type of reservoir is a model for a thermal source, while the second type is a model for a local temperature probe Engquist and Anderson (1981). The standard treatment of this electro-thermal device is based on the Landauer-Büttiker (LB) multi terminal formalism, in which the electric and thermal current are expressed in terms of transmission probabilities of electrons from one reservoir into another.

To treat this system employing Luttinger’s thermomechanical potential we have two possibilities Eich et al. (2014b). The first possibility (method I), which closely mimics the LB approach, is to start with the reservoirs in equilibrium with static thermomechanical fields . It is clear that these fields produce in each reservoirs the same populations that would be produced by the temperatures in the LB approach. At the initial time the potentials are turned off in all the reservoirs and the subsequent evolution of the system is tracked. This evolution takes place under the action of the unperturbed Hamiltonian () just as in the LB model, but with initial populations that were dictated by the temperatures . In Ref. Eich et al. (2014b) we have shown that the long-term or steady-state currents that flow in the reservoirs are identical to those obtained in the LB approach.

The second possibility (method II) is to start with all the reservoirs at the same temperature. At the initial time different thermomechanical potentials are turned on in the reservoirs, and the subsequent evolution of the system is calculated. The physical situation is now different from the LB setup, because the initial populations reflect equal temperatures in all reservoirs, but the system evolves under the action of persistent thermomechanical potentials in the leads. Long-term or steady-state currents are established, in the attempt to equilibrate with the now constant value of . The results are different from those of the LB method, which, as just described, is mimicked by method I. However, in the linear regime both methods are identical provided we identify . The differences that appear under strong bias arise entirely from the persistent modification of the Hamiltonian via . It is precisely in this situation that one must ask the question: does the thermomechanical potential have a genuine physical meaning, or is it just a formal device to calculate thermal responses? Are there any physical situations in which modeling the system in terms of a persistent potential, i.e., method II, would be more appropriate than method I?

The answers to these questions are not entirely clear at this time. It seems plausible that the application of a persistent thermomechanical potential would be an appropriate description for an adiabatic heating process, in which the temperature is changed not by heat transfer, i.e., repopulation of energy levels, but by an actual mechanical modification of the energy levels.

### ii.3 Equations of motion

Let us now consider more closely the form of the thermal or heat current. This is identified with the entropy current density times the temperature. One way to proceed is to calculate the time derivative of the entropy density operator, which we take to be

 β−1^s(r)=^q(r)−ω(r). (16)

The local thermodynamic potential density is just a number, which does not contribute to the equation of motion. It follows that

 β−1∂t^s(r)=∂t^q(r)≡−∇⋅ȷq(r), (17)

where the last equality is our definition of the heat current, which is given by the combination of energy current density and particle current density:

 ȷq(r)≡ȷh(r)−μȷn(r) (18)

Since we are now in the presence of the thermomechanical potential both the energy and thermodynamic potential density are multiplied by the renormalization factor in the following manner

 q(r) ≡ q0(r)[1+ψ(r)], ω(r) ≡ ω0(r)[1+ψ(r)]. (19)

Also the particle current density is modified in a similar way as can be established by computing the continuity equation in the presence of the coupling to . Since only the commutator of the density and the kinetic energy contributes, and the thermomechanical potential enters as a “mass renormalization” in the kinetic energy, this leads to

 ^ȷn(r)=[1+ψ(r)] ^ȷn,0(r) , (20)

The expression for the thermal current density is more convoluted. Because this current contains both an energy and a velocity, one would hope for it to scale with the local factor . This can be achieved for the current of kinetic energy in the form, and for the external potential energy current, but not for the interaction energy density current, which has a more complicated nonlocal structure. The complete energy current density, , is given by

 ^ȷh(r)=^ȷt(r)+^ȷv(r)+^ȷu(r)+^ȷf(r). (21)

The first term on the right hand side,

 (22)

is the kinetic energy contribution. A detailed derivation of this result is presented in Appendix A.1. The potential energy current is given by

 ^ȷv(r)≡[1+ψ(r)]2v(r)^ȷn,0(r) , (23)

and is derived in Appendix A.2. The remaining terms arise from the electron-electron interaction, and their expressions, are given in Eqs. (143) and  (147) of Appendix A.3. In particular, is a “convective” current of interaction energy carried by a moving volume element of the electron fluid: this current satisfies the local scaling with (see Eq. (143)). The last term, , describes the work done on the volume element by the surrounding medium and does not satisfy the local scaling (see Eq. (147)).

Let us consider, for orientation, the case of a homogeneous interacting electron gas of uniform density , energy density , and thermodynamic potential density . In equilibrium all currents vanish. Let us now perform a Galilean transformation to a reference frame with velocity so that, in the new reference frame, all the electrons are imparted a positive velocity . In the new reference frame it is a simple exercise to show that the currents are

 ȷn=n0v (24)

and

 ȷh=(ε0+p0)v (25)

where is the pressure. Then we can immediately verify that

 ȷq=ȷh−μȷn=(ε0−μn0−ω0)v=β−1n0s0v, (26)

where is the entropy density (thus, is entropy per particle). This confirms our surmise that the current (here calculated in the absence of the thermomechanical potential ) is indeed the entropy current.

Up to this point we have closely followed Luttinger’s approach to mechanically simulate a thermal gradient. This is a sufficient basis for the time-dependent DFT that we construct in the next section. Before closing this section, however, we want to briefly comment on a different, and interesting approach, that has recently been proposed by Tatara Tatara (2015), motivated by calculations of the Nernst effect (Hall effect driven by a thermal gradient) Qin et al. (2011).

The basic idea is that the Luttinger interaction term

 ^Qψ=∫d3r ^q0(r)ψ(r,t) (27)

can be eliminated, to first order in , in favor of an interaction with a thermal vector potential :

 ^QA=∫d3r ^ȷq(r)⋅Ath(r,t). (28)

The relation between the thermal vector potential and the Luttinger potential is

 ∂tAth(r,t)=−∇ψ(r,t). (29)

The reader will recognize the similarity between this transformation and the familiar gauge transformation of electrodynamics, in which the scalar electric potential , coupling to the charge density, can be eliminated in favor of a vector potential , coupling to the electric current density, such that . However, the gauge transformation of electrodynamics is exact, while the present transformation is approximate and becomes exact only to linear order in the Luttinger potential. To see that this is the case, it is sufficient to apply to the Luttinger Hamiltonian of Eq. (14) the unitary transformation

 ^U(t)=e−i∫d3r ^q0(r)χ(r,t). (30)

where and . Under this transformation the hamiltonian becomes

 ^Q′=^U†(t)^Q^U(t)−i^U†(t)∂t^U(t) (31)

The second term of this expression cancels the Luttinger interaction (Eq. (27)), while the first term, applied to first order in , generates the coupling (Eq. (28)) between the thermal current and the thermal vector potential. Up to this point, we have simply made a transformation that may present some technical advantages Tatara (2015). However, we can also treat the thermal vector potential as a driving field in its own right, in which case it will have not only a longitudinal component (equivalent to the Luttinger’s thermomechanical potential), but also a transverse component. The interest of such transverse components is that they give us control on the transverse components of the thermal current: in other words, an interaction of the form (28) could be used as the basis for a thermal current density functional theory, in which the full thermal current, with longitudinal and transverse components, becomes a basic variable. We will not pursue this possibility in this review. We will formulate our thermal density functional theory in terms of the scalar thermomechanical potential. The resulting theory determines, in principle exactly, the energy density and the longitudinal part of the thermal current (i.e. the divergence of the thermal current), from which the practically important thermal energy fluxes can be computed.

## Iii Structure of thermal DFT

### iii.1 Static thermal Density-Functional Theory

Thermal DFT aims at the description of thermoelectric transport phenomena and therefore is a time-dependent or nonequilibrium theory. There are, however, two scenarios in which a static thermal DFT is of interest: 1) For situation in which already the initial state of the system is exposed to nonuniform temperatures. 2) For the derivation of so-called adiabatic approximations. For both it is required to consider a generalized equilibrium theory. We will refer to this theory as quasi-equilibrium theory due to the fact that is time-independent, but allows for a non-vanishing static thermomechanical potential.

#### iii.1.1 Quasi-equilibrium grand potential

Let us consider a the quasi-equilibrium grand potential as functional of the density matrix, ,

 Ω[^D]=Tr{^D[^Qv,ψ+1βlog^D]} , (32)

with Hamiltonians of the form

 (33)

The chemical potential, , and the inverse temperature, , define the global zero and scale of the energy, respectively. Mermin Mermin (1965) showed that

 Ω[^D~v,ψ]<Ω[^D] , (34)

with the quasi-equilibrium grand-canonical-ensemble density matrix given by

 (35)

where we have introduced

 ^H~v,ψ (36a) ~v(r) =[1+ψ(r)][v(r)−μ] . (36b)

It is important to stress that in Mermin’s proof of Eq. (34) no assumption about any particular Hamiltonian is required.111For details we refer to the original proof given in the appendix of Mermin, 1965. Accordingly the proof is valid also for the quasi-equilibrium Hamiltonians defined in Eq. (33). Using the corresponding quasi-equilibrium density matrix, Eq. (35), we can write the grand potential as functional of the external potential, and the thermomechanical potential , i.e.,

 ~Ω[~v,ψ] =−1βlogTr{e−β∫d3r[(1+ψ(r))^h(r)+~v(r)^n(r)]}. (37)

Moreover, we know the density and energy density as functionals of the potentials,

 n(r)=Tr{D~v,ψ^n(r)} , (38a) h(r)=Tr{D~v,ψ^h(r)} . (38b)

#### iii.1.2 Uniqueness of the static potentials–density mapping

Let us now define the potentials

 ~vλ(r) =(1−λ)~v0(r)+λ~v1(r) , (39a) ψλ(r) =(1−λ)ψ0(r)+λψ1(r) , (39b)

and denote the corresponding quasi-equilibrium density matrices by . Due to the variational principle (34) we have

 ~Ω[(1−λ)~v0+λ~v1,(1−λ)ψ0+λψ1]=Tr{^Dλ[^H~vλ,ψλ+1βlog^Dλ]} =(1−λ)Tr{^Dλ[^H~v0,ψ0+1βlog^Dλ]}+λTr{^Dλ[^H~v1,ψ1+1βlog^Dλ]} >(1−λ)Tr{^D0[^H~v0,ψ0+1βlog^D0]}+λTr{^D1[^H~v1,ψ1+1βlog^D1]} =(1−λ)~Ω[~v0,ψ0]+λ~Ω[~v1,ψ1] , (40)

which proves that is a strictly concave functional of the potentials. Hence we can employ a Legendre transformation to obtain a convex functional , i.e., a functional of the density, , and energy density , which are the conjugate variables to the external potentials and ,

 n(r) =δ~Ω[~v,ψ]δ~v(r) , (41a) h(r) =δ~Ω[~v,ψ]δψ(r) . (41b)

More precisely: , the quasi-equilibrium free energy, is defined as the negative of the Legendre transform of the grand potential Eq. (37),

 F[n,h] =~Ω[v[n,h],ψ[n,h]]−∫d3r(ψ[n,h](r)h(r)+v[n,h](r)n(r)) . (42)

Equations. (41a) and (41b) define the densities and as functionals of the external potentials and . The concavity of guarantees that we can invert these functional mappings to yield the external potentials as functionals of the densities, which, in fact, is the analogue of Mermin’s finite temperature DFT (FT-DFT) mapping theorem Mermin (1965) for static thermal DFT. From Eq. (42) follows immediately that

 ~v(r) =−δF[n,h]δn(r) , (43a) ψ(r) =−δF[n,h]δh(r) . (43b)

The quasi-equilibrium free energy functional of static thermal DFT is connected to the free energy functional of FT-DFT by,

 Feqβ[n] =F[n,heqβ[n]] . (44)

In equation (44) we have implicitly defined the equilibrium energy density (at inverse temperature ) as functional of the density, , using Mermin’s FT-DFT mapping theorem.

#### iii.1.3 Constrained Search Formulation

We conclude this section by presenting an alternative way to define the free-energy functional of static thermal DFT. We start by rewriting the grand potential Eq. (37) as

 (45)

where the minimization is taken over all statistical operators normalized to . Now we employ the so-called constrained-search formalism, which is an alternative route to define functionals in static DFT.Levy (1982); Lieb (1983) The idea of the constrained-search procedure is to split the minimization process into two steps: 1) Minimize over all statistical operators yielding prescribed densities and . 2) Minimize over all and . This leads to

 ~Ω[~v,ψ]=minn(r),h(r)∫d3r(ψ(r)h(r)+~v(r)n(r))+F[n,h] , (46)

where

 F[n,h] =∫d3rh(r)−1βS[n,h], (47a) S[n,h] (47b)

Eq. (46) emphasizes that and are related by a Legendre transformation and the set of potentials and the set of densities form a pair of conjugated variables. Moreover, Eq. (47a) highlights that the universal free-energy functional consists of a (trivial) part linear in the energy density and the entropy (times the reference temperature) as a functional of the densities. The entropy functional is defined as the maximum entropy compatible with the prescribed densities and . This definition of the entropy should be contrasted with equilibrium FT-DFT, where the entropy, as functional of the density alone, cannot be defined in terms of a constrained search, because the energy contribution to the free energy is no trivially given in terms of the density. Instead, only the free energy can be defined via a constrained search procedure,

 Feqβ[n]=min^D→n(r)⟨∫d3r^h(r)+1βlog^D⟩ . (48)

### iii.2 Time-dependent thermal Density-Functional Theory

Transport phenomena are intrinsically out of equilibrium and therefore a time-dependent description is required Di Ventra (2008). In the previous section we have explained the formal framework of static thermal DFT. Here, we describe the full-fledged time-dependent thermal DFT, suitable to address transport phenomena including the effects of retardation or history dependence.

#### iii.2.1 The action functional in thermal DFT

In order to setup the basic functionals we start by promoting the quasi-equilibrium grand potential, which can be viewed as generating functional for the quasi-equilibrium density and energy density, to an action functional,

 (49)

Formally definitions Eq. (49) and (37) are quite similar. The biggest difference is that the potentials and are time-dependent potentials. Therefore we have an additional time integral in the exponent which runs along the so-called Keldysh contour , depicted in Fig. 2, in the complex time plane Keldysh (1964, 1965); Stefanucci and van Leeuwen (2013). In Eq. (49) we have parametrized the contour by a real parameter , i.e., , and the contour ordering is defined w.r.t. the parameter . Note that we formally write the potentials as functions of the parameter , which allows them to be different for the same physical time depending on whether is on the forward or backward running part of the Keldysh contour.

The class of physical potentials is given by and that are constant in time along the vertical branch of the contour and have an identical time-dependence on the part of the contour running forward and backward on the real time axis. In the following we will adapt the convention that we refer to physical potentials whenever the time argument of the potentials is . Moreover, if we evaluate for physical potentials the contribution from the horizontal branch actually cancel and we get

 ~Λ[~v,ψ]=−iℏβ~Ω[~v|,ψ|] , (50)

which means that the value of the action functional is proportional to the grand potential functional introduced in Sec. (III.1). The subscript “” in Eq. (50) denotes that the -independent potential along the verical axis of the Keldysh contour is plugged into the grand potential functional . The value of the action functional does not provide any further information.

However, the action functional is the generating functional for the time-dependent density and energy density,

 n(r,t) =δ~Λ[~v,ψ]δ~v(r,t) , (51a) h(r,t) =δ~Λ[~v,ψ]δψ(r,t) . (51b)

Accordingly serves as tool to compute physically interesting expectation values. Moreover, it has been shown van Leeuwen (1998) that writing a generating functional employing the Keldysh contour resolves the apparent causality paradox van Leeuwen (1998); Vignale (2008), which plagued initial action functionals for TDDFT.

Now, under the assumption that Eqs. (51) can be inverted, we define another action functional via Legendre transformation, i.e.,

 A[n,h]≡~Λ[~v[n,h],ψ[n,ψ]]−∫C∫d3r(ψ[n,h](r,τ)h(r,τ)+~v[n,h](r,τ)n(r,τ)) . (52)

Note that in Eq. (52) all functional dependencies imply a history dependence. This means that the potentials and depend on the density and energy density at earlier times. From Eq. (52) it is easy to see that

 ~v(r,t) =−δA[n,h]δn(r,t) , (53a) ψ(r,t) =−δA[n,h]δh(r,t) . (53b)

Its is crucial to keep in mind that incorporates the dependence on the evolution of both densities, and , on the Keldysh contour. Put differently: the variation of yields time-dependent potentials, and , which depend on the history of both densities. However, the similarity to the free energy functional of static thermal DFT (cf. Sec. III.1) suggests a first approximation: Neglecting the history dependence all together. This so-called adiabatic approximation reads

The definition of the adiabatic functional implies that the exact action functional can be written as . In addition to the adiabatic action, , the dynamical excess entropy, , contributes to the full action. Note that the excess entropy is formally defined as the difference of the full action functional and the adiabatic approximation . We have opted to refer to the remainder as excess entropy, since the adiabatic approximation contains a trivial dependence on the time-dependent energy density. Accordingly we decompose the action functional into

Only the last part, the dynamical excess entropy , leads to memory, i.e., a dependence of the potentials on the densities at previous times. It is this part that will introduce retardation and dissipation effect in thermal DFT. In App. B we address the issue of the invertibility of Eqs. (51) by explicitly deriving the condition under which the linear response functions can be inverted.

#### iii.2.2 The Kohn-Sham system

One key aspect of DFTs is the so-called Kohn-Sham (KS) system Kohn and Sham (1965), a system of fictitious non-interacting electrons, which reproduces the density of the interacting system. The construction of the KS system in thermal DFT exhibits subtle differences compared to the usual KS scheme. Virtually every incarnation of DFT requires the densities in the interacting and non-interacting system to be equal. While this is still true for the electronic density, the energy densities of the interacting and non-interacting system are not identical in thermal DFT.

Let us start by noting that the operators yielding the energy density differ for interacting and non-interacting systems,

 ^h(r) =^t(r)+^u(r) , (56a) ^hs(r) =^t(r) . (56b)

On the one hand, the operator for the interacting system, , contains kinetic and interaction contributions, and , respectively. The operator for the non-interacting system, on the other hand, is only the kinetic energy density. Intuitively it seems rather cumbersome to reproduce the interacting energy density, containing interaction and kinetic contributions, with the non-interacting energy density, which is purely kinetic. Instead, considering the standard formulation of DFT, it is natural to write the interacting energy density as the kinetic energy density of the KS system plus the interaction energy density,

 h(r,t)=hs(r,t)+ϵHxc(r,t) . (57)

The energy density contains the Hartree contribution and a contribution due to exchange-correlation (xc) effects,

 ϵHxc(r) =ϵH(r)+ϵxc(r) , (58a) ϵH(r) =12n(r,t)vH(r,t) , (58b) vH(r,t) =∫d3r′n(r,t)w(|r−r′) . (58c)

The formal definition of the interaction energy density employs the energy density as functional of the density for a given reference temperature, which we already introduced at the end of Sec. III.1.2, i.e.,

 ϵHxc,β[n(t)](r)=heqβ[n(t)](r)−heqs,β[n(t)](r) . (59)

We can rearrange Eq. (57) using definition (59) in order to introduce the excess energy density,

 ¯h(r,t)=h(r,t)−heqβ[n(t)](r)=hs(r,t)−heqs,β[n(t)](r) . (60)

Equation (60) implies that the connection between the interacting and the non-interacting energy density, given in Eq. (57), requires the excess energy densities to be identical in both systems. Physically this means that the KS system is chosen to be out–of–equilibrium by the same amount as the interacting system. The metric for being out–of–equilibrium is the excess energy density, i.e., the difference of the time-dependent energy density and the instantaneous equilibrium energy density. The instantaneous equilibrium energy density, in turn, is determined by the density common to the interacting and the KS system.

The action functional for the non-interacting KS system, can be decomposed analogously to the interacting functional, cf. Eq. (55). The difference between the action functional of the interacting and the KS system is the so-called Hartree-exchange-correlation action,

Note that we have not specified whether we choose or as our energy-density variable. However, in DFT one usually obtains the expectation values directly from the KS system, which strongly suggests to use the KS energy density, . This means that the xc entropy, i.e., the sum of the adiabatic xc entropy and the dynamic xc entropy, is defined by

 ¯Sxc[n,hs]=¯S[n,hs+ϵHxc[n]]−¯Ss[n,hs] . (62)

Now we can write out explicitly the functional derivative w.r.t.  for the xc action. From Eq. (53b) it follows immediately that

 −1βδ¯Sxc[n,hs]δhs(r,t)=ψs(r,t)−ψ(r,t)=¯ψxc(r,t) , (63)

where the functional derivative is taken at constant density .

Some care has to be taken when differentiating w.r.t. the electronic density at constant . In definition (62) the energy-density argument of is shifted by an density-dependent amount. From Eq. (53a) we obtain

 ~vs(r,t)−~v(r,t)=δAHxc[n,hs]δn(r,t)−∫d3r′δA[n,hs+ϵHxc[n]]δhs(r′,t)δϵHxc[n(t)](r′)δn(r,t) . (64)

The second term on the right hand side is a counter term required to shift the energy argument of the interacting functional inside the density derivative. We point out that this term only appears if we construct explicit approximations for . As we will see later, in the adiabatic local density approximation the potentials are directly approximated and therefore the counter term will not appear. In the present discussion we formally derive the xc potentials from and, hence, we will keep the counter term. Taking into account decomposition (61) we arrive at

 ~vs(r,t)−~v(r,t) =δEHxc[n]δn(r,t)−1βδ¯Sadiaxc[n,hs]δn(r,t)−1βδ¯Sdynxc[n,hs]δn(r,t) =+∫d3r′ψ(r′,t)δϵHxc[n(t)](r′)δn(r,t) . (65)

Furthermore, using Eq. (61) to decompose the xc thermomechanical potential, , we finally arrive at the expression for the KS potentials:

 ~vs(r,t) =~v(r,t)+veqHxc(r,t)+¯vadiaxc(r,t)+¯vdynxc(r,t)+∫d3r′ψ(r′,t)δϵHxc[n(t)](r′)δn(r,t) , (66a) ψs(r,t) =ψ(r,t)+¯ψadiaxc(r,t)+¯ψdynxc(r,t) . (66b)

Both the external potential, , and the thermomechanical field, , have contributions from the adiabatic excess entropy, , and the dynamical excess entropy, . The external potential has an additional contribution corresponding to the instantaneous equilibrium potential,

 veqHxc(r,t) =δEeqHxc[n]δn(r,t) . (67)

The additional contribution due to the density-dependent shift of the energy argument is the potential energy associated with the thermomechanical potential coupled to the change of the instantaneous equilibrium energy density. Also this contribution is an adiabatic contribution which is only determined by the instantaneous density. Therefore we combine this contribution with Eq. (67) to define

 ~veqHxc(r,t)=vH(r,t)+veqxc(r,t)+∫d3r′ψ(r′,t)δϵHxc[n(t)](r′)δn(r,t) , (68)

where we also explicitly extracted the Hartree potential for clarity. This allows us to write the time-dependent KS equation in the compact form:

 iℏ∂tϕα(r,t)=[−ℏ2∇r1+ψ(r,t)+¯ψxc(r,t)2m∇r+~v(r,t)+~veqHxc(r,t)+¯vxc(r,t)]ϕα(r,t) . (69)

In Eq. (69) we further combined the adiabatic and dynamical excess contributions of the potential and thermomechanical field for brevity,

 ¯vxc(r,t) =¯vadiaxc(r,t)+¯vdynxc(r,t) , (70a) ¯ψxc(r,t) =¯ψadiaxc(r,t)+¯ψdynxc(r,t) . (70b)

The self-consistent solution of the time-dependent Schrödinger equation (69) together with the initial occupations of the orbitals yield the time-evolution of the density and the KS energy density, . Note that the presence of a thermomechanical potential bears close resemblance with a position dependent mass, which appear, e.g., in the description of compositionally grated semi-conductors Geller and Kohn (1993).

We conclude the section by recapitulating the required approximations for the implementation of the KS scheme in thermal DFT: First of all an approximation of is needed in order to compute the interacting energy density from the KS system. Furthermore implies an approximation for . Note that is, in principle, a functional already defined in Mermin’s FT-DFT. The really new features in the KS equation of thermal DFT are the contributions and , which can be obtained from an approximation to the excess entropy. They contain all the retardation effects induced by the interactions.

## Iv Alternative formulations

In the previous sections we have treated the phenomenon of heat and energy currents by considering the local temperature as an “internal” parameter of the theory. In other words, although the presence of a temperature somehow requires tracing out degrees of freedom of a bath, or a set of baths, we do not consider such a problem as an open one 222Note that, due to the presence of a temperature, we cannot classify it as a closed quantum problem either.. However, alternative approaches have been proposed in the literature that aim at treating explicitly thermoelectricity as an open quantum system problem. Here, we briefly review two of them. One approach solves the stochastic Schrödinger equation for the state vector (or, an equation of motion for the density matrix) and reformulates a TDCDFT in such a context. The other employs the stochastic Schödinger equation to describe the coupling of black bodies, i.e., idealized sources of incoherent light, to a quantum system. When the system is coupled to multiple black bodies at different temperatures they induce a thermal gradient and therefore a heat flow in the system.

### iv.1 Stochastic density functional theory for open systems

There are two ways in which one could formulate a TDCDFT for open quantum systems: via an equation of motion for the density matrix Burke et al. (2005) or an equation of motion for the state vector D’Agosta and Di Ventra (2007). Since the Hamiltonian of DFT is dependent on the density and/or current density, it depends on the state of the system, and hence it is generally stochastic D’Agosta and Di Ventra (2013). In this case then, a closed form of the equation of motion for the density matrix is not available, and one is forced to start with the equation of motion for the state vector D’Agosta and Di Ventra (2007). In the absence of memory, such an equation takes the form ()

 d|Ψ⟩=[−i^H|Ψ⟩−12^V†^V|Ψ⟩]dt+^V|Ψ⟩dW, (71)

where is a infinitesimal Wiener process (Itô calculus is employed), and is an operator that describes the interaction of the system with an environment (generalization to many environments simply adds a sum over different Wiener processes in Eq. (71)). The Hamiltonian is a general many-body Hamiltonian in the presence of an external vector potential .

The theorem of stochastic TDCDFT D’Agosta and Di Ventra (2007) then states that there is a one-to-one correspondence between the ensemble-averaged current density and the external vector potential, thus leading to the set of Kohn-Sham equations

 d|ΨKS⟩=(−i^HKS−12^V†^V)|ΨKS⟩dt+^V|ΨKS⟩dW (72)

where is a Slater determinant of single-particle wave-functions and

 ^HKS=N∑i=1[^pi+eA