Theory Department
Max Planck Institute of Microstructure Physics

Theory Department >> Research >> Development of functional theories for static and time-dependent quantum systems

Development of functional theories for static and time-dependent quantum systems

Research topics

Potential Funtional theory A. Cangi, C. Proetto

Functionals for non-collinear magnetism S. Sharma, K. Dewhurst, F. Eich, F. Essenberger

Basic rules for functional construction of finite temperature C. Bersier, C. Proetto, A. Sanna, A. Cangi, S. H. Mirhosseini

Self-consistent RPA for bond dissociation M. Hellgren, D. R. Rohr

Discontinuities of the exchange-correlation kernel M. Hellgren

Electronic structure via potential functional approximations

In the original form of density functional theory (DFT), suggested by Thomas and Fermi (TF) and made formally exact by Hohenberg and Kohn, the energy of a many-body quantum system is minimized directly as a functional of the density. Its modern incarnation uses the Kohn-Sham (KS) scheme, which employs the orbitals of a fictitious non-interacting system, such that only a small fraction of the total energy need be approximated. With the goal to simulate ever larger molecular and solid state systems, interest is rapidly reviving in finding an orbital-free approach to DFT. The major bottleneck in modern calculations is the solution of the KS equations, which can be avoided with a pure density functional for the kinetic energy of non-interacting fermions, T_s. Despite decades of effort, no generally applicable approximation for T_s has been found.

We consider the potential as a more natural variable to use in deriving approximations to quantum systems. In exact potential functional theory (PFT) we obtain the ground-state (gs) energy from

(1)

where

is the kinetic energy operator,

the inter-electronic repulsion,

the one-body potential, and Ψ are N-particle wavefunctions. Defining

(2)

with Ψ_ν denoting the gs wavefunction of potential ν(r), the exact gs energy is given by

(3)

n denoting the gs density. In practice, this requires two separate approximations, F^A[ν] and n^A[ν] yielding a direct approximation

In a recent letter ^[1] we go beyond those results by considering explicit potential functional approximations to interacting and non-interacting systems of electrons. We show that (i) the universal functional, F[ν], is determined entirely from knowledge of the density as a functional of the potential, such that only one approximation is required, namely n^A[ν], (ii) the variational principle imposes a condition relating energy and density approximations, (iii) a simple condition guarantees satisfaction of the variational principle, (iv) with an orbital-free approximation to the non-interacting density as a functional of the potential, T_s is automatically determined, i.e., there is no need for a separate approximation, and (v) satisfaction of the variational principle improves accuracy of approximations.

We deduce an approximation to F from any n^A[ν](r) by introducing a coupling constant in the one-body potential, ν^λ(r) = (1-λ)ν₀(r) +λν(r), where ν₀(r) is some reference potential. Via the Hellmann-Feynman theorem (and for the choice ν₀=0) we obtain

(4)

where

This establishes that the universal functional is determined solely by the knowledge of the density as functional of the potential. Moreover, insertion of n^A(r) on the right defines an associated approximate F^cc[n^A[ν]], where cc denotes coupling constant.

Furthermore, consider the variational principle. In PFT, this yields^[2]

(5)

with a possibly different value from E_ν^A,dir for a given pair of approximations {F^A, n^A}. We also show that this minimization over trial potentials

can be avoided, such that F^cc[n^A[ν]] does yield E_ν^A,var, if, and only if,

(6)

where

This condition is satisfied by the exact response function, but not necessarily by an approximation n^A[ν]; furthermore, enforcing Eq. (6) on n^A[ν] improves upon the accuracy.

Of much more practical use is the application of these results to the non-interacting electrons in a KS potential, mimicking inter-electronic interactions via the Hartree (ν_H) and exchange-correlation (ν_XC) contribution:

(7)

For a given

which determines

this equation can be easily solved by standard iteration techniques, bypassing the need to solve the KS equations via a given

