{'fraction_non_alphanumeric': 0.034596245750603535, 'fraction_numerical': 0.007380401044489333, 'mean_word_length': 4.422129614788695, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spectral Density Functionals for Electronic Structure Calculations

We introduce a spectral density functional theory which can be used to compute energetics and spectra of real strongly--correlated materials using methods, algorithms and computer programs of the electronic structure theory of solids. The approach considers the total free energy of a system as a functional of a local electronic Green function which is probed in the region of interest. Since we have a variety of notions of locality in our formulation, our method is manifestly basis--set dependent. However, it produces the exact total energy and local excitational spectrum provided that the exact functional is extremized. The self--energy of the theory appears as an auxiliary mass operator similar to the introduction of the ground--state Kohn--Sham potential in density functional theory. It is automatically short--ranged in the same region of Hilbert space which defines the local Green function. We exploit this property to find good approximations to the functional. For example, if electronic self--energy is known to be local in some portion of Hilbert space, a good approximation to the functional is provided by the corresponding local dynamical mean--field theory. A simplified implementation of the theory is described based on the linear muffin--tin orbital method widely used in electronic strucure calculations. We demonstrate the power of the approach on the long--standing problem of the anomalous volume expansion of metallic plutonium.

which is able to interpolate between well separated and well overlapping atomic orbitals. The solutions of many-body equations have to be carried out on the level of the Green functions which contain necessary information about the total energy and the spectrum of the solid.
The development of such techniques has a long history in condensed matter physics. Studies of strongly correlated systems have traditionally focused on model Hamiltonians using techniques such as diagrammatic methods 3 , Quantum-Monte Carlo simulations 4 , exact diagonalizations for finite-size clusters 5 , density matrix renormalization group methods 6 and so on. Model Hamiltonians are usually written for a given solid-state system based on physical grounds. In the electronic-structure community, the developments of LDA+U 7 and self-interaction corrected (SIC) 8 methods , many-body perturbative approaches based on GW and its extensions 2 , as well as time-dependent version of the density functional theory 9 have been carried out. Some of these techniques are already much more complicated and time-consuming comparing to the standard LDA based algorithms, and the real exploration of materials is frequently performed by its simplified versions by utilizing such, e.g., approximations as plasmon-pole form for the dielectric function 15 , omitting self-consistency within GW 2 or assuming locality of the GW self-energy 16 .
In general, diagrammatic methods are most accurate if there is a small parameter in the calculation, say, the ratio of the on-site Coulomb interaction U to the band width W . This does not permit the exploration of real strongly correlated situations, i.e. when U/W ∼ 1. Systems near Mott transition is one of such examples, where strongly renormalized quasiparticles and atomic-like excitations exist simultaneously. In these situations, self-consistent methods based on the dynamical mean-field based theory (DMFT) 10 , and its cluster generalizations such as dynamical cluster approximation (DCA) 17 , or cellular dynamical mean field theory (C-DMFT) 18,19 , are the minimal many body techniques which have to be employed for exploring real materials.
Thus, a combination of the DMFT based methods with the electronic structure techniques is promising, because a realistic material-specific description where the strength of correlation effects is not known a priori can be achieved. This work is in its beginning stages of development but seems to have a success. The development was started 20 by introducing so-called LDA+DMFT method and applying it to the photoemission spectrum of La 1−x Sr x TiO 3 . Near Mott transition, this system shows a number of features incompatible with the one-electron description 21 . The LDA++ method 22 has been discussed, and the electronic structure of Fe has been shown to be in better agreement with experiment than the one based on LDA. The photoemission spectrum near the Mott transition in V 2 O 3 has been studied 23 , as well as issues connected to the finite temperature magnetism of Fe and Ni were explored 24 . LDA +DMFT was recently generalized to allow computations of optical properties of strongly correlated materials 25 . Further combinations of the DMFT and GW methods have been proposed 12,26,27 and a simplified implementation to Ni has been carried out 27 .
Sometimes the LDA+DMFT method 11 omits full self-consistency. In this case the approach consists in deriving a model Hamiltonian with parameters such as the hopping integrals and the Coulomb interaction matrix elements extracted from an LDA calculation. Tight-binding fits to the LDA energy bands or angular momentum resolved LDA densities of states for the electrons which are believed to be correlated are performed. Constrained density functional theory 28 is used to find the screened on-site Coulomb U and exchange parameter J. This information is used in the downfolded model Hamiltonian with only active degrees of freedom to explore the consequences of correlations. Such technique is useful, since it allows us to study real materials already at the present stage of development. A more ambitious goal is to build a general method which treats all bands and all electrons on the same footing, determines both hoppings and interactions internally using a fully self-consistent procedure, and accesses both energetics and spectra of correlated materials.
Several ideas to provide a theoretical underpinning to these efforts have been proposed. The effective action approach to strongly correlated systems has been used to give realistic DMFT an exact functional formulation 29 . Approximations to the exact functional by performing truncations of the Baym-Kadanoff functional have been discussed 30 . Simultaneous treatment of the density and the local Green function in the functional formulation has been proposed 12 . Total energy calculations using LDA+DMFT have recently appeared in the literature 31,32,33,34 . DMFT corrections have been calculated and added to the LDA total energy in order to explain the isostructural volume collapse transition in Ce 31 . Fully self-consistent calculations of charge density, excitation spectrum and total energy of the δ phase of metallic Plutonium have been carried out to address the problem of its anomalous volume expansion 32 . The extensions of the method to compute phonon spectra of correlated systems with the applications to Mott insulators 33 and high-temperature phases of Pu 34 have been also recently developed.
In this paper we discuss the details of this unified approach which computes both total energies and spectra of materials with strong correlations and present our applications for Pu. We utilize the effective action free energy approach to strongly correlated systems 29,30 and write down the functional of the local Green function. Thus, a spectral density functional theory (SDFT) is obtained. It can be used to explore strongly correlated materials from ab inito grounds provided useful approximations exist to the spectral density functional. One of such approximations is described here, which we refer to as a local dynamical mean field approximation. It is based on extended 35 and cluster 17,18,19 versions of the dynamical mean-field theory introduced in connection with the model-Hamiltonian approach 10 .
Implementation of the theory can be carried out on the basis of the energy-dependent analog for the one-particle wave functions. These are useful for practical calculations in the same way as Kohn-Sham particles are used in density functional based calculations. The spectral density functional theory in its local dynamical mean field approximation, requires a self-consistent solution of the Dyson equations coupled to the solution of the Anderson impurity model 36 either on a single site 10 or on a cluster 17,18 . Since it is the most time-consuming part of all DMFT algorithms, we are carrying out a simplified implementation of it based on a slave boson Gtuzwiller 37,38,39 and Hubbard I 14,40 methods. This is described in detail in a separate publication 41 . We illustrate the applicability of the method addressing the problem of δ−Pu. Various aspects of the present work have appeared already 12,32 .
Our paper is organized as follows. In Section II we describe the spectral density functional theory and discuss local dynamical mean field approximation which summarizes the ideas of cluster and extended 35 versions of the DMFT. We show that such techniques as LDA+DMFT 11 , LDA+U 7 , and local GW 12,16 methods are naturally seen within the present method. Section III describes our implementation of the theory based on the energy-resolved one-particle description 20 and linear-muffin-tin orbital method 42,43,44 for electronic structure calculation. Section IV discusses application of the method to the volume expansion of Pu. Section V is the conclusion.

II. SPECTRAL DENSITY FUNCTIONAL THEORY
Here we discuss the basic postulates and approximations of spectral density functional theory. The central quantity of our formulation is a "local" Green function G loc (r, r ′ , z), i.e. a part of the exact electronic Green function which we are interested to compute. This is by itself arbitrary since we can probe the Green function in a portion of a certain space such, e.g., as reciprocal space or real space. These are the most transparent forms where the local Green function can be defined. We can also probe the Green function in a portion of the Hilbert space. If a function can be expanded in some basis set {χ ξ } our interest can, e.g, be associated with diagonal elements of the matrix F ξξ ′ (z).
As we see, the locality is a basis set dependent property. Nevertheless, it is a very useful property because a most economical description of the function can be achieved. This is true when the basis set which leads to such description of the function is known. The choice of the appropriate Hilbert space is therefore crucial if we would like to find an optimal description of the system with the accuracy proportional to the computational cost. In spectral density functional theory that has a meaning of finding good approximations to the functional. Therefore we always rely on a physical intuition when choosing a particular representation which should be tailored to a specific physical problem.
At the beginning we formulate spectral density functional theory in completely real space but keep in mind that such formulation is not unique. Thus, we are interested in finding a part of the electronic Green function restricted within a certain cluster area. Due to translational invariance of the Green function on the original lattice given by primitive translations {R}, i.e. G(r + R, r ′ + R, z) = G(r, r ′ , z), it is always sufficient to consider r lying within a primitive unit cell Ω c positioned at R = 0. Thus, r ′ travels within some area Ω loc centered at R = 0. We set the local Green function to be the exact Green function G(r, r ′ , z) within a given cluster Ω loc and zero outside. In other words, where the theta function is a unity when vector r ∈ Ω c , r ′ ∈ Ω loc and zero otherwise. It is schematically illustrated on Fig. 1. This construction can be translationally continued onto entire lattice by enforcing the property θ loc (r + R, r ′ + R) = θ loc (r, r ′ ).
We will now discuss the free energy of a system as a functional of the local Green function.

A. Functional of Local Green Function
We consider full many-body Hamiltonian describing the electrons moving in the periodic ionic potential V ext (x) = V ext (r)δ(τ ) and interacting among themselves according to the Coulomb law: v C (x − x ′ ) = e 2 /|r − r ′ |δ(τ − τ ′ ) [we use imaginary time-frequency formalism, where x = (r, τ )]. This is the formal starting point of our all-electron first-principles calculation. So, the theory of everything is summarized in the action S:.
FIG. 1: Illustration of the area in real space where the local Green function G loc is defined. Note that r is restricted by the unit cell at the origin while r ′ and r ′′ travel within the crystal.
(atomic Rydberg units, = 1, m e = 1/2, are used throughout). We will ignore relativistic effects in this action for simplicity but considering our applications to Pu, these effects will be included later in the implementation. In addition, the effects of electron-phonon interaction will not be considered. We will take the effective action functional approach to describe our correlated system 30 . The approach allows to obtain the free energy of the solid from a functional Γ evaluated at its stationary point. The main question is the choice of the variable of the functional which is to be extremized. This question is highly non-trivial because the exact form of the functional is unknown and the usefulness of the approach, depends on our ability to construct good approximations to it, which in turn depends on the choice of variables. The Baym-Kadanoff (BK) Green function theory considers exact Green function G(x, x ′ ) = − T τ ψ(x)ψ + (x ′ ) as a variable, i.e. Γ BK [G]. Density functional theory considers density ρ(r) = G(r, r, τ = 0) of the solid as a variable, i.e. Γ DF T [ρ]. Spectral density functional theory will consider local Green function G loc (x, x ′ ) = G(x, x ′ )θ loc (r, r ′ ) as a variable, i.e. Γ SDF [G loc ].
Notice on the variety of choices we can make, in particular in the functional Γ SDF [G loc ] since the definition of the locality is up to us. The usefulness of a given choice is dictated by the existence of good approximations to the functional, as, for example, the usefulness of the DFT is the result of the existence of the LDA or GGA, which are excellent approximations for weakly correlated systems. Here we will argue that the usefulness of SDFT is the existence of the local dynamical mean field approximations.
Any of the discussed functionals can be obtained by a Legendre transform of the effective action. The electronic Green function of a system can be obtained by probing the system by a source field and monitoring the response. To obtain Γ BK [G] we probe the system with time-dependent two-variable source field J(x, x ′ ) or its imaginary frequency transform J(r, r ′ , iω) defined in all space. If we restrict our consideration to saddle point solutions periodic on the original lattice, we can assume that the field obeys the periodicity criteria J(r + R, r ′ + R, iω) = J(r, r ′ , iω). This restricts the electronic Green function to be invariant under lattice translations. In order to obtain a theory based on the density as a physical variable, we probe the system with a static periodical field J(r)δ(τ ). This delivers 45,46,47 the density functional theory Γ DF T [ρ]. In order to obtain Γ SDF [G loc ] we will probe the system with a local field J loc (x, x ′ ) = J loc (r, r ′ , τ − τ ′ ) restricted by θ loc (r, r ′ ).
Introduction of the time dependent local source J loc (x, x ′ ) modifies the action of the system (3) as follows Due to translational invariance, the integral over r variable here is the same for any unit-cell Ω c and the integral over r ′ should be restricted by the area where J loc = 0, i.e. by the cluster area Ω loc .The average of the operator ψ(x)ψ + (x ′ ) probes the local Green function which is precisely defined by expression (2). The partition function Z, or equivalently the free energy of the system F, becomes a functional of the auxiliary source field The effective action for the local Green function, i.e., spectral density functional, is obtained as the Legendre transform of F with respect to the local Green function G loc (x, x ′ ), i.e.
where we use the compact notation TrJ loc G loc for the integrals Using the condition: J loc = −δΓ SDF /δG loc to eliminate J loc in (6) in favor of the local Green function we finally obtain the functional of the local Green function alone. The source field sets the degree of locality of the object of interest. Considering its definition by expanding the cluster till entire solid, we obtain the Baym-Kadanoff functional which determines the Green function in all space. Shrinking its definition to a singe point r and assuming its frequency (time) independence, i.e. J(r)δ(r − r ′ )δ(τ − τ ′ ), we obtain density functional theory. In its extremum, all functionals always reach the total free energy of the system regardless the choice of the variable. This situation is similar 46 to classical thermodynamics where the thermodynamic potential is either the Helmholtz free energy, or the Gibs free energy or the entalpy depending on which variables, temperature, pressure, volume are used. Note also that due to assumed time-dependence of the source field, away from the extremum the Green function functionals cannot be interpreted as energies.
Having repeated a formal derivation of the existence 29 of the functional Γ SDF [G loc ] as well as of the functionals Γ BK [G] and Γ DF T [ρ] we now come to the problem of writing separately various contributions to it. This development parallels the well known decomposition of the total energy into kinetic energy of a non interacting system, potential energy, Hartree energy and exchange-correlation energy. The strategy consists in performing an expansion of the functional in powers of the charge of the electron 29,45,47,48,49 . The lowest order term is the kinetic part of the action, and the energy associated with the external potential V ext . In the Baym Kadanoff Green function theory this term has the form (3): The G 0 (r, r ′ , iω) is the non-interacting Green function, which is given by where µ is a chemical potential. Note that since finite temperature formulation is adopted we did not obtain simply K BK [G] = Tr( − ∇ 2 + V ext )G but also have got all entropy based contributions. Let us now turn to the density functional theory. In principle, it does not have a closed formula to describe fully interacting kinetic energy as the density functional. However, it solves this problem by introducing a non-interacting part of the kinetic energy. It is described by its own Green function G KS (r, r ′ , iω), which is related to the Kohn-Sham (KS) representation. An auxiliary set of non-interacting particles is introduced which is used to mimic the density of the system. These particles move in some effective one-particle Kohn-Sham potential V ef f (r) = V ext (r) + V int (r). This potential is chosen merely to reproduce the density and does not have any other physical meaning at this point. The Kohn-Sham Green function is defined in the entire space by the relation is adjusted so that the density of the system ρ(r) can be found from G KS (r, r ′ , iω). Since the exact Green function G and the local Green function G loc can be also used to find the density, we can write a general relationship: where the sum over iω assumes the summation on the Matsubara axis at given temperature T . With the introduction of G KS the non-interacting kinetic portion of the action plus the energy related to V ext can be written in complete analogy with (8) as follows In order to describe the different contributions to the thermodynamical potential in the spectral density functional theory, we introduce a notion of the energy-dependent analog of Kohn-Sham representation. These auxiliary particles are interacting so that they will describe not only the density but also a local part of the Green function of the system, and will feel a frequency dependent potential. The latter is a field described by some effective mass operator M ef f (r, r ′ , iω) = V ext (r)δ(r − r ′ ) + M int (r, r ′ , iω). We now introduce an auxiliary Green function G(r, r ′ , iω) connected to our new "interacting Kohn-Sham" particles so that it is defined in the entire space by the relationship G −1 (r, r ′ , iω) = G −1 0 (r, r ′ , iω) − M int (r, r ′ , iω). Thus, M int (r, r ′ , iω) is a function which has the same range as the source that we introduce: it is adjusted until the auxiliary G(r, r ′ , iω) coincides with the local Green function inside the area restricted by θ loc (r, r ′ ), i.e We illustrate the relationship between all introduced Green functions in Fig. 2. Note that G(r, r ′ , iω) also delivers the exact density of the system. With the help of G the kinetic term in the spectral density functional theory can be represented as follows Since G KS is a functional of ρ, DFT considers the density functional as the functional of Kohn-Sham wave functions, i.e. as Γ DF T [G KS ]. Similarly, since G is a functional of G loc , it is very useful to view the spectral density functional Γ SDF as a functional of G: where the unknown interaction part of the free energy Φ SDF [G loc ] is the functional of G loc . If the Hartree term is explicitly extracted, this functional can be represented as where E H [ρ] is the Hartree energy depending only on the density of the system, and where Φ xc SDF [G loc ] is the exchangecorrelation part of the free energy. Notice that the density of the system can be obtained via G loc or G, therefore the Hartree term can be also viewed as a functional of G loc or G. Notice also, that since the kinetic energies (8), (12), (14) This leads to the equations for the corresponding Green functions in all theories: as well as By using (9) for G −1 0 and by multiplying both parts by the corresponding Green functions we obtain familiar Dyson equations and The stationarity condition brings the definition of the auxiliary mass operator M int (r, r ′ , iω) which is the variational derivative of the interaction free energy with respect to the local Green function: It plays the role of the effective self-energy which is short-ranged (local) in the space. The corresponding expressions hold for the interaction parts of the exact self-energy of the electron Σ int (r, r ′ , iω) and for the interaction part of the Kohn-Sham potential V int (r).
If the external potential is added to these quantities we obtain total effective self-energies/potentials of the SDF, BK and DF theories: M ef f (r, r ′ , iω), Σ ef f (r, r ′ , iω), V ef f (r) respectively. If the Hartree potential V H (r) is separated we obtain the exchange-correlation parts: M xc (r, r ′ , iω), Σ xc (r, r ′ , iω), V xc (r). Note that strictly speaking the substitution of variables, G KS vs. ρ, in the density functional as well as the substitution of variables, G vs. G loc , in the spectral density functional is only possible under the assumption of the so-called V -representability (or M -representability), i.e. the existence of such effective potential (mass operator) which can be used to construct the exact density (local Green function) of the system via the non-interacting Kohn-Sham particles of the DFT or its energy-dependent generalization in SDFT.
Note also that the effective mass-operator of spectral density functional theory is local by construction, i.e. it is nonzero only within the cluster area Ω loc restricted by θ loc (r, r ′ ). It is an auxiliary object which cannot be identified with the exact self-energy of the electron Σ ef f (r, r ′ , iω). This is similar to the observation that the Kohn-Sham potential of the DFT cannot be associated with the exact self-energy as well. Nevertheless, the SDFT always delivers local Green function and the total free energy exactly (at least in principle) as long as the exact functional is used. In the limit when the exact self-energy of the electron is indeed localized within Ω loc , the SDFT becomes the Baym-Kadanoff functional which delivers the full Green function of the system, i.e. we can immediately identify M ef f (r, r ′ , iω) with Σ ef f (r, r ′ , iω) and the poles of G(r, r ′ , iω) with exact poles of G(r, r ′ , iω) where the information about both k and energy dependence as well as life time of the quasiparticles is contained. We thus see that, at least formally, increasing the size of Ω loc in the SDF theory leads to a complete description of the many-body system, the situation quite different from the DFT which misses such scaling.
From a conceptual point of view, the spectral density functional approach constitutes a radical departure from the DFT philosophy. The saddle-point equation (23) is the equation for a continuous distribution of spectral weight and the obtained local spectral function G loc can now be identified with the observable local (roughly speaking, k-integrated) one-electron spectrum. This is very different from the Kohn-Sham quasiparticles which are the poles of G KS not identifiable rigorously with any one-electron excitations. While the SDFT approach is computationally more demanding than DFT, it is formulated in terms of observables and gives more information than DFT.
On one side, spectral density functional can be viewed as approximation or truncation of the full Baym Kadanoff theory where Φ BK [G] is approximated by Φ SDF [G loc ] by restricting G to G loc 30 and the kinetic functionals K BK [G] and K SDF [G] are thought to be the same. Such restriction will automatically generate a short-ranged self-energy in the theory. This is similar to the interpretation of DFT as approximation which would generate the DFT potential as the self-energy. However, SDFT can be thought as a separate theory whose manifestly local self-energy is an auxiliary operator introduced to reproduce the local part of the Green function of the system, exactly like the Kohn-Sham ground state potential is an auxiliary operator introduced to reproduce the density of the electrons in DFT.
Spectral density functional theory contains the exchange-correlation functional Φ SDF [G loc ]. An explicit expression for it involving a coupling constant λ = e 2 integration can be obtained in complete analogy with the Harris-Jones formula 50 of density functional theory 49 . One considers Γ SDF [G, λ] at an arbitrary interaction λ and expresses Here the first term is simply the kinetic part K SDF [G] as given by (14) which does not depend on λ. The second part is thus the unknown functional Φ SDF [G loc ]. The derivative with respect to the coupling constant in (3) is given by the is the density-density correlation function at a given interaction strength λ computed in the presence of a source which is λ dependent and chosen so that the local Greens function of the system is G. Since ψ + (x)ψ(x) = ρ(r)δ(τ ), we can obtain : Establishing the diagrammatic rules for the functional Φ SDF [G loc ] while possible 29 , is not as simple as for the functional Φ BK [G]. The latter is formally represented as a sum of two-particle diagrams constructed with G and v C . It is well known that instead of expanding Φ BK [G] in powers of the bare interaction v C and G, the functional form can be obtained by introducing the dynamically screened Coulomb interaction W (r, r ′ , iω) as a variable. In the effective action formalism 30 this was done by introducing an auxiliary Bose variable coupled to the density, which transforms the original problem into a problem of electrons interacting with the Bose field. W is the connected correlation function of the Bose field.
Our effective action is now a functional of G, W and of the expectation value of the Bose field. Since the latter couples linearly to the density it can be eliminated exactly, a step which generates the Hartree term. After this elimination, the functional takes the form The entire theory is viewed as the functional of both G and W. Here, Ψ BK [G, W ] is the sum of all two-particle diagrams constructed with G and W with the exclusion of the Hartree term, which is evaluated with the bare Coulomb interaction. An additional stationarity condition δΓ BK /δW = 0 leads to the equation for the screened Coulomb interaction W where the function Π(r, r ′ , iω) = −2δΨ BK /δW (r, r ′ , iω) is the exact interacting susceptibility of the system, which is already discussed in connection with representation (30). A similar theory is developed for the local quantities 30 , and this generalization represents the ideas of extended dynamical mean field theory 35 , now viewed as an exact theory. Namely, one constructs an exact functional of the local Greens function and the local correlator of the Bose field coupled to the density which can be identified with the local part of the dynamically screened interaction. The real-space definition of it is the following which is non-zero within a given cluster Ω loc . Note that formally this cluster can be different from the one considered to define the local Green function (2) but we will not distinguish between them for simplicity. An auxiliary interaction W(r, r ′ , iω) is introduced which is the same as the local part of the exact interaction within non-zero area of θ loc (r, r ′ ) The interaction part of the spectral density functional is represented in the form similar to (32) and the spectral density functional is viewed as a functional Γ SDF [G loc , W loc ] or alternatively as a functional formally not a sum of two-particle diagrams constructed with G loc and W loc , but in principle a more complicated diagrammatic expression can be derived. Alternatively, a more explicit expression involving a coupling constant integration can be given. Examining stationarity δΓ SDF /δW = 0 yields a saddle-point equation for W(r, r ′ , iω) where the effective susceptibility of the system is the variational derivative Notice again a set of parallel observations for P as for M ef f , Eq. (26). The effective susceptibility of spectral density functional theory is local by construction, i.e. it is non-zero only within the cluster restricted by θ loc (r, r ′ ). Formally, it is an auxiliary object and cannot be identified with the exact susceptibility of the electronic system Π(r, r ′ , iω). However, if the exact susceptibility Π(r, r ′ , iω) is sufficiently localized, this identification becomes possible. If cluster Ω loc includes physical area of localization, we can immediately identify P(r, r ′ , iω) with Π(r, r ′ , iω) and W(r, r ′ , iω) with W (r, r ′ , iω) in all space. However, both W and W are always the same within Ω loc regardless its size, as it is seen from (34) and (35). At the stationarity point, Γ SDF [G, W] is the free energy F of the system. If one inserts (20) into (14) and (37) into (36) we obtain the formula: Similar formulae hold for the Baym-Kadanoff and density functional theories where the first two terms in all expressions (39), (40), (41) are interpreted as corresponding kinetic energies, the third term is the energy related to the external potential V ext which is in fact TrV ext ρ in all cases. The other terms represent the interaction parts of the free energy. Note that all entropy originated contributions are included in both kinetic and interaction parts. If temperature goes to zero, the entropy part disappears and the total energy formulae will be recovered. For example, in spectral density functional theory we obtain: We will also discuss this limit later in more details in Section III. The SDFT approach is so far not very useful since a tractable expression for the functional form of Φ SDF [G loc ] or Ψ SDF [G loc , W loc ] has not been given yet. This is quite similar to the unknown exchange-correlation functional of the DFT. As we have learned from the developments of the dynamical mean-field methods, a very useful approximation exists to access these functionals. This is based on a full many-body solution of a finite-size cluster problem treated as an impurity embedded into a bath subjected to a self-consistency condition. Such local dynamical mean field theory will be discussed below.

B. Local Dynamical Mean Field Approximation
The spectral density functional theory, where an exact functional of certain local quantities is constructed in the spirit of Ref. 29 uses effective self-energies and susceptibilities which are local by construction. This property can be exploited to find good approximations to the interaction energy functional. For example, if it is a priori known that the real electronic self-energy is local in a certain portion of the Hilbert space, a good approximation is the corresponding local dynamical mean field theory obtained for example by a restriction or truncation of the full Baym-Kadanoff functional or its generalization to use W and G as natural variables, to local quantities in the spirit of Ref. 30.
The local DMFT approximates the functional Φ SDF [G loc ] (or Ψ SDF [G loc , W loc ]) by the sum of all two-particle diagrams evaluated with G loc and the bare Coulomb interaction v C (or screened local interaction W loc ). In other words, the functional dependence of the interaction part Φ BK [G] in the Baym-Kadanoff functional for which the diagrammatic rules exist is now restricted by G loc and is used as Obviously that the variational derivative of such restricted functional will generate the self-energy confined in the same area as the local Green function itself.
Remarkably the summation over all local diagrams can be performed exactly via introduction of an auxiliary quantum impurity model subjected to a self-consistency condition 10,51 . If this impurity is considered as a cluster C, the cellular DMFT (C-DMFT) can be used which breaks the translational invariance of the lattice to obtain accurate estimates of the self energies. The C-DMFT approximation, can also be motivated using the cavity construction. The solid should be separated onto large cells which circumscribe the areas Ω loc . Considering the effective action S, Eq. (3), the integration volume is separated onto the cellular area Ω C and the rest bath area Ω − Ω C = Ω bath . The action is now represented as the action of the cluster cell, Ω C plus the action of the bath, Ω bath , plus the interaction between those two. We are interested in the local effective action S C of the cluster degrees of freedom only, which is obtained conceptually by integrating out the bath in the functional integral: where Z C and Z are the corresponding partition functions. This integration is carried out approximately, keeping only a charge-charge interaction as quartic terms and neglecting all the higher order terms generated in this process to arrive to a cavity action of the form 18,26,30,35 : where the integration over the spatial variables is performed over Ω C . Here G 0 (x, x ′ ) or its Fourier transform G 0 (r, r ′ , iω) is identified as the bath Green function appeared in the Dyson equation for the local mass operator M int (r, r ′ , iω) and for the local Green function G loc (r, r ′ , iω) of the cluster, and V 0 (r, r ′ , iω) is the "bath interaction" appeared in the Dyson equation for the local susceptibility P(r, r ′ , iω) and local interaction W loc (r, r ′ , iω), i.e Note that neither G 0 nor V 0 can be associated with non-interacting G 0 and bare interaction v C , respectively. Note also that both r and r ′ indexes in G 0 (r, r ′ , iω) and in V 0 (r, r ′ , iω) vary within the cellular area Ω C . The same should be assumed for the local quantities G loc (r, r ′ , iω) and W loc (r, r ′ , iω). Since these functions are translationally invariant on the original lattice, this property can be used to set up these functions within Ω C . An interesting observation can be made on the role of the impurity model which in the present context appeared as an approximate way to extract the self-energy of the lattice using input bath Green function and bath interaction. Alternatively, the impurity problem can be thought itself as the model which delivers exact mass operator of the spectral density functional 29 . If the latter is known, there should exist such bath Green function and such bath interaction which can be used to reproduce it. In this respect, the local interaction W loc appeared in our formulation can be thought as an exact way to define the local Coulomb repulsion "U ", i.e. such interaction which delivers exact local self-energy.
To summarize, the effective action for the cluster cell (44), the Dyson equations (45), (46) connecting local and bath quantities as well as the original Dyson equations (20), (37) constitute a self-consistent set of equations as the saddle-point conditions extremizing the spectral density functional Γ SDF (G, W). They combine cellular and extended versions of DMFT and represent our philosophy in the ab initio simulation of a strongly correlated system. Since M int and P are unknown at the beginning, the solution of these equations assumes self-consistency. First, assuming some initial M int ,and P the original Dyson equations (20), (37) are used to find Green function G and screened interaction W. Second, the Dyson equations for the local quantities (45), (46) are used to find G 0 , V 0 . Third, quantum impurity model with the cluster action S loc after (44) is solved by available many-body technique to give new local M int and P. The process is repeated till self-consistency is reached. This is schematically illustrated in Fig. 3. Note here that while single-site impurity problem has a well-defined algorithm to extract the lattice self-energy, this is not generally true for the cluster impurity models 19 . The latter provides the self-energy of the cluster, and an additional prescription such as implemented within cellular DMFT or using DCA should be given to construct the self-energy of the lattice. Unfortunately, writing down the precise functional form for Φ SDF [G loc , W loc ] or Ψ SDF [G loc , W loc ] is still a problem because the evaluation of the entropy requires the evaluation of the energy as a function of temperature and an additional integration over it. In general, the free energy F = E − T S, where E is the total energy and S is the entropy. Since Γ SDF [G] = K SDF [G] + Φ SDF [G loc ], both energy and entropy terms exist in the kinetic and interaction functionals. The energy part of K SDF [G] = Tr( − ∇ 2 + V ext )G and the energy part of Φ SDF [G loc , W loc ] can be written explicitly as 1 2 TrM int G loc . The entropy correction is a more difficult one. In principle, it can be evaluated by performing calculations of the total energy E SDF [G] = Tr( − ∇ 2 + V ext )G + 1 2 TrM int G loc at several temperatures and then taking the integral 10 : The infinite temperature limit S(∞) for a well defined model Hamiltonian can be worked out. This program was implemented for the Hubbard model 52 and for Ce 31 .
Two well separate problems are now seen. For a given material using the formulae (20), (37), (2), (34), (45), (46) G, W, G loc , W loc , G 0 , V 0 should be computed using the methods and algorithms of the electronic structure theory. This procedure will in part be described in Section III. With given input G 0 and V 0 , the solution of the impurity model constitutes a well separated problem which can be carried out either using the QMC method or other impurity solver. Some of the techniques are discussed in Refs. 10,11. In Section IV, while applying a simplified version of the theory to plutonium, we will briefly describe an impurity solver used in that calculation. A full description of this method will appear elsewhere 41 .
The described algorithm is quite general, totally ab initio and allows to determine all quantities, such as the oneelectron local Green functions G loc and the dynamically screened local interactions W loc . Unfortunately, its full implementation is a very challenging project which so far has only been carried out at the level of model Hamiltonians 26 . There are several simplifications which can be made, however. The screened Coulomb interaction W(r, r ′ , iω) can be treated on different levels of approximations. In many cases used in practical calculations with the LDA+DMFT method, this interaction W is assumed to be static and parametrized by a set of some optimally screened on-site parameters, such as Hubbard U and exchange J. These parameters can be fixed by constrained density functional calculations, extracted from atomic spectra data or adjusted to fit the experiment. Since the described theory can perform a search in a constrained space with fixed interaction W, this justifies the use of U and J as input numbers. A more refined approximation, can use a method such as GW to generate an energy-dependent W 53 which is then treated using extended DMFT 26 . Alternatively we can envision that W is already so short ranged that we can ignore the EDMFT self consistency condition, and we treat W as W f ix (x, x ′ ) .This leads to performing a partial self-consistency with respect to the Green function only. The procedure is reduced to solving Dyson equations (20), (45) as well as to finding M int via the solution of the impurity problem. A full self-consistency can finally be restored by including a second loop to relax W.
A methodological comment should be made in order to make contact with the literature of cluster extensions of single site DMFT within model Hamiltonians. We adopted a less restrictive notion of locality by defining an effective action of the one particle Green function (and of the interaction) whose arguments are in nearby unit cells. This maintains the full translation invariance of the lattice. At the level of the exact effective action , this is an exact construction, and its extremization will lead to portions of the exact Greens function which obeys causality. Notice however that it has been proved recently 19 that generating approximations to the exact functional by restricting the Baym Kadanoff functional to non local Greens functions leads to violations of causality. For this reason, we propose to use techniques such as CDMFT which are manifestly causal for the purpose of realizing approximations to the local Greens functions.
Our final general comment concerns the optimal choice of local representation or, precisely, the definition of the local Green function. This is because the local dynamical mean-field approximation is likely to be accurate only if we know in which portion of the Hilbert space the real electronic self-energy is well localized. Unfortunately, this is not known a priori, and in principle, only a full cluster DMFT calculation is capable to provide us some hints in attempts to answer this question. However, considerable empirical evidence can be used as a guide for choosing a basis for DMFT calculations, and we discuss these issues in the following sections.

C. Choice of Local Representation
We have already pointed out that spectral density functional theory is a basis set dependent theory since it probes the Green function locally in a certain region determined by a choice of basis functions in the Hilbert space. Provided the calculation is exact, the free energy of the system and the local spectral density in that Hilbert space will be recovered regardless the choice of it. We have developed the theory assuming that the basis in the Hilbert space is indeed the real space which gives us the choice (2) for the local Green function, i.e. the part of the real Green function restricted by θ loc (r, r ′ ). While this is most natural choice for the purpose of formulating locality in r and r ′ variables, it is also very useful to discuss a more general choice, connected to some space of orbitals χ ξ (r) which can be used to represent all the relevant quantities in our calculation. As we have in mind to utilize sophisticated basis sets of modern electronic structure calculations, we will sometimes waive the orthogonality condition and will introduce the overlap matrix O ξξ ′ = χ ξ |χ ξ ′ especially in cases when we discuss a practical implementation of the method.
We note that the space χ ξ (r) can in principle be interpreted as the reciprocal space plane wave representation χ ξ (r) = e i(k+G)r , ξ = k + G with k being the Brillouin zone vector and G being the reciprocal lattice vector. Thus the Green function can be probed in the region of the reciprocal space. It can be interpreted as the real space representation if χ ξ (r) = δ(ξ − r) where the sums over ξ are interpreted as the integrals over the volume, and the locality in this basis is precisely exploited in (2). A tremendous transparency of the theory will also arrive if we interpret the orbital space {χ ξ } as a general non-orthogonal tight-binding basis set when index ξ combines the angular momentum index lm, and the unit cell index R, i.e. χ ξ (r) = χ lm (r − R) = χ α (r − R). Note that we can add additional degrees of freedom to the index α such, for example, as multiple kappa basis sets of the linear muffin-tin orbital based methods, Gaussian decay constants in the Gaussian orbital based methods, and so on. If more than one atom per unit cell is considered, index α should be supplemented by the atomic basis position within the unit cell, which is currently omitted for simplicity. For spin unrestricted calculations α accumulates the spin index σ and the orbital space is extended to account for the eigenvectors of the Pauli matrix.
Let us now introduce the representation for the exact Green function in the localized orbital representation Assuming the single-site impurity case, we can separate local and non-local parts G loc (r, r ′ , iω) + G non−loc (r, r ′ , iω) as follows where we denoted the site-diagonal matrix elements δ RR ′ G αβ (R − R ′ , iω) as G loc,αβ (iω). Note that this definition is different from the real-space definition (2). For example, (2) contains the information about the density of the system. The formula (49) does not describe the density since R = R ′ elements of the matrix G αβ (R − R ′ , iω) are thrown away. The locality of (49) is controlled exclusively by the decay of the orbitals χ α (r) as a function of r, not by θ loc (r, r ′ ) The local part of the Green function, G loc (r, r ′ , iω), which is just defined with respect to the Hilbert space {χ α } can be found by developing the corresponding spectral density functional theory. Since the basis set is assumed to be fixed, the matrix elements G loc,αβ (iω) appear only as variables of the functional. As above, we introduce an auxiliary Green function G αβ (k, iω) to deal with kinetic energy counterpart. Stationarity yields the matrix equation: where the non-interacting Green function (9) is the matrix of non-interacting one-electron Hamiltonian The self-energy M int,αβ (iω) is the derivative δΦ SDF [G loc,αβ (iω)]/δG loc,αβ (iω) and takes automatically the kindependent form. While formally exact, this theory would have at least one undesired feature since, for example, the density of the system can no longer be found from the definition (49) of G loc (r, r ′ , iω). As a result the Hartree energy cannot be simply recovered. If treated exactly Φ SDF [G loc,αβ (iω)] should contain the Hartree part. However, we see that the theory delivers k-independent M int,αβ (iω) including the Hartree term. There seems to be a paradox since modern electronic structure methods calculate the matrix element of the Hartree potential within a given basis exactly, i.e. χ k α |V H |χ k β . The k-dependence is trivial here and is connected to the known k-dependence of the basis functions used in the calculation. Therefore, while formulating the spectral density functional theory for electronic structure calculation, we need to keep in mind that in many cases, the k-dependence is factorizable and can be brought into the theory without a problem. This warns us that the choice of the local Green function has to be done with care so that useful approximations to the functional can be worked out. It also shows that in many cases the k-dependence is encoded into the orbitals. It is not that non-trivial k-dependence of the self-energy operator, which is connected to the fact that M int (r, r ′ , iω) may be long-range, i.e. decay slowly when r departs from r ′ . It may very well be proportional to δ(r − r ′ ) like the LDA potential and still deliver the k-dependence.
It turn out that the desired k dependence with the choice of the Green function after (49) can be quickly reinstated if we add the density of the system as another variable to the functional. This is clear since the density is a particular case of the local Green function in (2 ) taken at r = r ′ and summed over iω. Therefore combination of definition (49) and ρ is another, third possibility of defining G loc . This will allow treatment of all local Hartree-like potentials without a problem. Moreover, as we discuss below, this may allow to design better approximations to the functional since the Hilbert space treatment of locality is more powerful: it may allow us to treat more long-ranged self-energies than the ones restricted by θ loc (r, r ′ ), and the basis sets can be optimally adjusted to specific self-energies exactly as the basis sets used in electronic structure calculations are tailored to the LDA potential.
We have noted earlier that the mass operator M int (r, r ′ , iω) is an auxiliary object of the spectral density functional theory. It has the same meaning as the DFT Kohn-Sham potential: it is local self-mass operator that needs to be added to the non-interacting Green function in order to reproduce the local Green function of the system, as the DFT potential is added to the non-interacting Green function to reproduce the density of the system. Roughly speaking, SDFT provides the exact energy and exact one-electron density of states which is advantageous compared to the DFT which provides the energy and the density only. However, we obtain the full k-dependent one-particle spectra as the poles of auxiliary Green function G(r, r ′ , z). Can these poles be interpreted as the exact k-dependent one-electron excitations? This question is similar to the question of the DFT: can the Kohn-Sham spectra be interpreted as the physical one-electron excitations? To answer both questions we need to know something about exact self-energy of the electron. If it is energy-independent, totally local, i.e. proportional to δ(r − r ′ ) and well-approximated by the DFT potential, the Kohn-Sham spectra represent real one-electron excitations. The exact SDFT waives most of the restrictions: if the real self-energy is localized within the area R loc , the exact SDFT calculation with the cluster Ω loc including R loc will find the exact k-dependent spectrum. If we pick Ω loc larger than R loc , the SDFT equations themselves will choose physical localization area for the self-energy during our self-consistent calculation. However, these statements become approximate if we utilize the local dynamical mean field approximation instead of extremizing the exact functional. Even if the real self-energy of the electron is sufficiently short-ranged, this approximation will introduce some error in the calculation, the situation similar to LDA within DFT. However, the local dynamical mean field theory does not necessarily have to be formulated in real space. The assumption of localization for self-energy can be done in some portion of the Hilbert space. In that portion of the Hilbert space the cluster impurity model needs to be solved.
The choice of the appropriate Hilbert space, such, e.g., as atomic-like tight-binding basis set is crucial, to obtain an economical solution of the impurity model. Let us for simplicity discuss the problem of optimal basis in some orthogonal tight-binding (Wannier-like) representation for the electronic self-energy We can separate our orbital space {χ α } onto the subsets describing light {χ A } and heavy {χ a } electrons. Assuming either off-diagonal terms between them are small or we work with exact Wannier functions, the self-energy Σ(r, r ′ , iω) can be separated onto contributions from the light, Σ L (r, r ′ , iω), and from the heavy, Σ H (r, r ′ , iω), electrons. Σ αβ (k, iω) is expected to be k-dependent but largely ω independent for the light block, i.e Σ L (r, The k-dependency here should be well-described by LDA-like approximations, therefore we expect Σ L (r, r ′ , iω) ∼ V ef f (r)δ(r − r ′ ). A different situation is expected for the heavy block where we would rely on the result The first term here gives the k-dependence coming from an LDA-like potential. It describes the dispersion in the heavy band. The second term is the energy dependent correction where site-diagonal approximation R = R ′ is imposed. What is the best choice of the basis to use in connection with Σ ′ ab (iω) in (53)? Here the decay of the orbitals χ α (r) as a function of r is now entirely in charge of the self-energy range. In light of the spectral density functional theory, the answer is the following: the local dynamical mean field approximation would work best for such basis χ α (r) whose range approximately corresponds to a self-energy localization R loc of the real electron. Even though R loc is not known a priori, something can be learned about its value based on a substantial empirical evidence. It is, for example, known that LDA energy bands when comparing to experiments at first place miss the energy dependent Σ ′ ab (iω) like corrections. This is the case for bandwidths in transition metals (and also in simple metals), the energy gaps of semiconductors, etc. It is also known that many-body based theories work best for massively downfolded model Hamiltonians where only active low-energy degrees of freedom at the region around the Fermi level E F remain. The many-body Hamiltonian is obtained as a fit to the bands near E F . This can always be done by long ranged Wannier functions. It is also clear that the correlation effects are important at first place for the partially occupied bands since only these bring various configurational interactions in the many-body electronic wave functions. For example, the well-known oneband Hamiltonian for CuO 2 plane of high-T c materials considers an antibonding combination of Cu x 2 −y 2 and O x,y orbitals which crosses E F . Also, the calculations based on the LDA+DMFT method usually obtain reliable results when treating only the bands crossing the Fermi level as the correlated one-electron states. This is, for example, the case of Pu or our 25 and previous 54 calculation for LaTiO 3 where t 2g three band Hamiltonian is considered. All this implies that the range for Σ ′ ab (iω) term in (53) should correspond to the properly constructed Wannier orbitals, which is fairly long-ranged. What happen if we instead utilize mostly localized representation which, for example, can be achieved by tight-binding fits to the energy bands at higher energy scale. For the case of CuO 2 this would correspond to a three band Hamiltonian with Cu x 2 −y 2 and O x,y orbitals treated separately. For LaTiO 3 system this is a Hamiltonian derived from Ti t2g and O p orbitals. The answer here can be given as a practical matter of most economic way to solve the impurity problem: provided Cu and O levels are well separated, provided both approaches use properly downfolded for each case Coulomb interaction matrix elements V RR ′ R"R ′′′ αβγδ , and provided correlations are treated on all orbitals, the final answer should be similar regardless the choice of the basis. A faster algorithm will be obtained by treating the one-band Hamiltonian with antibonding Cu x 2 −y 2 -O x,y orbital. If indeed the self-energy is localized on the scale of the distance between Cu and O, it is clear where the inefficiency of the three-band model appears: the second term in (53) needs to be extended within the cluster involving both Cu and nearest O sites and should involve both Cu and O centered orbitals simply to reach the cluster boundary. In the one-band case this is encoded into the decay of the properly constructed Wannier state.
The previous discussion is merely a conjecture. It does not imply that the localization range for the real self-energy of correlated electron at given frequency ω is directly proportional to the size of Wannier states located in the vicinity of ω + µ. It may very well be that in many cases this range is restricted by a single atom only (atomic sphere of Cu in the example above). Clearly more experience can be gained by studying a correlation between the decay of the Coulomb matrix element V RRR ′ R ′ as a function of R − R ′ and the obtained matrix Σ(R − R ′ , ω) using a suitable cluster DMFT technique. These works are currently being performed and will be reported elsewhere 55 . The given discussion however warns that in general the best choice of the basis for single-site dynamical mean field treatment may not be the case of mostly localized representation. In this respect the area restricted by θ loc (r, r ′ ) which is used to formulate SDFT in the real space may need to be extended up to a cluster. However, alternative formulation with the choice of local Green function after (49) may be more economical since a single-site approximation may still deliver good results. As we have argued, such spectral density functional theory will also need the density of the system to complete the definition of local Green function. The local dynamical mean field approximation can be applied to the interaction functional Φ SDF which is viewed as Φ SDF [ρ, G loc ]. This idea is used by the LDA+DMFT method described below.

D. LDA+DMFT Method
Various methods such as LDA+U 7 , LDA+DMFT 11 and local GW 12,16 which appeared recently for realistic calculations of properties of strongly correlated materials can be naturally understood within spectral density functional theory. Let us, for example, explore the idea of expressing the energy as the density functional. Local density approximation prompts us that a large portion of the exchange-correlation part Φ xc [ρ] can be found easily. Indeed, the charge density is known to be accurately obtained by the LDA. Why not think of LDA as the most primitive impurity solver, which generates manifestly local self-energy with localization radius collapsed to a single r point? It is tempting to represent where the new functionalΦ SDF [G loc ] needs in fact to take care of those electrons which are strongly correlated and heavy, thus badly described by LDA. Conceptually, that means that the solution of the cluster impurity model for the light electrons is approximated by LDA and does not need a frequency resolution for their self-energies.
Unfortunately, the LDA has no diagrammatic representation, and it is difficult to separate the contributions from the light and heavy electrons. The E LDA xc [ρ] is a non-linear functional and it already includes the contribution to the energy from all orbitals in some average form. Therefore we need to take care of a non-trivial double counting, encoded in the functional Φ DC [G loc ]. The precise form of the double counting is related to the approximation imposed forΦ[G loc ]. We postpone this discussion until establishing the connection to the LDA+U method in the following subsection.
The LDA+DMFT approximation considers both the density and the local Green function G loc,αβ (iω) defined in (49) as the parameters of the spectral density functional 56 . A further approximation is made to accelerate the solution of a single-site impurity model: the functional dependence comes from the subblock of the correlated electrons only. If localized orbital representation {χ α } is utilized, a subspace of the heavy electrons {χ a } can be identified. Thus, the approximation can be written asΦ SDF [G loc,ab (iω)] , where G loc,ab (iω) is the heavy block of the local Green function. The double counting correction depends only on the average density of the heavy electrons. Its precise form will be discussed below, but for now we assume that Φ DC [G loc ] = Φ DC [n c ] withn c = T iω a G loc,aa (iω)e iω0 + , where index a runs within a correlated l c shell only. We can write the LDA+DFMT approximation for the interaction energy as follows: The kinetic energy part is treated as usual with introducing the auxiliary Green function G(r, r ′ , iω).
The full functional Γ LDA+DMF T [G] is considered as a functional of the matrix G αβ (R − R ′ , iω) or its Fourier transformed analog G αβ (k, iω). The stationarity is examined with respect to G αβ (k, iω) and produces the saddlepoint equation similar to (20). It has the following matrix form where the non-interacting Green function (9) is the matrix of non-interacting one-electron Hamiltonian The self-energy M int,αβ (k, iω) is the variational derivative of Φ LDA+DMF T [ρ, G loc ]. Its precise form depends on the basis set used in the LDA+DMFT calculation.
In general, it can be split onto several contributions including Hartree, LDA exchange-correlation, DMFT and the double-counting correction. In orthogonal tight-binding, both DMFT,M ab (iω), and double counting, δ ab V DC aa , matrices do not depend on k. These matrices are non-zero within the heavy block only. The Dyson equation (56) can be rewritten by separating from M int,αβ (k, iω) the total LDA potential V LDA (r) = V ext (r) + V H (r) + V LDA xc (r): The Green function G αβ (k, iω) obtained from (9) is used to find G loc,αβ (iω) = k G αβ (k, iω) which is then used in another Dyson equation to compute the bath Green function: In Section III we will also describe an accurate procedure to solve the real space form (45) of the Dyson equation using the LMTO basis set. The LDA+DMFT bath Green function G 0,ab (iω) is the only essential input to the auxiliary impurity model. Thus, the procedure of self-consistency within LDA+DMFT is reduced to the following steps. First, some self-energy matrix of the heavy orbitalsM ab (iω) is guessed. Then, the Dyson equation (56) is solved in the entire Hilbert space and delivers the Green function G αβ (k, iω). After that, the local Green function of the correlated electrons is constructed, which is then used in the equation (59) to deliver the bath Green function G 0,ab (iω). This matrix is the input to the impurity model. Solution of this model delivers the new self-energyM ab (iω) and the process is iterated towards self-consistency. Notice that once the DMFT self-consistency is reached, the process can either be stopped or continued since the Green function G αβ (k, iω) delivers new charge density of the system which modifies the Hartree and LDA exchangecorrelation potentials in the expression (58). In this respect, the LDA+DMFT method assumes a double iterational loop, the internal one over the self-energy and the external one over the density. This is precisely dictated by the spectral density functional stationarity condition. We illustrate such loop on Fig. 4. Note that in order to access accurate total energies and remove ambiguity that the LDA Green function (and not any other one) is used as an input to the DMFT calculation, this density self-consistency loop needs to be carried out. Our application to the volume expansion in Pu described later in this paper involves solution of the SDFT equations allowing the full relaxation of the charge density.
Iterations over the density are not complicated to incorporate in the programs for electronic structure calculations. The spherical part of the density at a given site can be written approximately using the atomic sphere approximation as an integral over the partial density of states N l (E) till the Fermi level E F where ϕ l (r, E) are the solutions of the radial Schrödinger equation. Usually these are constructed using spherical part of the LDA potential but in the present context the non-local self-energy operator needs to be utilized 60 : Expression (60) can be simplified further if we assume a Taylor expansion of ϕ l (r, E)=ϕ l (r, E νl ) + (E − E νl )φ l (r, E νl ) around some linearization energies E νl taken at the centers of gravities of the occupied energy bands, i.e.
During the iterations, DMFT modifies the density of states, ∆N l (E), which leads to the change ∆ρ(r) of the density. The latter has a feedback onto the change ∆ϕ l (r, E) [or changes ∆ϕ l (r, E νl ), ∆φ l (r, E νl )]. If we assume that these changes are small comparing to the original LDA values, we can work out a linear response relationship for the density and for the LDA potential Thus, to first order, these are the quantities which are iterated in the external density loop of the spectral density functional as shown on Fig. 4. The main physical point of the LDA+DMFT methodology is identification of a subset of the correlated orbitals {χ a } which is separated from the full Hilbert space {χ α }. In the case like Pu, this is the subset of f-electron orbitals. In other situations, this subset can be isolated based on physical grounds. If {χ a } is appropriately constructed Wannier representation, this subset may describe the bands crossing the Fermi level. We expect the dynamical self-energy corrections to appear at first place only within the subset {χ a }. However, changes in the electronic densities of states, ∆N l (E), will appear for all light and heavy electrons.
We did not discuss so far the relaxation of the screened Coulomb interaction W(r, r ′ , iω), which, in principle, needs to be done during the self-consistency in parallel to G(r, r ′ , iω). We stress that the short-range behavior appears only for the local effective susceptibility P(r, r ′ , iω) in (33) and means its fast decay when r departs from r ′ . Contrary, the function W(r, r ′ , iω) can be as long range as the bare Coulomb interaction if necessary. This is dictated by Eq. (33) and is similar to the relationship (20) between G and M int .
The locality assumption for P should simplify the self-consistency over W. This should be faster than the one employed in the full GW method which formally tries to compute full Π(r, r ′ , iω). In the language of local orbital representation {χ α }, this means computation of all matrix elements for Π αβ (R − R ′ , iω) or its Fourier transform Π αβ (k, iω) as compared to the site-diagonal (δ RR ′ ) or small cluster cases of SDFT. This will be discussed below in connection to the recently proposed 12,16 local version of the GW method.
So far we did not mention the problem of the optimal choice of the double counting corrections. This is discussed below in connection to the LDA+U method.

E. Double Counting and LDA+U Method
Historically, the LDA+U method has been introduced 57 as an extension of the local spin density approximation (LSDA) to treat the ordered phases of Mott insulating solids. In this respect it is a natural extension of LSDA. However, this method was first to recognize that a better energy functional can be constructed if not only the density but the density matrix of correlated orbitals is brought into the density functional. We have discussed the correlated subset {χ a } and local Green functions G loc,αβ (iω) in connection to the LDA+DMFT method. The density matrix n ab is related to the correlated subblock of the local Green function Therefore, the LDA+U method can be viewed as an approximation (Hartree-Fock approximation) to the spectral density functional within LDA+DMFT. The correct interaction energy among the correlated electrons can be written down explicitly using the Hartree-Fock approximation. In our language the LDA+DMFT interaction energy functional (55) is rewritten in the form where the functional formΦ U [n ab ] is known explicitly: Here, indexes a, b, c, d involve fixed angular momentum l of the heavy orbitals and run over magnetic m and spin σ quantum numbers. The on-site Coulomb interaction matrix U abcd is the on-site Coulomb interaction matrix element V RRRR α=aβ=bγ=cδ=d appeared in (54) which is again taken for the subblock of the heavy orbitals. Note that sometimes U abcd is defined as V RRRR α=aβ=cγ=bδ=d . The double counting term Φ DC [n ab ] needs to be introduced since both the L(S)DA and U terms account for the same interaction energy between the correlated orbitals. This includes in first place the Hartree part. However, the precise form of the double counting is unclear due to non-linear nature of the LDA exchange-correlation energy. In practice, it was proposed 7 that the form for Φ DC is ab (U abab − U abba )and wheren σ c = a∈lc n aa δ σ a σ ,n c =n ↑ c +n ↓ c . Some other forms of the double countings have also been discussed in Ref. 58.
The minimization of the functional Γ LDA+U [ρ, n ab ] is now performed. The self-energy correction in (58) appears as the orbital dependent correctionM ab − V DC ab : While the correction is static, it is best viewed as the Hartree-Fock approximation to the self-energy M ab (iω) within the LDA+DMFT method. Note that such interpretation allows us to utilize double counting forms within LDA+DMFT asM(r, r ′ , i∞) orM(r, r ′ , i0). Note also that the solution of the impurity problem collapses in the LDA+U method since the self-energy is known analytically by formula (68). From a practical point of view, despite the great success of the LDA+U theory in predicting materials properties of correlated solids 7 there are obvious problems with this approach when applied to metals or to systems where the orbital symmetries are not broken. They stem from the well-known deficiencies of the Hartree-Fock approximation. The most noticeable is that it only describes spectra of magnetically ordered systems which have Hubbard bands. We have however argued that a correct treatment of the electronic structure of strongly correlated systems has to treat both Hubbard bands and quasiparticle bands on the same footing. Another problem occurs in the paramagnetic phase of Mott insulators: in the absence of any broken symmetry the LDA+U method reduces to the LDA, and the gap collapses. In systems like NiO where the gap is of the order of eV, but the Neel temperature is a few hundred Kelvin, it is unphysical to assume that the gap and the magnetic ordering are related. For this reason the LDA+U predicts magnetic order in cases that it is not observed, as, e.g., in the case of Pu 59 .

F. Local GW Approximation
We now discuss the relaxation of the screened Coulomb interaction W(r, r ′ , iω) which appeared in the spectral density functional formulation of the problem. Both LDA+DMFT and LDA+U methods parametrize the interaction W with optimally screened set of parameters, such, e.g, as the matrix U abcd appeared in (66). This matrix is supposed to be given by an external calculation such, e.g., as the constrained LDA method 28 . To determine this interaction self consistently an additional self-consistency loop described by the equations (33) and (46) has to be switched on together with calculation of the local susceptibility P(r, r ′ , iω) by the impurity solver. This brings a truly self-consistent ab initio method without input parameters and the double counting problems.
A simplified version of this method has been recently proposed 12,16 which is known as a local version of the GW method (LGW). Within the spectral density functional theory, this approximation appears as approximation to the functional Ψ SDF [G loc , W loc ] taken in the form As a result, the susceptibility P(r, r ′ , iω) is approximated by the product of two local Green functions, i.e.P = −2δΨ LGW /δW loc = G loc G loc , and the exchange-correlation part of our mass operator is approximated by the local GW diagram, i.e. M xc = δΨ LGW /δG loc = −G loc W loc . Thus, the impurity model is solved and the procedure can be made self-consistent: For a given M int and P, the Dyson equations (20), (37) for G and W are solved. Then, the local quantities G loc, W loc are generated and used to find new M int and P thus avoiding the computation of the bath Green function G 0 after (45), and the interaction V, after (46).
Note that since the local GW approximation (70) is relatively cheap from computational point of view, its implementation on a cluster and for all orbitals should not be a problem. The results of the single-site approximation for the local quantities have been developed independently and reported in the literature. 16 . The cluster extension is currently being performed and the results will be reported elsewhere 60 .
Note finally that the local GW approximation is not the only one which can be implemented as the simplified impurity solver. For example, another popular approximation known as the fluctuational exchange approximation (FLEX) can be worked out along the same lines. Note also that the combination of the DMFT and full GW diagram has been recently proposed 12,27 and a simplified implementation for Ni 27 , and for a model Hamiltonian 26 have been carried out. This procedure incorporates full k-dependence of the self-energy known diagrammatically within GW together with the additional local DMFT diagrams.

III. CALCULATION OF LOCAL GREEN FUNCTION
The solution of the Dyson equations described in the previous section for a given strongly correlated material requires the calculation of the local Green function during the iterations towards self-consistency. This is very similar to the procedure in the density functional theory, when the charge density is computed. A big advantage of DFT is the use of Kohn-Sham orbitals which reduces the equation (22) for the Kohn-Sham Green function to a set of one-particle Schrödinger's like equations for the wave functions. As a result the kinetic energy contribution is calculated directly and the evaluation of the total energy of a solid is not a problem. Here, a similar algorithm will be described for the energy-dependent Dyson equation, the solution in terms of the linear-muffin-tin orbital basis set will be discussed, and the formula for evaluating the total energy will be given.

A. Energy Resolved One-Particle Representation
We introduced the auxiliary Green function G(r, r ′ , iω) to deal with the kinetic part of the action in SDFT. It satisfies to the Dyson equation (9). Let us now introduce the representation of generalized energy-dependent oneparticle states where the left ψ L kjω (r) and right ψ R kjω (r) states satisfy to the following Dyson equations: [we dropped the imaginary unit for simplicity in the notation ψ kjω (r) which shall be thought as a shortened version of ψ kj (r, iω)]. These equations should be considered as the eigenvalue problems with complex non-hermitian self-energy. As a result, the eigenvalues E kjω [a shortened form for E kj (iω)] being the same for both equations are complex in general. The explicit dependency on the frequency iω in both eigenvectors and eigenvalues comes from the self-energy. Note that left and right eigenfunctions are orthonormal and can be used to evaluate the charge density of a given system using the Matsubara sum and the integral over the k-space: where We have cast the notation of spectral density theory in a form similar to DFT. The function g kjω is the Green function in the orthogonal left/right representation which plays a role of a "frequency dependent occupation number". It needs to be pointed out that the frequency dependent energy bands E kjω represent an auxiliary set of complex eigenvalues. These are not the true poles of the exact one-electron Green function G(r, r ′ , z) considered at complex z plane. However, they are designed to reproduce the local spectral density of the system. Note also that these bands E kjz are not the true poles of the auxiliary Green function G(r, r ′ , z). The latter ones still need to be located by solving a non-linear equation corresponding to the singularities in the expression (71) after analytic continuation to real frequency. For a one-band case this equation is simply: z + µ− E kz = 0 , whose solution delivers the quasiparticle dispersion Z k . General knowledge of the poles positions Z kj will allow us to write an alternative expression for G which is similar to (71), but with the eigenvectors found at Z kj thus carrying out no auxiliary frequency dependence. These poles are the real one-electron excitational spectra in case G is a good approximation to G. However, the use of (71) is advantageous, since it avoids additional search of poles and allows direct evaluation of the local spectral and charge densities the system.
The energy-dependent representation allows us to obtain a very compact expression for the total energy. As we have argued, the entropy terms are more difficult to evaluate. However, they are generally small as long as we stay at low temperatures. The pure kinetic part of the free energy expressed via [see, Eq.(39)] needs to be separated onto the energy and entropy terms. Both contributions can be evaluated without a problem, but in light of neglecting the entropy correction in the interaction part, we concentrate on evaluating the kinetic energy only: The SDFT total energy formula is now arrived by utilizing the relationship If the self-energy is considered as input to the iteration while the Green function is the output, near stationary point, it should have a convergency faster than the convergency in the Green function.
It is instructive to consider the non-interactive limit when the self-energy represents a local energy-independent potential, say, the ground-state Kohn Sham potential of the density functional theory. This provides an important test for our many-body calculation. It is trivial to see that in the DFT limit, we obtain the Kohn-Sham eigenfunctions and the one-electron energy bands are no longer frequency dependent. The sum over Matsubara frequencies in the expression for the charge density (76) can be performed analytically using the expression for the Fermi-Diraq occupation numbers: and the formula (76) collapses to the standard expression for the density of non-interacting fermions. The total energy expression (80) is converted back to the DFT expression for the total energy since the eigenvalue E kjω becomes the DFT band structure E kj , and the summation over Matsubara frequencies T iω e iω0 + g kjω gives according to (84) the Fermi-Diraq occupation number f kj . The standard DFT expression is recovered: where E kj = ψ kj | − ∇ 2 + V ef f |ψ kj = ψ kj | − ∇ 2 + V ext + V H + V xc |ψ kj .

B. Use of Linear Muffin-Tin Orbitals
The next problem is to solve the Dyson equation for the eigenvalues. The sophisticated basis sets developed to solve the one-electron Schrödinger equation can be directly used in this case. We utilize the LMTO method described extensively in the past literature 42,43,44 as it provides a minimal atom-centered local orbital basis set ideally suited for the electronic structure calculation. Within the LMTO basis, the full Green function is represented as a sum and, as we have argued in the previous section, the matrix G αβ (k, iω) needs to be considered as a variable in the spectral density functional. The stationarity yields the equation for the Green function where the overlap matrix O αβ (k) = χ k α |χ k β , the non-interacting Hamiltonian matrix h (0) αβ (k) = χ k α |−∇ 2 +V ext (r)|χ k β and the self-energy formally comes as a matrix element over the LMTOs. Again, it is worth to point out that the self-energy here depends on k via the orbitals even if the single-impurity case is considered. In calculations performed on a cluster, the self-energy will also pick its non-trivial k-dependence coming from the nearest sites. While formally valid, the present approach is not very efficient since the Green function G(r, r ′ , iω) has to be evaluated via (86). This is the k-integral which has poles in a complex frequency plane, and integrating singular functions need to be performed with care. In this respect, we adopt the eigenvalue representation (71). We expand the energy-dependent states in terms of the LMTO basis {χ k α } as follows: The unknown coefficients A kjω,R(L) α are now the quantities which have to be considered as variables in the spectral density functional. The stationarity yields the equations These are the non-hermitian eigenvalue problems solved by standard numerical methods. The orthogonality condition involving the overlap matrix is Note that the present algorithm just inverts the matrix (87) with help of the "right" and "left" eigenstates. The Green function (87) in the basis of its eigenvectors becomes This formula can be safely used to compute the Green function as the integral over the Brillouin zone, Eq. (86), because the energy denominator can be integrated analytically using the tetrahedron method 61 . Our next topic here is the evaluation of the bath Green function G 0 (r, r ′ , iω). It can be found from the integral equation where r and r ′ run over Ω loc In order to solve this equation, it is useful to represent r = ρ + R, r ′ = ρ ′ + R ′ , and redenote G 0 (r, Considering one atom per unit cell let us see how this can be solved using single-kappa LMTO-ASA method. The generalization to multiatomic systems with multiple kappa basis sets as well as inclusion of full potential terms in the calculation can be done along the same lines. The form of the LMTO basis function inside the sphere centered at R is where S k Lα are the structure constants of the LMTO method and where Φ H,J L (ρ) are linear combinations of the solutions of the radial Schrödinger equation taken at spherical part of the potential and matched to the spherical Hankel (H) and Bessel (J) functions at the muffin-tin sphere boundary as well as their energy derivatives taken at some set of energies E νl at the center of interest. The local Green function can be represented in this basis set as follows where the matrices G HH loc,αRβR ′ (iω), G HJ loc,αRL ′ R ′ (iω), G JH loc,LRβR ′ (iω), G JJ loc,LRL ′ R ′ (iω) (indexes R and R ′ are restricted to a cluster) are given by the following Brillouin zone integrals Here are the convolutions of the eigenvectors with the LMTO structure constants. We now utilize a similar representation for the bath Green function where the matrices G HH 0,αRβR ′ (iω), can be found from the following Dyson equation (where the matrices sizes have been enlarged by a factor of 2) with the self-energy matrices are defined as follows The solution of the impurity model with G (µν) 0,LRL ′ R ′ (iω) delivers new matrix elements (100). The k-dependent selfenergy (88) to be used in constructing the new Green function in (87) is found first by restoring the k-dependence from the cluster M int,LRL ′ R ′ (iω)e ik(R−R ′ ) and second, restoring the k-dependence of the LMTO basis as follows (101) In practical calculations performed with the LDA+DMFT method for Pu, only the subset of orbitals {χ a } is treated as correlated (f electrons of Pu) and a single-impurity case is considered. It is useful to separate the Hartree and LDA exchange-correlations terms. Instead of dealing with the non-interacting Hamiltonian in (87), we can rearrange the contributions to arrive where M αβ (k, iω) = δ αa δ βbMab (k, iω) and V DC αβ (k) = δ αa δ βb V DC ab (k) represent the DMFT correction and double counting term described by (69). These matrices are non-zero within the correlated subset. To accelerate the calculation of the impurity model, we can parametrize the self-energy matrix asM . With such parametrization, the local Green function which enters the Dyson equation should be defined as follows G loc,ab (iω) = kj µν A kjω,R . This represents the generalization of a partial-density-of-state formula of the LMTO method. The bath Green function can be found from the equation: ab (iω) and can be passed to the impurity solver. The latter delivers a new self-energỹ M and used to reconstruct new k-dependent self-energy after (101). Such procedure preserves all k-dependent information coming from the orbitals.

IV. APPLICATIONS TO PLUTONIUM
This section describes the application of the theory to Plutonium. Pu is known to be an anomalous metal 62 . It has six crystallographic structures. Starting from the low temperature α phase (0 to 100 C) with 16 atoms per unit cell it shows a series of phase transitions and ends up in relatively simple fcc δ (300-450 C) and bcc ε phases (500-650 C) just before it melts. The temperature dependence of atomic volume in Pu is anomalous 63 . It shows an enormous volume expansion between α and δ phases which is about 25%. Within the δ phase, the metal shows negative thermal expansion. Transition between δ and higher-temperature ε phase occurs with a 5% volume collapse. Also, Pu shows anomalous resistivity behavior 64 characteristic for the heavy fermion systems, but neither of its phases is magnetic. The susceptibility is small and relatively temperature independent. The photoemission 65 shows a strong narrow Kondo-like peak at the Fermi level consistent with large values of the linear specific heat coefficient.
Density functional based LDA and GGA calculations describe the properties of Pu incorrectly. They predict magnetic ordering 66 . They underestimate 67 equilibrium volume of the δ and ε phase by as much as 30% , the largest discrepancy known in LDA, which usually predicts the volume of solids within a few % accuracy even for such correlated systems as high temperature superconductors. Despite this, the volume of the α phase is predicted correctly by LDA 67 . Since the transport and thermodynamic properties of α and δ Pu are very similar, the nature of the α phase and the reason why LDA predicts accurately its structure and volume is by itself is another puzzle.
To address these questions several approaches have been developed. The LDA+U method was applied to δ Pu 59,68 . It is able to produce the correct volume of the δ phase for values of the parameter U˜4 eV consistent with atomic spectral data and constrained density functional calculations. Similar calculation has been performed by a so-called orbitally ordered density functional method 69 . However, both methods predict Pu to be magnetic, which is not observed experimentally. The LDA+U method is unable to predict the correct excitation spectrum. Also, to recover the α phase within LDA+U the parameter U has to be set to zero which is inconsistent with its transport properties and with microscopic calculations of this parameter. Another approach proposed 70 in the past is the constrained LDA approach in which some of the 5f electrons, are treated as core, while the remaining are allowed to participate in band formation. Results of the self-interaction-corrected LDA calculations have been reported 71 , as well as qualitative discussion of the bonding nature across the actinides series has been given 72 .
Thus, the problem of Pu is challenging because its f-electrons are close to the Mott transition 73 . It provides us a crucial test for our quantitative theory of strong correlations. A short version of this work has appeared already 32 . Our implementation is based on the self-consistent LDA+DMFT method and uses the LMTO method in its tightbinding representation 43 . Spin-orbit coupling effects are important for actinide compounds and have been included in the calculation for Pu. The "full potential" terms have been neglected in the calculation through the use of the atomic sphere approximation with a one-kappa LMTO basis set. The necessary k-space integrals for evaluating Green functions and charge densities have been carried out using the tetrahedron method using (8,8,8) grid in the Brillouin zone.
We study in detail the total energy as a function of the parameter U and give our predictions for the volumes in α, δ, and ε phases. A comparative analysis of the one-electron spectra in both α and δ phases is also given and the comparison with the photoemission experiment 65 is made. Since the dynamical mean field theory requires the solution of the Anderson impurity model for the multi-orbital f-shell of Pu, we have developed a method which, inspired by the success of the iterative perturbation theory 10 , interpolates the self-energy between small and large frequencies. At low frequencies, the exact value of the self-energy and its slope is extracted from the Friedel sum rule and from a slave-boson treatment 37,38,39 . This approach is accurate as it has been shown recently to give the exact critical value of U in the large degeneracy limit 74 At high frequencies the self-energy behavior can be computed based on high-frequency moments expansions 14,40 . The result of interpolation can be encoded into a simple form for the self-energy like the following continuous fraction expansion where the unknown coefficients are determined to satisfy known conditions in the low-and high-frequency limits. This kind of self-energy fits Quantum Monte Carlo (QMC) data in large region of parameters, such as U and doping, and where this comparison is at all possible (small degeneracy and high temperature). Thus, our approach interpolates between four major limits: small and large iω's valid for any U as well as small and large U ′ s (band vs. atom) valid for any iω. The analytical continuation to the real frequency axis is not a problem with the present method and avoids the use of the Pade 75 and maximum entropy 4 based techniques.
Complete details of the method can be found in Ref. 41. Here we only mention technicalities connected to the f-electrons of Pu where we deal with the impurity Green functions which are the matrices 14x14. However, for the relativistic f-level in cubic symmetry, the matrices can be reduced to 5x5 with 4 non-zero off-diagonal elements. The solution of such impurity problem is still a formidable numerical problem. We therefore make some simplifications. First, the off-diagonal elements are in general small and will be neglected. We are left with the 5f 5/2 state split into 2 levels which are 2-fold (Γ 7 ), and 4-fold (Γ 8 ) degenerate, and with the 5f 7/2 state split into 3 levels which are 2-fold (Γ 6 ), 2-fold (Γ 7 ) and 4-fold (Γ 8 ) degenerate. Second, since in Pu the intermultiplet spin-orbit splitting is much larger than the intramultiplet crystal field splitting (> 5:1), we reduce the problem of solving AIM for the levels separately by treating the 5f 5/2 Γ 7 and Γ 8 levels as one 6-fold degenerate level, and the 5f 7/2 Γ 6 , Γ 7 and Γ 8 levels as another 8-fold degenerate level.

A. Calculation of Volume
We perform our calculations for δ and ǫ phases of Pu having simple fcc and bcc structures respectively. Since our method does not yet allow us to treat complicated lattices, only a simplified study of the α phase which formally has 16 atoms per unit-cell will be reported. The total energy as a function of volume is evaluated self-consistently using formula (80). It contains all the electrons including the core electrons. The local density approximation includes generalized gradient corrections after Ref. 76. Since the LDA+DMFT approximation (55) is used, we subtract from the LDA the average interaction energy of the f electrons in the form (67) of the double counting term and then add improved estimates of these quantities using the self-consistent solution of the impurity model.
To illustrate the importance of correlations, we discuss the results for various strengths of the on-site Coulomb interaction U . Fig. 5 reports our theoretical predictions. First, the total energy as a function of volume of the fcc lattice is computed. The temperature is fixed at 600K, i.e. in the vicinity of the region where the δ phase is stable. U =0 GGA curve indicates a minimum at V /V 0 =0.7. This volume is in fact close to the volume of the α phase. Certainly, we expect that correlations should be less important for the compressed lattice in general, but there is no sign whatever of the δ phase in the U =0 calculation. The total energy curve is dramatically different for U >0. The details depend sensitively on the actual value of U . The behavior at U =4 eV shows the possibility of a double minimum; it is actually realized for a slightly smaller value of U . We find that for U = 3.8 eV, the minimum occurs near V /V δ =0.80 which corresponds to the volume of the α phase if we allow for monoclinic distortions and a volume-dependent U . When U increases by 0.2 eV the minimum occurs at V /V δ =1.05 which corresponds to the volume of the δ phase, in close agreement with experiment. Since the energies are so similar, we may expect that as temperature decreases, the lattice undergoes a phase transition from the δ phase to the α phase with the remarkable decrease of the volume by 25%. We repeated our calculations for the bcc structure using the temperature T = 900K where the ǫ phase is stable. Fig. 5 shows these results for U = 4 eV with a location of the minimum at around V /V δ =1.03. While the theory has a residual inaccuracy in determining the δ and ǫ phase volumes by a few percent, a hint of volume decrease with the δ → ǫ transition is clearly reproduced. Thus, our first-principles calculations reproduce the main features of the experimental phase diagram of Pu.
Note that the values of U ∼4 eV which are needed in our simulation to describe the α → δ transition, are in good agreement with the values of on-site Coulomb repulsion between f-electrons estimated by atomic spectral data 77 , constrained density functional studies 78 , and our previous LDA+U studies 59 .
The double-well behavior in the total energy curve is unprecedented in LDA or GGA based calculations but it is a natural consequence of the proximity to a Mott transition. Indeed, recent studies of model Hamiltonian systems 10,79 in the vicinity of the Mott transition show that two DMFT solutions which differ in their spectral distributions can coexist. It is very natural that allowing the density to relax in these conditions can give rise to the double minima as seen in Fig. 5.

B. Calculation of Spectra
We now report our calculated spectral density of states for the fcc structure using the volume V /V δ =0.8 and V /V δ =1.05 corresponding to our theoretical studies of α and δ phases. To compare the results of the dynamical mean-field calculations with the LDA method as well as with the experiment, we discuss the results presented in Fig.  6. Fig. 6(a) shows the density of states calculated using LDA+DMFT method in the vicinity of the Fermi level. Solid black line corresponds to the δ phase and solid grey line corresponds to the α phase. We predict the appearance of a strong quasiparticle peak near the Fermi level which exists in the both phases. Also, the lower and upper Hubbard bands can be clearly distinguished in this plot. The width of the quasiparticle peak in the α phase is found to be larger by 30 per cent compared to the width in the δ phase. This indicates that the low-temperature phase is more metallic, i.e. it has larger spectral weight in the quasiparticle peak and smaller weight in the Hubbard bands. Recent advances have allowed the experimental determination of these spectra, and our calculations are consistent with these measurements 65 . Fig. 6(b) shows the measured photoemission spectrum for δ (black line) and α (gray line) Pu. We can clearly see a strong quasiparticle peak. Also a smaller peak located at 0.8 eV for the δ-phase can be found. We interpret it as the lower Hubbard band.
The result of the local density approximation is shown on Fig. 6(a) by dashed line. The LDA produces two peaks near the Fermi level corresponding to 5f 5/2 and 5f 7/2 states separated by the spin-orbit coupling. The Fermi level falls into the dip between these states and cannot reproduce the features seen in photoemission. We should also mention that LDA+U fails completely 59,68 to reproduce the intensity of the f -states near the Fermi level as it pushes the f -band 2-3 eV below the Fermi energy. This is the picture expected from the static Hartree-Fock theory such as the LDA+U. Only full inclusion of the dynamic effects within the DMFT allows to account for both the quasiparticle resonance and the Hubbard satellites which explains all features of the photoemission spectrum in δ Pu.
The calculated by LDA+DMFT densities of states at E F equal to 7 st./[eV*cell] are consistent with the measured values of the linear specific heat coefficient. We still find a residual discrepancy by about factor of 2 due to either inaccuracies of the present calculation or due to the electron-phonon interactions. However, these values represent an improvement as compared to the LDA calculations which appear to be 5 times smaller. Similar inaccuracy has been seen in the LDA+U calculation 59 .
A simple physical explanation drawn from these studies suggests that in the δ phase the f electrons are slightly on the localized side of the interaction-driven localization-delocalization transition with a sharp and narrow Kondo-like peak and well-defined upper and lower Hubbard bands. It therefore has the largest volume as has been found by previous LDA+U calculations 59,68 which take into account Hubbard bands only. The low-temperature α phase is more metallic, i.e. it has larger spectral weight in the quasiparticle peak and smaller weight in the Hubbard bands. It will therefore have a much smaller volume that is eventually reproduced by LDA/GGA calculations which neglect both Coulomb renormalizations of quasiparticles and atomic multiplet structure. The delicate balance of the energies of the two minima may be the key to understanding the anomalous properties of Pu such as the great sensitivity to small amounts of impurities (which intuitively would raise the energy of the less symmetric monoclinic structure, thus stabilizing the δ phase to lower temperature) and the negative thermal expansion. Notice however, that the α phase is not a weakly correlated phase: it is just slightly displaced towards the delocalized side of the localizationdelocalization transition, relative to the δ phase. This is a radical new viewpoint in the theoretical literature on Pu, which has traditionally regarded the α phase as well understood within LDA. However, the correlation viewpoint is consistent with a series of anomalous transport properties in the α phase reminiscent of heavy electron systems. For example, the resistivity of α-Pu around room temperature is anomalously large, temperature independent and above the Mott limit 64 (the maximum resistivity allowed to the conventional metal). Strong correlation anomalies are also evident in the thermoelectric power 80 .

V. CONCLUSION
In conclusion, this work describes a first principles method for calculating the electronic structure of materials where many-body correlation effects between the electrons are not small and cannot be neglected. It allows simultaneous evaluation of the total free energy and the local electronic spectral density. The approach is based on the effective action functional formulation of the free energy and is viewed as spectral density functional theory. An approximate form of the functional exploits a local dynamical mean-field theory of strongly correlated systems accurate in the situations when the self-energy is short ranged in a certain portion of space. The localization is defined with reference to some basis in Hilbert space. It does not necessarily imply localization in real space and is treated using a general basis set following the ideology of the cellular dynamical mean field theory. Further approximations of the theory, such as LDA+DMFT and local GW are discussed. Implementation of the method is described in terms of the energy-dependent one-particle states expanded via the linear muffin-tin orbitals. Application of the method in its LDA+DMFT form is given to study the anomalous volume expansion in metallic Plutonium. We obtain equilibrium volume of the δ phase in good agreement with experiment with no magnetic order imposed in the calculation. The calculated one-electron densities of states are consistent with the results of the photoemission. Our most recent studies 34 of the lattice dynamical properties of Pu address the problem of the δ → ε transitions and show good agreement with experiment 81 . Alternative developments of the LDA+DMFT approach by several groups around the world discuss other applications of the dynamical mean field theory in electronic structure calculations. The results obtained are promising. Volume collapse transitions, materials near the Mott transition, systems with itinerant and local moments, as well as many other exciting problems are beginning to be explored using these methods.
The authors would like to thank for many enlightening discussions the participants of the research school on Realistic Theories of Correlated Electron Materials, Kavli Institute for Theoretical Physics, fall 2002, where the part of this work was carried out. Many helpful discussions with the participants of weekly condensed matter seminar at Rutgers University are also acknowledged. Conversations with A. Georges, A Lichtenistein, N. E. Zein related to various aspects of dynamical mean-field and GW theories are much appreciated. The authors are indebted to E. Abrahams, A. Arko, J. Joyce, H. Ledbetter, A. Migliori, and J. Thompson for discussing the issues related to the work on Pu.