where the subscript S denotes single-particle quantities. However, once self-consistency is achieved, we need T_s to extract the total energy of the interacting electronic system. All our derivations apply equally to the non-interacting problem, so we deduce:

(8)

which is the analog of Eq. (4) for non-interacting electrons in the external potential ν(r) (which is called ν_s(r) when it is the KS potential of some interacting system). This defines a kinetic energy approximation determined solely by the density approximation, T_s^cc[n_s^A[v]]. This result eliminates the need for constructing approximations to the non-interacting kinetic energy, T_s.

To illustrate the power of these results, we consider a simple example of non-interacting, spinless fermions in a one-dimensional box with potential ν(x). In Fig. (1) we plot exact and approximate kinetic energy densities (KED) demonstrating that our coupling-constant approach from Eq.(8) evaluated on n_s^A[ν] from Ref. [3] yields both locally and globally more accurate results, and required no separate approximation for the KED.

: Figure 1: Top panel: Exact (ex) KED of the coupling-constant (cc) and the direct (dir) approach (black and blue) with the corresponding semiclassical (sc) approximations (red and green) for one particle in ν(x) = -5 sin²(πx), 0 < x < 1 . Bottom panel: Absolute errors.

Potential functional approximations to the density, n_s^A[ν], are presently being developed via a systematic asymptotic expansion in terms of the potential, which has already been found in simple cases ^[3],
[4]. The leading terms are local approximations of the TF type, and the leading corrections, which are relatively simple, nonlocal functionals of the potential, greatly improve over the accuracy of local approximations in a systematic and understandable way. However, a practical realization of the presented approach awaits general-purpose approximations to n_s^A[v](r) for an arbitrary three-dimensional case. We have shown here that it no longer awaits the corresponding kinetic energy approximations.

References

[1]: A. Cangi, D. Lee, P. Elliott, K. Burke and E. K. U. Gross, Phys. Rev. Lett. 106, 236404 (2011).
[2]: E.K.U. Gross and C.R. Proetto, J. Chem. Theory Comput. 5, 844 (2009).
[3]: P. Elliott, D. Lee, A. Cangi, and K. Burke, Phys. Rev. Lett. 100, 256406 (2008).
[4]: A. Cangi, D. Lee, P. Elliott, and K. Burke, Phys. Rev. B 81, 235128 (2010).

To top

Functionals for non-colinear magnetism

Most formulations of spin density functional theory (SDFT) restrict the magnetization vector field to have global collinearity. Nevertheless, there exists a wealth of strong non-collinearity in nature, for example molecular magnets, spin-spirals, spin-glasses and all magnets at finite temperatures.

The local spin density approximation (LSDA) can be extended to these non-collinear cases^[1] but this extension has the undesirable property of having the exchange-correlation (xc) field parallel to the magnetization density at each point in space. When used in conjunction with the equation of motion for the spin magnetization in the absence of spin currents and external fields^[2],
[3] this local collinearity eliminates the torsional term, resulting in no time evolution. This severe shortcoming of LSDA, where the physical prediction is qualitatively wrong, opens up an important new direction for the development of functionals where this time evolution is correctly described.

Towards this goal we have extended the Kohn-Sham optimized effective potential (OEP) method to the non-collinear case and derive the corresponding integral equations, applicable to both finite and extended systems ^[3],

[4] To derive the OEP equations in the general non-collinear case, we start with the Kohn-Sham (KS) equation for two-component spinors Φ_i which has the form of a Pauli equation. For non-interacting electrons moving in an effective scalar potential v_s and a magnetic vector field B_s it reads as (atomic units are used throughout)

(1)

This equation can be derived by minimizing the total energy which, in SDFT, is given as a functional of the density

and the magnetization density

For a given external scalar potential v_ext and magnetic field B_ext this total energy reads

(2)

where

is the Hartree energy. The xc potential and xc magnetic field are given by

(3)

respectively. The exact functional form of E_xc[ρ,m] is unknown and has to be approximated in practice.

Assuming that the densities [ρ,m] are non-interacting (v,B)-representable one may, equivalently, minimize the total-energy functional (2) over the effective scalar potential and magnetic field. Thus the conditions

(4)

must be satisfied.

If the functional derivatives in Eq. 4 are evaluated for an xc functional that depends explicitly on the KS spinors, one obtains the natural extension of the OEP equations to non-collinear magnetism. By the usage of spinor valued wavefunctions we can stay within a single global reference frame, in contrast to the case where functionals originally designed for collinear magnetism are used in a non-collinear context by introducing a local reference frame at each point in space. The most commonly used orbital functional is the EXX energy given by

(5)

where the label occ indicates that the summation runs only over occupied states. In the following we restrict ourselves to an exchange-only treatment although generalization to other orbital functionals is straightforward.

For the energy functional Eq.(2) using the EXX approximation to E_xc one obtains the following coupled integral equations for the exchange potential and magnetic field

(6)

and

(7)

where

and j runs only over the unoccupied states. The matrix Λ is given by

(8)

where

(9)

are the non-local matrix elements of the Coulomb interaction between states i and j.

In order to explore the impact of treating non-collinear magnetism in the way outlined above and at the same time to ensure that our numerical analysis be as accurate as possible, we implemented the OEP equations for the fully non-collinear case within the full-potential linearized augmented plane wave (FP-LAPW) method implemented within the Elk code ^[5]. This method is then applied to study the non-collinear spin magnetism in an unsupported Cr (111) monolayer. In Fig. 1 are shown the magnetization density and B field calculated using both the LSDA and the OEP method.

: Fig. 1: Fully non-collinear magnetization density and B field obtained using the LSDA and exchange-only EXX functionals for an unsupported

A striking feature of the OEP B field is that, unlike its LSDA counterpart, it is not locally parallel to the magnetization density and this will produce manifestly different spin-dynamics. This is because the equation of motion for the spin magnetization reads

(10)

where J_s is the spin current and γ the gyromagetic ratio. In the time-independent LSDA and conventional GGA, m(r) and B_xc(r) are locally collinear, as is clear from Fig. 1, and therefore m(r) x B_xc(r) vanishes. This also holds true in the adiabatic approximation of time dependent SDFT which, by Eq. (10), implies that these functionals cannot properly describe the dynamics of the spin magnetization. In contrast, already at the static level, for the EXX functional m(r) x B_x(r) does not vanish. In fact, in the ground state of a non-collinear ferromagnet without external magnetic field, m(r) x B_xc(r) exactly cancels the divergence of the spin current, ∇. J_s, i.e. these terms are equally important, and it is essential to have a proper description of m(r) x B_xc(r). These results indicate that a time-dependent generalization of our method could open the way to an ab-initio description of spin dynamics. How well this functional really performs in describing the spin dynamics remains a question for future investigations.

To top

References

[1]: J. Kuebler, K.-H. Hoeck, J. Sticht and A. R. Williams, J. Phys. F18, 469 (1993).
[2]: K. Capelle, G. Vignale and B. L. Gyoerffy, Phys. Rev. Lett.87, 206403 (2001).
[3]: S. Sharma, J. K. Dewhurst, C. Ambrosch-Draxl, S. Kurth, N. Helbig, S. Pittalis, S. Shallcross, L. Nordstroem and E.K.U. Gross Phys. Rev. Lett.98, 196405 (2007)
[4]: S. Sharma, S. Pittalis, S. Kurth, S. Shallcross, J. K. Dewhurst and E.K.U. Gross Phys. Rev. B76, 100401 (Rapid Comm.) (2007)
[5]: http://elk.sourceforge.net

To top

Basic rules for functional construction at finite temperature

Because of the small mass ratio between electrons and nuclei, standard electronic structure calculations treat the former as being in their ground state, but routinely account for the finite temperature of the latter, as in ab initio molecular dynamics^[1]. But as electronic structure methods are applied in ever more esoteric areas, the need to account for the finite temperature of electrons increases. Phenomena where such effects play a role include rapid heating of solids via strong laser fields^[2], dynamo effects in giant planets ^[3], magnetic^[4] and superconducting phase transitions ^{[5],
[6]} , shock waves ^[7] , warm dense matter ^[8] , and hot plasmas ^[9].

Within density functional theory, the natural framework for treating such effects was created by Mermin ^[10]. Application of that work to the Kohn-Sham (KS) scheme at finite temperature also yields a natural approximation: treat KS electrons at finite temperature, but use ground-state exchange-correlation (XC) functionals. This works well in recent calculations ^{[7]
,
[8]}, where inclusion of such effects is crucial for accurate prediction. This assumes that finite-temperature effects on exchange-correlation are negligible relative to the KS contributions, which may not always be true.

In the present work we establish the basic rules that will allow the building of finite-temperature exchange-correlation functionals neyond these standard approximations ^[11].

Central to the thermodynamic description of many-electron systems is the grand-canonical potential, defined as the statistical average of the grand-canonical operator

(1)

where Ĥ, Ŝ, ^N, τ and μ are the Hamiltonian, entropy, and particle-number operators, temperature and chemical potential, respectively. In detail,

where ^T and ^V_ee are the kinetic energy and the Coulomb electron-electron interaction operators, and ^V represents an external scalar potential v(r). The entropy operator is given by

where k is the Boltzmann constant and

is a statistical operator, with | Ψ_N,i⟩ and w_N,i being orthonormal N-particle states and statistical weights, respectively, with the latter satisfying the (normalization) condition Σ_N,iw_N,i = 1. The statistical average of an operator ^A is obtained as

(2)

To create a DFT at finite temperature, Mermin ^[10] rewrites this as (in modern parlance)

(3)

where the minimizing n(r) is the equilibrium density n⁰(r), and

(4)

is the finite-temperature analog of the universal Hohenberg-Kohn functional, defined through a constrained search^[15]. This depends only on τ and not on μ. We denote

as the minimizing statistical operator in Eq. (4), and define the density functionals:

i.e., each density functional is the trace of its operator over the minimizing ^Γ for the given τ and density. Because it arises so often in this work, we have defined the "kentropy" K^τ[n]. Next consider a system of non-interacting electrons at the same temperature τ and denote its one-body potential as v_s(r).

All the previous considerations apply, and we choose v_S(r) to make its density match that of the interacting problem. This defines the KS system at finite temperature. The corresponding functionals will be labeled with the subscritp s.

We have found that the exchange (x) and correlation (c) parts of the grand-canonical potential, as well as the kentropy and the potential energy (U) are bound to fulfil the exact inequalities ^[11] :

(5)

and no approximation should violate these basic rules.

Some of the most important results in ground-state DFT come from uniform scaling of the coordinates ^[17]. In the following considerations, when we refer explicitly to wavefunctions, we shall restrict to wavefunctions having finite norm on their entire domain of definition. Under norm-preserving homogeneous scaling of the coordinate r→γr, with γ > 0, to the scaled wave function ^[17].

(6)

corresponds the scaled density n_γ(r) = γ³n(γr). A simple 1D sketch of the scaling procedure is given in this figure:

Writing Ψ^γ(r₁,....r_N) = ⟨r,...,r_N|Ψ^γ⟩ in terms of the (representation-free) element |Ψ^γ⟩ of Hilbert space, the scaled statistical operator is defined as^[11]

(7)

where the statistical weights are hold fixed, i.e., the scaling only acts on the states.

With the above definition, the statistical average of an operator whose pure-state expectation value scales homogeneously , scales homogeneously ^[17] as well. In particular, we have:

The scaling behavior of the density functionals is, however, more subtle. First consider the non-interacting functionals in some detail. We observe that

(8)

which implies

(9)

For non-interacting electrons, the statistical operator at a given temperature that is the minimizer for a given scaled density is simply the scaled statistical operator, but at a quadratically scaled temperature, an effect that is obviously absent in the ground-state theory.

Next, we consider the interacting case. The analysis of the scaling relations allows to derive a set of inequalities that provide tight constraints on the functionals, and are used in non-empirical functional construction in the ground state. We obtain the following:

(10)

(11)

(12)

which are reversed if γ < 1.

Lastly, we consider the adiabatic coupling constant for finite temperature, its relationship to scaling, and derive the adiabatic connection formula. Define

(13)

with

being the corresponding minimizing

Of course, non-interacting functionals are not affected by a coupling constant modification. Eq. (10) implies that the exchange and Hartree density functionals have a linear dependence on λ. Employing the minimization property of Eq. (13) and the Hellmann-Feynman theorem, we find^[11]

(14)

where

(15)

Eq.(14) is the finite-temperature adiabatic connection formula, whose zero-temperature limit played a central role in ground-state DFT.

(Eq. 5), and the scaling inequalities can be combined, analogously to Ref. [17] to show that U_xc^τ[n](λ) is monotonically decreasing in λ , as shown in the schematic drawing below.

: Fig. 1. Geometrical interpretation of the adiabatic connection formula. The shaded area between the curve and the λ axis represents Ω_xc^τ[n]. All magnitudes are negative, in agreement with the inequalities in Eq. (5).

So far, all results presented have been exact. To see them in practice, consider the finite-temperature local density approximation (LDA) to Ω_xc^τ[n]

(16)

where ω_xc^unifτ(n) is the XC grand canonical potential density of a uniform electron gas of density n. Because a uniform electron gas is a quantum mechanical system, its energies satisfy all our conditions, guaranteeing by construction that LDA satisfies all the exact conditions listed here. In the Jacob's ladder of functional construction^[13], more sophisticated approximations should also satisfy these conditions. To give one simple example of the usefullness of our results, Eq. (10) implies

(17)

where

is a dimensionless measure of the local temperature. Thus the largest fractional deviations from ground-state results should occur (in LDA) in regions of lowest density, but these contribute less in absolute terms.

In summary, there is a present lack of approximate density functionals for finite temperature. We have derived many basic relations needed to construct such approximations, and expect future approximations to either build these in, or be tested against them. In principle, such approximations should already be implemented in high-temperature DFT calculations, at least at the LDA level, as a check that XC corrections due to finite temperature do not alter calculated results. If they do, then more accurate approximations than LDA will be needed to account for them.

To top

References

[1]: R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985).
[2]: J. Gavnholt, A. Rubio, T. Olsen, K. S. Thygesen, and J. Schiøtz, Phys. Rev. B 79, 195405 (2009).
[3]: R. Redmer, T. R. Mattsson, N. Nettelmann, and M. French, Icarus 21, 798 (2011).
[4]: N. M. Rosengaard and B. Johansson, Phys. Rev. B 55, 14975 (1997).
[5]: G. Profeta, C. Franchini, N. N. Lathiotakis, A. Floris, A. Sanna, M. A. L. Marques, M. Lüders, S. Massidda, E. K. U. Gross, and A. Continenza, Phys. Rev. Lett. 96, 047003 (2006).
[6]: P. Cudazzo, G. Profeta, A. Sanna, A. Floris, A. Continenza, S. Massidda, and E. K. U. Gross, Phys. Rev. Lett. 100, 257001 (2008).
[7]: S. Root, R. J. Magyar, J. H. Carpenter, D. L. Hanson, and T. R. Mattsson, Phys. Rev. Lett. 105, 085501 (2010).
[8]: A. Kietzmann, R. Redmer, M. P. Desjarlais, and T. R. Mattsson, Phys. Rev. Lett. 101, 070401 (2008).
[9]: M. W. C. Dharma-wardana and M. S. Murillo, Phys. Rev. E 77, 026401 (2008).
[10]: N. D. Mermin, Phys. Rev. 137, A1441 (1965).
[11]: S. Pittalis, C. R. Proetto, A. Floris, A. Sanna, C. Bersier, K. Burke, and E. K. U. Gross, Phys. Rev. Lett. **** , ****** (2011).
[12]: R. G. Dandrea, N. W. Ashcroft, and A. E. Carlsson, Phys. Rev. B 34, 2097 (1986).
[13]: C. Fiolhais, F. Nogueira, and M. A. Marques, eds., A Primer in Density Functional Theory (Springer, 2003).
[14]: D. C. Langreth and J. P. Perdew, Solid State Commun. 31, 567 (1979).
[15]: R. Parr and X. Yang, Density-functional Theory of Atoms and Molecules (Oxford University Press, 1989).
[16]: M. Greiner, P. Carrier, and A. Görling, Phys. Rev. B 81, 155119 (2010).
[17]: M. Levy and J. P. Perdew, Phys. Rev. A 32, 2010 (1985).

To top

Correlation potentials for molecular bond dissociation within the self-consistent random phase approximation

The random phase approximation (RPA) for the total energy in density functional theory (DFT) incorporates exchange effects exactly and the correlation energy is treated non-perturbatively by summing a subset of polarization diagrams to infinite order. The RPA has the further advantage of allowing for a systematic improvement using the adiabatic connection fluctuation dissipation (ACFD) theorem in conjunction with the exchange-correlation (XC) kernel of time-dependent density functional theory.

For a unique solution the RPA should be solved self-consistently. This means that also the XC potential has to be obtained. In this work we have calculated self-consistent RPA potentials for diatomic molecules and investigated the RPA correlation potentials at different inter-atomic distances^[1]. By increasing the inter-atomic distance correlation effects become more dominant and the XC potential can give a clear image of the performance of the functional in this respect. In particular, the LiH molecule allows us to study the correlation effects behind the derivative discontinuity as well as the so-called fractional spin error. Within the ACFD framework the exact correlation energy is written as

(1)

where χ_s is the non-interacting KS density response function and χ_λ is the scaled interacting density response function. The scaling refers to a system with a linearly scaled Coulomb interaction λv(r,r') plus a fictitious potential which keeps the density fixed as λ is changed. The parameter λ runs between 0 (non-interacting KS case) and 1 (fully interacting case). We have used the short hand notation

for any two-point functions f and g. Within TDDFT the function χ^λ reads

(2)

The scaled XC kernel f_xc^λ is defined as the functional derivative of the scaled XC potential v_xc^λ with respect to the density n, evaluated at the ground state density.

In the RPA, f_xc^λ = 0 and thus corresponds to the simplest approximation within the ACFD formalism. Within RPA the λ-integral in Eq. (1) can be evaluated analytically with the result

(3)

Diagrammatically Eq. (3) is equal to an infinite summation of ring-diagrams.

The RPA correlation potential v_c can be obtained as the functional derivative of Eq. 3 with respect to the density. If we let V_s signify the total KS potential, G_s the non-interacting KS Green function and χ_s= -iG_sG_s, the functional derivative can be obtained via the chain rule

(4)

The result is the well-known linearized Sham-Schlüter (LSS) equation

(5)

Here, we have used the notation (r₁,t₁) = 1 etc. and introduced Λ(3,2;1)=-iG_s(3,1)G_s(1,2). The correlation part of the self-energy Σ_c in the RPA is given by

(6)

where

(7)

In Fig. 1 we show the RPA correlation potential of 1D LiH (Z_Li=3.6,Z_H=1.2). The solid, dotted and dashed curves represent the correlation potentials for 2, 3, and 8 Bohr interatomic distances, respectively. The build up of the peak at the bond midpoint is apparent. However, a step due to the derivative discontinuity, as is present in the exact KS correlation potential, is not observed.

: Fig. 1: The correlation potential along the bond with the bond midpoint at zero. The system is 1D LiH with a soft coulomb potential. The Li atom is located at -1, -1.5, and -4 Bohr, respectively. The H atom is located at 1, 1.5, and 4 Bohr, respectively.

Fig. 2 shows the correlation potential of the three dimensional LiH for bond distances 3 (solid), 6 (dotted), and 8 Bohr (dashed). The same features as in the 1D case are also found in the results for the 3D correlation potentials. In the region of the Li atom (-10 to 0) the potential qualitatively resembles that of the 1D system in Fig. 1. We see a well with a peak close to the nucleus. Also in the region of the H atom (0 to 10) figures 1 and 2 show the same structure. A well emerges with increasing bond distance. The difference between 1D and 3D is found only at the bond midpoint. A peak emerges only for the 1D system (Fig. 1). In contrast, the 3D system (Fig. 2) exhibits a peak at zero only for small and intermediate bond distances (3 and 5 Bohr). The orbital and potential basis functions are simply not diffuse enough to extend to the bond midpoint.

: Fig. 2: The correlation potential of LiH along the bond axis. The bond midpoint is at zero. The Li atom is located at -1.5, -2.5, and -4 Bohr, respectively. The H atom is located at 1.5, 2.5, and 4 Bohr, respectively.

We have further analyzed the RPA functional in the context of fractional charge and spin by calculating the total energy at different particle numbers. Fig. 3 shows the total energy as a function of the number of electrons at the H atom. The number of electrons at the Li atom will then be four minus the x value. The exact functional (black) has a minimum at 1.0, because it dissociates LiH into a neutral H atom and a neutral Li atom.

: Fig. 3: The total energy of dissociated LiH as a function of number of electrons at the H atom for EXX (blue), RPA (red), and exact (black).

In Fig. 3 we see that EXX (blue) and RPA (red) do not dissociate LiH into the neutral atoms. In both cases there is a surplus of electronic charge at the H atom. However, in the case of RPA the surplus (0.16 electrons) is much smaller than in the case of EXX (0.4 electrons). The large improvement may be related to the peak that is present in the RPA correlation potential. The total energy in the dissociation limit is shown to be largely improved in the RPA showing only a small fractional spin error.

References

[1]: M. Hellgren, D.R. Rohr, E.K.U. Gross, arXiv:1110.6062 (2011).

To top

Discontinuities of the exchange-correlation kernel and charge-transfer excitations in time-dependent density functional theory

Density functional theory (DFT) owes its great success to the fact that the highly intricate many-body wave function is replaced by the much simpler electronic density. In addition, the density is calculated from a system of non-interacting electrons known as the Kohn-Sham (KS) system. The effective KS potential is composed of the external, the Hartree and the so-called exchange-correlation (XC) potential v_xc, where the latter is given by the functional derivative of the XC energy E_xc with respect to the density, and incorporates all the effects of the Coulomb interaction beyond the Hartree level. The difficulties in using DFT lies in finding good approximations to E_xc. Many works have therefore been devoted to study its exact properties. One of the most intriguing features is the existence of derivative discontinuities at integral particle numbers^[1], a property which turns out to have several important physical implications.

It is not hard to see that a derivative discontinuity in E_xc gives rise to a singularity in v_xc in the form of a jump by a constant Δ_xc, i.e., v_xc⁺(r)=v_xc ^-(r)+Δ_xc, where ± refers to N→ N₀^±, and N₀ is an integer. A classic example where the jump in v_xc becomes important is in the dissociation of closed-shell molecules composed of open-shell atoms. In order to get the correct dissociation limit a step in v_xc has to develop over the atom with the lower ionization energy (a similar situation but with closed-shell fragments is depicted in Fig. 1).

: Fig. 1: Left: f_xc and v_xc of a dissociating He-Be²⁺ system. Right: Discontinuity of the Be²⁺ atom.

With the advent of time-dependent DFT (TDDFT) the variation of v_xc with respect to the density, i.e., the XC kernel f_xc(r,r',ω), has become another central quantity of interest. The reason is the possibility of extracting excitation energies from the linear density response function, which in TDDFT can be expressed in terms of the KS response function χ_s, the Coulomb interaction v, and f_xc

(1)

This equation represents an enormous simplification to the problem of calculating neutral excitation energies ^[2] Again, however, the difficulty lies in finding approximations to f_xc.

In this work we have evaluated the discontinuities of f_xc in order to increase the understanding of its exact behavior ^[3]. The discontinuities of f_xc turned out to be of a more complex nature than those of v_xc , where the general structure is an r -and ω-dependent function g_xc(r,ω), which may diverge. A novel idea to obtain an exact expression for the discontinuity is to study the limiting behavior of E_xc as a function of particle number N around an integer N₀, making explicit the local nature of the KS system described by v_xc compared to the actual nonlocal density dependence in E_xc. The constant jump in v_xc can then formally be written as

(2)

where f(r)=δn(r)/ δN is the Fukui function.The same idea can then be used for the static XC kernel and we find an expression for g_xc

(3)

After performing a common denominator approximation to the response functions it is easy to see that the second term exhibits a diverging behavior

which depends on the difference between the ionization energy I and the KS affinity A_s. This divergency can be seen in f_xc of the dissociating system in Fig. 1 where the calculations have been made in the exact-exchange approximation, which is known to have a discontinuity. The right panel of Fig. 1 displays g_xc, calculated in the same approximation but defined on an ensemble which allows for fractional charges. A striking similarity can be seen and we can conclude that the peaked structure in the XC kernel of a dissociating system is just the discontinuity.

A diverging behavior of the kernel is exactly what is needed in order to describe charge-transfer (CT) excitations as has been discussed in the literature for more than a decade. The reason is that Eq. (1) involves only matrix elements of f_xc between so-called excitation functions, i.e., products of occupied and unoccupied KS orbitals. As the distance between the fragments increases these products vanish exponentially and thus there is no correction to the KS eigenvalue differences unless the kernel diverges. The formulation of TDDFT for fractional charges has to allow the particle number to change in time and we have therefore proposed an ensemble of the form

where Ψ_k(t) is the k-particle many-body state at time t and α_k(t) are given TD probabilities. A Runge-Gross theorem can be proved also in this case allowing for an analysis of the same kind as in the static case. The function g_xc acquires anω-dependence and Fig. 2 illustrates the difference compared to the static case showing mainly an increase of the exponential growth. This increase is, however, crucial for the description of CT excitations, as seen in Fig. 3, where AEXX uses f_x(ω=0) and TDEXX uses f_x(ω_KS).

: Fig. 2: v_xc and f_xc at ω=0 and ω_KS of a dissociating He-Be system.

: Fig. 3: CT excitation energies of the system in Fig. 2 in different approximations.

In conclusion we have determined the discontinuities of the XC kernel and shown that they are a key feature in describing CT excitations. This somewhat unphysical feature of the KS quantities v_xc and f_xc reflects the fact that we are dealing with electrons far from a non-interacting description. Almost all existing approximate functionals miss the derivative discontinuity and therefore the construction of such approximations remains as one of the biggest challenges in DFT and TDDFT.

References

[1]: J. P. Perdew et al., PRL 49, 1691 (1982)
[2]: E. K. U. Gross et al., PRL 55, 2850 (1985)
[3]: M. Hellgren and E. K. U. Gross, arXiv:1108.3100 (2011)

To top