format: text

text:
# A Rigorous Extension of the Kohn-Sham Equation for Strongly Correlated Electron Systems

## 1 Introduction

The density functional theory (DFT) is a general theory to describe non-uniform interacting electron gas at zero temperature$`^{\text{?}\text{}\text{?}\text{}\text{?}\text{}\text{?}\text{}\text{?}\text{)}}`$ or in equilibrium at finite temperatures.$`^{\text{?}\text{, ?, ?)}}`$ This general formalism provides variational principle to search the ground state or the density matrix of the many-body system through minimization of a functional of the single-particle density. It would be allowed to say that the modern theory of the electronic band-structure calculations would not be established, if DFT were not known. DFT is a rigorous theory stating that the ground state energy is obtained as a minimum of an energy density functional.$`^{\text{?}\text{}\text{?}\text{}\text{?}\text{)}}`$

To perform realistic calculation, we have to utilize a few steps of reformulation with assumptions and approximations starting from the rigorous description. The first strategy is the Kohn-Sham (KS) scheme which provides us with an effective one-body problem.$`^{\text{?, }\text{?}\text{)}}`$ Namely, one can obtain an effective equation called the Kohn-Sham equation (KSE), which looks like the Schrödinger equation for a non-interacting system. Utilizing the so-called local-density approximation (LDA)$`^{\text{?, }\text{?}\text{)}}`$ or the generalized gradient approximation (GGA),$`^{\text{?}\text{}\text{?}\text{)}}`$ KSE becomes a tractable self-consistent equation with a Hartree term and a so-called exchange-correlation potential.

The DFT-LDA scheme has been regarded as a major strategy of the band-structure calculation. A lot of studies of the electronic states have been done to reveal that the KS scheme supplies a well-controlled theoretical description of various materials including semiconductors, ionic crystals and metals,$`^{\text{?)}}`$ as well as the large scale systems by virtue of the iterative minimization technique$`^{\text{?}\text{}\text{?}\text{)}}`$ and description by the localized orbitals.$`^{\text{?}\text{}\text{?}\text{}\text{?}\text{)}}`$

However, it has been recognized that the band calculation using DFT-LDA does not reproduce several kinds of correlation effects,$`^{\text{?, }\text{?}\text{)}}`$ e.g. the Mott metal-insulator transition and the Kondo effect. Usage of LDA (or GGA) has thought to be a serious reason for this difficulty. But, the KS scheme itself might not be suitable to describe these correlation effects. This is because the density is given by an auxiliary non-interacting system. For non-interacting Fermions, the ground state is given by a single Slater determinant. This wave function cannot describe the correlation effects. Especially, for correlation effects originating from short range repulsion, the difficulty can lead to qualitative errors.

Since DFT tells us only the single-particle density of the ground state, one might say that it is by no means hopeful to expect DFT as a general theory to reproduce every properties of materials. However, the present KS scheme has been practically proven to be a useful skill for materials scientists, even when they study single-particle excitations, i.e. the band structure. Thus, we intended to find an extension of the KS approach keeping its convenience but enlarging a range of applications.

We expected that a natural extension of the KS scheme is to utilize a density given by a many-body state of a simplified auxiliary problem. If this scheme becomes possible, not only single-particle excitations but also collective excitations might be well-described with reasonable accuracy. At present, many theoretical techniques have been developed to solve models with short-range interactions. The Hubbard model$`^{\text{?}\text{}\text{?}\text{}\text{?}\text{)}}`$ is a typical example, which is known to show various kinds of phase transitions and critical phenomena. Thus, if we can improve our technical scheme to unify DFT and knowledge of the strong correlation effects, our theoretical understanding of materials would be wider and richer.

We can say that the so-called LDA+U (and LDA++) approach has been developed along this line.$`^{\text{?}\text{}\text{?}\text{}\text{?}\text{, ?)}}`$ In the LDA+U, the KS orbitals are used as a basis for an additional many-body problem. Indeed, idea of LDA++ is defined as a unified strategy in which the dynamical mean-field approximation (DMF)$`^{\text{?}\text{)}}`$ is used to solve a many-body Hamiltonian, which is essentially an extended Hubbard model given by the KS orbitals.$`^{\text{?)}}`$ As well-known, DMF is a method to deal with short-range dynamical correlation correctly. However, the LDA+U energy functional is an artificial functional obtained by adding extra U-terms and subtracting a counter term, which is valid only in an atomic limit. Thus, it is desired to formulate a strategy to derive a proper energy functional by a mathematically reliable approach. As we will show, it is possible to determine the energy functional using a variational method which is same as that used in the Kohn-Sham theory.

In this paper, we present a version of our formalism to obtain an extension of the Kohn-Sham equation. In this version, we consider a non-relativistic many-electron system. Thus, we will start from the energy functional given by Levy.$`^{\text{?)}}`$ We will introduce model functionals, which are composed of a kinetic energy functional, an interaction energy functional and a new exchange-correlation functional. A single particle density is given by a many-body state which is a solution of an auxiliary many-body problem. The energy functionals proposed in the present article are essentially extensions of the KS functionals.

The original Kohn-Sham scheme has derived based on a next assumption. Namely, the true single-particle density was assumed to be reproduced by a single Slater determinant and the Slater determinant was assumed to be a solution of an effective single-particle equation. This assumption is a so-called v-representability in KSE. To avoid usage of the v-representability, Hadjisavvas and Theophilou (HT) have introduced functionals on the space of Slater determinants instead of density functionals. Their strategy succeeded in bypassing the v-representability problem. However, their argument is a bit complicated for intuitive argument. Thus we will start our derivation using a variational method similar to the Kohn-Sham scheme from §2. A mathematically rigorous definition of our new functionals in accordance with the HT argument will be also presented in §6.

## 2 Definition of the Problem

We consider a non-uniform electron gas in an external scalar potential $`v_{\mathrm{ext}}(𝐫)`$. Here the position vector $`𝐫`$ is defined in a three-dimensional space. The problem is given by the next non-relativistic Hamiltonian,

$$\widehat{H}=\widehat{T}+\widehat{V}+d^3rv_{\mathrm{ext}}(𝐫)\widehat{n}(𝐫).$$
(2.1)
Using electron field operators with spin $`\sigma `$, which are denoted as $`\psi _\sigma (𝐫)`$ and $`\psi _\sigma ^{}(𝐫)`$, the kinetic energy operator $`\widehat{T}`$ and the 2-body interaction $`\widehat{V}`$ are given as follows.

$$\widehat{T}=\frac{\mathrm{}^2}{2m}d^3r\underset{\sigma }{}\underset{𝐫^{}𝐫}{lim}\psi _\sigma ^{}(𝐫^{})\mathrm{\Delta }_𝐫\psi _\sigma (𝐫).$$
(2.2)
$$\widehat{V}=\frac{1}{2}d^3rd^3r^{}\frac{e^2}{|𝐫𝐫^{}|}\underset{\sigma ,\sigma ^{}}{}\psi _\sigma ^{}(𝐫)\psi _\sigma ^{}^{}(𝐫^{})\psi _\sigma ^{}(𝐫^{})\psi _\sigma (𝐫).$$
(2.3)
The operator, $`\widehat{n}(𝐫)_\sigma \psi _\sigma ^{}(𝐫)\psi _\sigma (𝐫)`$, gives the electron density, $`n(𝐫)`$, at $`𝐫`$.

The many-body problem given by $`\widehat{H}`$ is difficult to solve in general. But we know that the system is stable, i.e. the energy is bounded below and that the ground state exists for any $`v_{\mathrm{ext}}(𝐫)`$. Thus we introduce $`|\mathrm{\Phi }_{\mathrm{GS}}`$ representing the ground state of $`\widehat{H}`$. The single particle density of $`|\mathrm{\Phi }_{\mathrm{GS}}`$ is given by,

$$n_{\mathrm{GS}}(𝐫)=\mathrm{\Phi }_{\mathrm{GS}}|\widehat{n}(𝐫)|\mathrm{\Phi }_{\mathrm{GS}}.$$
(2.4)
The rigorous version of DFT has been established by introduction of the next energy functional by Levy.$`^{\text{?)}}`$

$$F[n]=\underset{\mathrm{\Phi }n}{\mathrm{min}}\mathrm{\Phi }|\widehat{T}+\widehat{V}\mathrm{\Phi }.$$
(2.5)
Here, $`\mathrm{min}_{\mathrm{\Phi }n}`$ represents the constraint minimization, where $`\mathrm{\Phi }`$ is optimized under a constraint that $`\mathrm{\Phi }|\widehat{n}(𝐫)|\mathrm{\Phi }=n(𝐫)`$. We know that the minimum for the Coulomb problem under discussion exists and that $`F[n]`$ is a well-defined functional of $`n(𝐫)`$.$`^{\text{?, ?, ?)}}`$ It is recognized that the strategy of the KS approach should be based on the functional, eq. (2.5).$`^{\text{?, ?, ?)}}`$

The central assumption to derive our extension of the KS scheme is the following. We assume that there exists an auxiliary many-body system and that the ground state of the model system, $`|\mathrm{\Psi }_{\mathrm{GS}}`$, satisfies, $`\mathrm{\Psi }_{\mathrm{GS}}|\widehat{n}(𝐫)|\mathrm{\Psi }_{\mathrm{GS}}=n_{\mathrm{GS}}(𝐫)`$. Our purpose is to find a systematic strategy to determine the model system. This model should be described by a set of auxiliary equations, which we introduce now.

The first set of equations is to define single-particle orbitals.

$$\left\{\frac{\mathrm{}^2}{2m}\mathrm{\Delta }_𝐫+v_{\mathrm{eff}}(𝐫)\right\}\varphi _{i,\sigma }(𝐫)=\epsilon _i\varphi _{i,\sigma }(𝐫)$$
(2.6)
Here $`v_{\mathrm{eff}}(𝐫)`$ is an effective potential, which will be determined by a variational method. Spin-dependence of $`v_{\mathrm{eff}}(𝐫)`$ does not appear in our non-relativistic version, since we consider only the density functional and do not deal with a spin-density functional in the present paper. However, as we discuss in the last section, $`|\mathrm{\Psi }_{\mathrm{GS}}`$ can be magnetic.

We assume existence of $`v_{\mathrm{eff}}(𝐫)`$ in the same manner as Kohn and Sham did. Here, we should note that the existence of $`v_{\mathrm{eff}}(𝐫)`$ assumed in the original KS scheme has not been proved nor disproved, although the v-representability assumed in the Hohenberg-Kohn theory was disproved.$`^{\text{?, ?, }\text{?}\text{)}}`$ In §6, we will utilize the HT technique to redefine our scheme avoiding the v-representability problem.

Since eq. (2.6) determines a complete set of single-particle wave functions, we can introduce Fermion operators, $`c_{i,\sigma }^{}`$ and $`c_{i,\sigma }`$. They create (or annihilate) the $`i`$-th state, $`\varphi _{i,\sigma }(𝐫)`$, with a spin $`\sigma `$ and satisfy canonical commutation relations for the Fermion operators, $`\{c_{i,\sigma },c_{j,\sigma ^{}}^{}\}=\delta _{i,j}\delta _{\sigma ,\sigma ^{}}`$. Note that we may use another basis set comprised of localized orbitals e.g. the Wannier basis, to consider an auxiliary many-body problem. However, the kinetic energy term becomes a tight-binding Hamiltonian with hopping terms between every pair of orbitals. Thus, we utilize the eigen states of eq. (2.6) to show the formal strategy. We will soon explain the reason why the formulation is easily generalized using localized electron picture.

Using $`c_{i,\sigma }^{}`$ and $`c_{i,\sigma }`$, the effective Hamiltonian of the auxiliary many-body problem is written as follows.

$`\widehat{H}_{\mathrm{eff}}`$ $`=`$ $`{\displaystyle \underset{i,\sigma }{}}\epsilon _ic_{i,\sigma }^{}c_{i,\sigma }+{\displaystyle \underset{i,j,\sigma }{}}V_{\mathrm{eff}}^{(1)}(i,j,\sigma )c_{i,\sigma }^{}c_{j,\sigma }`$ (2.7)
$`+`$ $`{\displaystyle \underset{i,j,k,l,\sigma ,\sigma ^{}}{\overset{(2)}{}}}V_{\mathrm{eff}}^{(2)}(i,j,k,l;\sigma ,\sigma ^{})c_{i,\sigma }^{}c_{j,\sigma ^{}}^{}c_{k,\sigma ^{}}c_{l,\sigma }`$
$`+`$ $`\widehat{H}_{\mathrm{rxc}}.`$
Thus the second equation of the auxiliary problem is,

$$\widehat{H}_{\mathrm{eff}}|\mathrm{\Psi }=E|\mathrm{\Psi }.$$
(2.8)
Again, effective interaction parameters, $`V_{\mathrm{eff}}^{(1)}(i,j,\sigma )`$ and $`V_{\mathrm{eff}}^{(2)}(i,j,k,l;\sigma ,\sigma ^{})`$, should be determined by a variational method. The following sections is devoted to explanation of a method to determine these effective interactions and the effective potentials in a self-consistent manner. The last term, $`\widehat{H}_{\mathrm{rxc}}`$, in eq. (2.7) represents residual exchange-correlation interactions not described by two-particle interactions, $`V_{\mathrm{eff}}^{(2)}(i,j,k,l;\sigma ,\sigma ^{})`$. Definition of $`\widehat{H}_{\mathrm{rxc}}`$ will be given also using the variational method.

An important point of our method is that a summation denoted by $`^{(2)}`$ is arbitrary as far as $`\widehat{H}_{\mathrm{eff}}`$ is an Hermitian operator. Namely, we have degrees of freedom to optimize the auxiliary effective problem. In other words, we may consider Hubbard-type localized interactions at any time by considering a proper unitary transformation from extended states, $`\varphi _{i,\sigma }(𝐫)`$, to localized orbitals.

The many-body state, $`|\mathrm{\Psi }`$, should be given as a summation of Slater determinants. We use a term, “a multi Slater determinant” as a summation of Slater determinants below. In the present notation, it reads,

$$|\mathrm{\Psi }=\underset{l=\{i_1,\sigma _1,\mathrm{},i_N,\sigma _N\}}{}f_lc_{i_1,\sigma _1}^{}c_{i_2,\sigma _2}^{}\mathrm{}c_{i_N,\sigma _N}^{}|0,$$
(2.9)
for an $`N`$-particle system, where $`l`$ represents a multi-index, $`\{i_1,\sigma _1,\mathrm{},i_N,\sigma _N\}`$, and coefficients $`f_l`$ should be determined by solving eq. (2.8). Namely, we assume that the ground state of the second auxiliary problem, $`|\mathrm{\Psi }_{\mathrm{GS}}`$, is obtainable, when we derive the formal theory.

We introduce the reduced density matrices.

$`\rho (i,j,\sigma )`$ $`=`$ $`\mathrm{\Psi }|c_{i,\sigma }^{}c_{j,\sigma }|\mathrm{\Psi },`$ (2.10)
$`\rho ^{(2)}(i,j,k,l;\sigma ,\sigma ^{})`$ $`=`$ $`\mathrm{\Psi }|c_{i,\sigma }^{}c_{j,\sigma ^{}}^{}c_{k,\sigma ^{}}c_{l,\sigma }|\mathrm{\Psi },`$ (2.11)
$`\overline{\rho }^{(2)}(i,j,k,l;\sigma ,\sigma ^{})`$ $`=`$ $`\rho ^{(2)}(i,j,k,l;\sigma ,\sigma ^{})`$ (2.12)
$``$ $`\rho (i,l,\sigma )\rho (j,k,\sigma ^{}).`$
The single particle density of $`|\mathrm{\Psi }`$ is given as,

$`n(𝐫)`$ $`=`$ $`\mathrm{\Psi }|\widehat{n}(𝐫)|\mathrm{\Psi }`$ (2.13)
$`=`$ $`\mathrm{\Psi }|{\displaystyle \underset{\sigma }{}}\psi _\sigma ^{}(𝐫)\psi _\sigma (𝐫)|\mathrm{\Psi }`$
$`=`$ $`\mathrm{\Psi }|{\displaystyle \underset{i,j,\sigma }{}}\varphi _{i,\sigma }^{}(𝐫)\varphi _{j,\sigma }(𝐫)c_{i,\sigma }^{}c_{j,\sigma }|\mathrm{\Psi }`$
$`=`$ $`{\displaystyle \underset{i,j,\sigma }{}}\varphi _{i,\sigma }^{}(𝐫)\varphi _{j,\sigma }(𝐫)\rho (i,j,\sigma ).`$
Following the central assumption, $`n(𝐫)=n_{\mathrm{GS}}(𝐫)`$, if the auxiliary problem is properly chosen and if $`|\mathrm{\Psi }=|\mathrm{\Psi }_{\mathrm{GS}}`$.

The total energy of the auxiliary model is given as,

$`E`$ $`=`$ $`\mathrm{\Psi }|\widehat{H}_{\mathrm{eff}}|\mathrm{\Psi }`$ (2.14)
$`=`$ $`{\displaystyle \underset{i,\sigma }{}}\epsilon _i\rho (i,i,\sigma )+{\displaystyle \underset{i,j,\sigma }{}}V_{\mathrm{eff}}^{(1)}(i,j,\sigma )\rho (i,j,\sigma )`$
$`+`$ $`{\displaystyle \underset{i,j,k,l,\sigma ,\sigma ^{}}{\overset{(2)}{}}}V_{\mathrm{eff}}^{(2)}(i,j,k,l;\sigma ,\sigma ^{})\rho ^{(2)}(i,j,k,l;\sigma ,\sigma ^{})`$
$`+`$ $`\mathrm{\Psi }|\widehat{H}_{\mathrm{rxc}}|\mathrm{\Psi }`$
## 3 Model Energy Functional

To find an extension of KSE, the second important step is determination of energy functionals for the auxiliary model. We define a kinetic energy functional $`T_0[n]`$ and an interaction energy functional $`V_{\mathrm{ee}}[n]`$ as follows.

$`T_0[n]`$ (3.1)
$`=`$ $`\mathrm{\Psi }|\left[{\displaystyle \underset{i,j,\sigma }{}}{\displaystyle d^3r\varphi _{i,\sigma }^{}(𝐫)\left(\frac{\mathrm{}^2}{2m}\mathrm{\Delta }_𝐫\right)\varphi _{j,\sigma }(𝐫)c_{i,\sigma }^{}c_{j,\sigma }}\right]|\mathrm{\Psi }`$
$`=`$ $`{\displaystyle \underset{i,j,\sigma }{}}{\displaystyle d^3r\varphi _{i,\sigma }^{}(𝐫)\left(\frac{\mathrm{}^2}{2m}\mathrm{\Delta }_𝐫\right)\varphi _{j,\sigma }(𝐫)\rho (i,j,\sigma )}`$
$`=`$ $`{\displaystyle \underset{i,j,\sigma }{}}\epsilon _i\delta _{i,j}\rho (i,j,\sigma )`$
$``$ $`{\displaystyle \underset{i,j,\sigma }{}}{\displaystyle d^3r\varphi _{i,\sigma }^{}(𝐫)v_{\mathrm{eff}}(𝐫)\varphi _{j,\sigma }(𝐫)\rho (i,j,\sigma )}`$
$`=`$ $`{\displaystyle \underset{i,\sigma }{}}\epsilon _i\rho (i,i,\sigma ){\displaystyle \underset{i,j,\sigma }{}}{\displaystyle d^3rv_{\mathrm{eff}}(𝐫)n(𝐫)}.`$
$`V_{\mathrm{ee}}[n]`$ (3.2)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle d^3rd^3r^{}\frac{e^2}{|𝐫𝐫^{}|}n(𝐫)n(𝐫^{})}`$
$`+`$ $`{\displaystyle \underset{i,j,k,l,\sigma ,\sigma ^{}}{\overset{(2)}{}}}{\displaystyle \frac{1}{2}}{\displaystyle d^3rd^3r^{}\frac{e^2}{|𝐫𝐫^{}|}\varphi _{i,\sigma }^{}(𝐫)\varphi _{j,\sigma ^{}}^{}(𝐫^{})}`$
$`\times \varphi _{k,\sigma ^{}}(𝐫^{})\varphi _{l,\sigma }(𝐫)\overline{\rho }^{(2)}(i,j,k,l;\sigma ,\sigma ^{})`$
$`=`$ $`{\displaystyle \underset{i,j,k,l,\sigma ,\sigma ^{}}{}}{\displaystyle \frac{1}{2}}(\varphi _{i,\sigma }\varphi _{j,\sigma ^{}}|{\displaystyle \frac{e^2}{r}}|\varphi _{k,\sigma ^{}}\varphi _{l,\sigma })\rho (i,l,\sigma )\rho (j,k,\sigma )`$
$`+`$ $`{\displaystyle \underset{i,j,k,l,\sigma ,\sigma ^{}}{\overset{(2)}{}}}{\displaystyle \frac{1}{2}}(\varphi _{i,\sigma }\varphi _{j,\sigma ^{}}|{\displaystyle \frac{e^2}{r}}|\varphi _{k,\sigma ^{}}\varphi _{l,\sigma })\overline{\rho }^{(2)}(i,j,k,l;\sigma ,\sigma ^{})`$
$`=`$ $`{\displaystyle \underset{i,j,k,l,\sigma ,\sigma ^{}}{\overset{(2)}{}}}{\displaystyle \frac{1}{2}}(\varphi _{i,\sigma }\varphi _{j,\sigma ^{}}|{\displaystyle \frac{e^2}{r}}|\varphi _{k,\sigma ^{}}\varphi _{l,\sigma })\mathrm{\Psi }|c_{i,\sigma }^{}c_{j,\sigma ^{}}^{}c_{k,\sigma ^{}}c_{l,\sigma }|\mathrm{\Psi }`$
$`+`$ $`{\displaystyle \underset{i,j,k,l,\sigma ,\sigma ^{}}{\overset{(2)^{}}{}}}{\displaystyle \frac{1}{2}}(\varphi _{i,\sigma }\varphi _{j,\sigma ^{}}|{\displaystyle \frac{e^2}{r}}|\varphi _{k,\sigma ^{}}\varphi _{l,\sigma })\mathrm{\Psi }|c_{i,\sigma }^{}c_{l,\sigma }|\mathrm{\Psi }`$
$`\times \mathrm{\Psi }|c_{j,\sigma ^{}}^{}c_{k,\sigma ^{}}|\mathrm{\Psi }.`$
Here an interaction parameter is defined as,

$`(\varphi _{i,\sigma }\varphi _{j,\sigma ^{}}|{\displaystyle \frac{e^2}{r}}|\varphi _{k,\sigma ^{}}\varphi _{l,\sigma })`$
$`=`$ $`{\displaystyle d^3rd^3r^{}\frac{e^2}{|𝐫𝐫^{}|}\varphi _{i,\sigma }^{}(𝐫)\varphi _{j,\sigma ^{}}^{}(𝐫^{})\varphi _{k,\sigma ^{}}(𝐫^{})\varphi _{l,\sigma }(𝐫)},`$
and a summation $`^{(2)^{}}`$ represents $`^{(2)^{}}=^{(2)}`$.

We have to caution that in both eq. (3.1) and eq. (3.2) $`|\mathrm{\Psi }`$ is implicitly assumed to be $`|\mathrm{\Psi }_{\mathrm{GS}}`$. Namely, $`|\mathrm{\Psi }`$ is optimized with a condition, $`\mathrm{\Psi }|\widehat{n}(𝐫)|\mathrm{\Psi }=n(𝐫)`$. Then, the reason why $`T_0[n]`$ and $`V_{\mathrm{ee}}[n]`$ are functionals of the density is explained as follows. Since orbitals, $`\varphi _{i,\sigma }(𝐫)`$, are decided by a KS-like equation, eq. (2.6), these orbitals are a functional of $`n(𝐫)`$. Since the second equation, eq.(2.8), is determined by $`\varphi _{i,\sigma }(𝐫)`$, $`|\mathrm{\Psi }_{\mathrm{GS}}`$ and resulting density matrices are also a functional of $`n(𝐫)`$. Details of the simulation process are also discussed in §5.2. However, the above reasoning, which accords with the original KS argument, is not accurate enough. Namely, we have to specify a minimization procedure in detail. As clarified in §6, this process can be regarded as a search in the space of $`n(𝐫)`$. Thus, let us utilize notations $`T_0[n]`$ and $`V_{\mathrm{ee}}[n]`$.

Another energy functional which we call the residual exchange-correlation energy functional, $`E_{\mathrm{rxc}}[n]`$, is defined as the difference between the true energy functional and $`T_0[n]+V_{\mathrm{ee}}[n]`$ as,

$$F[n]=T_0[n]+V_{\mathrm{ee}}[n]+E_{\mathrm{rxc}}[n].$$
(3.3)
Thus the total energy functional for an electron system in an external potential $`v_{\mathrm{eff}}(𝐫)`$ is given as follows.

$$E[n]T_0[n]+V_{\mathrm{ee}}[n]+E_{\mathrm{rxc}}[n]+d^3rv_{\mathrm{eff}}(𝐫)n(𝐫).$$
(3.4)
Here we should comment the next point. In the last expression of eq. (3.2), selected interaction processes are explicitly kept as 2-body interactions in $`^{(2)}`$. But, for other processes, a corresponding Hartree term is used and corrections from remaining exchange-correlation effects are included in $`E_{\mathrm{rxc}}[n]`$.

## 4 Variational Method

To determine effective interactions and effective potentials, we utilize a variational method. Consider variation of $`E[n]`$ with respect to variation of $`n(𝐫)`$. The variation, $`\delta n(𝐫)`$, is separated into two parts.

$`\delta n(𝐫)`$ $`=`$ $`\delta n^{(o)}(𝐫)+\delta n^{(d)}(𝐫),`$ (4.1)
$`\delta n^{(o)}(𝐫)`$ $`=`$ $`{\displaystyle \underset{i,j,\sigma }{}}\left(\delta \varphi _{i,\sigma }^{}(𝐫)\varphi _{j,\sigma }(𝐫)+\varphi _{i,\sigma }^{}(𝐫)\delta \varphi _{j,\sigma }(𝐫)\right)`$ (4.2)
$`\times \rho (i,j,\sigma ),`$
$`\delta n^{(d)}(𝐫)`$ $`=`$ $`{\displaystyle \underset{i,j,\sigma }{}}\varphi _{i,\sigma }^{}(𝐫)\varphi _{j,\sigma }(𝐫)\delta \rho (i,j,\sigma ).`$ (4.3)
These equations mean that the density variation can be taken with respect to two types of independent directions in the functional space: one is by $`\delta \varphi _{i,\sigma }(𝐫)`$ and the other is by $`\delta \rho (i,j,\sigma )`$. To be more precise, we note that our model problem is solved by performing two successive steps alternately. The first one is to determine the single-particle orbitals given by eq. (2.6) fixing all density matrices and/or coefficients $`f_l`$ of the many-body state $`|\mathrm{\Psi }_{\mathrm{GS}}`$. After redefining $`\varphi _{i,\sigma }(𝐫)`$, the many-body problem given by eq. (2.8) is solved and density matrices are rebuilt. Thus, every variation of each physical quantity is taken with respect to $`\delta \varphi _{i,\sigma }(𝐫)`$ and $`\delta \rho (i,j,\sigma )`$. Our system of equations, i.e. eq. (2.6) and eq. (2.8), has to be solved self-consistently as a whole, so that $`\delta n^{(o)}(𝐫)`$ and $`\delta n^{(d)}(𝐫)`$ are mutually affected with each other. However, to derive determining equations we can consider the first order variation for each quantity. Therefor, we can classify variation into two types, i.e. those related to orbital variation and others related to variation of the density-matrices which are denoted by $`(o)`$ and $`(d)`$.

### 4.1 Variation of the Energy Functional with respect to the Single-Particle Orbitals

In this subsection, we consider $`\delta n^{(o)}(𝐫)`$ which is directly given by $`\delta \varphi _{i,\sigma }(𝐫)`$. Before presenting the result, we derive some relations. Let us consider variation of eq. (2.6), where $`\epsilon _i`$ and $`v_{\mathrm{eff}}(𝐫)`$ are also dependent variables of $`\delta n(𝐫)`$.

$`\left\{{\displaystyle \frac{\mathrm{}^2}{2m}}\mathrm{\Delta }+v_{\mathrm{eff}}(𝐫)\right\}\delta \varphi _{i,\sigma }(𝐫)+\delta v_{\mathrm{eff}}(𝐫)\varphi _{i,\sigma }(𝐫)`$
$`=\epsilon _i\delta \varphi _{i,\sigma }(𝐫)+\delta \epsilon _i\varphi _{i,\sigma }(𝐫)`$ (4.4)
Multiplying eq. (4.4) by $`\varphi _{i,\sigma }^{}(𝐫)`$ (or $`\varphi _{j,\sigma }^{}(𝐫)`$) from the left, integrating the result with respect to $`𝐫`$, we obtain the following useful identities.

$$\delta \epsilon _i=d^3r\delta v_{\mathrm{eff}}(𝐫)|\varphi _{i,\sigma }(𝐫)|^2,$$
(4.5)
$`(\epsilon _j\epsilon _i){\displaystyle d^3r\varphi _{j,\sigma }^{}(𝐫)\delta \varphi _{i,\sigma }(𝐫)}`$ (4.6)
$`=`$ $`{\displaystyle d^3r\delta v_{\mathrm{eff}}(𝐫)\varphi _{j,\sigma }^{}(𝐫)\varphi _{i,\sigma }(𝐫)}.`$
Variation of the total energy functional due to $`\delta n^{(o)}(𝐫)`$ is obtained as follows.

$`\delta E^{(o)}[n]`$ $`=`$ $`\delta T_0^{(o)}[n]+\delta V_{\mathrm{ee}}^{(o)}[n]+\delta E_{\mathrm{rxc}}^{(o)}[n]`$ (4.7)
$`+{\displaystyle d^3rv_{\mathrm{ext}}(𝐫)\delta n^{(o)}(𝐫)},`$
where $`\delta T_0^{(o)}[n]`$ and $`\delta V_{\mathrm{ee}}^{(o)}[n]`$ are given as,

$`\delta T_0^{(o)}[n]`$ (4.8)
$`=`$ $`{\displaystyle \underset{i,\sigma }{}}\delta \epsilon _i\rho (i,i,\sigma ){\displaystyle d^3r\delta v_{\mathrm{eff}}(𝐫)n(𝐫)}`$
$`{\displaystyle d^3rv_{\mathrm{eff}}(𝐫)\delta n^{(o)}(𝐫)}`$
$`=`$ $`{\displaystyle \underset{ij,\sigma }{}}{\displaystyle d^3r\delta v_{\mathrm{eff}}(𝐫)\varphi _{i,\sigma }^{}(𝐫)\varphi _{j,\sigma }(𝐫)\rho (i,j,\sigma )}`$
$`{\displaystyle d^3rv_{\mathrm{eff}}(𝐫)\delta n^{(o)}(𝐫)}`$
$`=`$ $`{\displaystyle \underset{ij,\sigma }{}}(\epsilon _i\epsilon _j){\displaystyle d^3r\varphi _{i,\sigma }^{}(𝐫)\delta \varphi _{j,\sigma }(𝐫)\rho (i,j,\sigma )}`$
$`{\displaystyle d^3rv_{\mathrm{eff}}(𝐫)\delta n^{(o)}(𝐫)},`$
$`\delta V_{\mathrm{eff}}^{(o)}[n]`$ (4.9)
$`=`$ $`{\displaystyle \underset{i,j,k,l,\sigma ,\sigma ^{}}{\overset{(2)}{}}}{\displaystyle \frac{1}{2}}\{(\delta \varphi _{i,\sigma }\varphi _{j,\sigma ^{}}|{\displaystyle \frac{e^2}{r}}|\varphi _{k,\sigma ^{}}\varphi _{l,\sigma })`$
$`+(\varphi _{i,\sigma }\delta \varphi _{j,\sigma ^{}}|{\displaystyle \frac{e^2}{r}}|\varphi _{k,\sigma ^{}}\varphi _{l,\sigma })+(\varphi _{i,\sigma }\varphi _{j,\sigma ^{}}|{\displaystyle \frac{e^2}{r}}|\delta \varphi _{k,\sigma ^{}}\varphi _{l,\sigma })`$
$`+(\varphi _{i,\sigma }\varphi _{j,\sigma ^{}}|{\displaystyle \frac{e^2}{r}}|\varphi _{k,\sigma ^{}}\delta \varphi _{l,\sigma })\}\overline{\rho }^{(2)}(i,j,k,l;\sigma ,\sigma ^{})`$
$`+{\displaystyle d^3rd^3r^{}\frac{e^2}{|𝐫𝐫^{}|}n(𝐫^{})\delta n^{(o)}(𝐫)}.`$
Thus we obtain the next expression for $`\delta E^{(o)}`$.

$`\delta E^{(o)}[n]`$ (4.10)
$`=`$ $`{\displaystyle \underset{ij,\sigma }{}}(\epsilon _i\epsilon _j){\displaystyle d^3r\varphi _{i,\sigma }^{}(𝐫)\delta \varphi _{j,\sigma }(𝐫)\rho (i,j,\sigma )}`$
$`+`$ $`{\displaystyle \underset{i,j,k,l,\sigma ,\sigma ^{}}{\overset{(2)}{}}}{\displaystyle \frac{1}{2}}\delta (\varphi _{i,\sigma }\varphi _{j,\sigma ^{}}|{\displaystyle \frac{e^2}{r}}|\varphi _{k,\sigma ^{}}\varphi _{l,\sigma })\overline{\rho }^{(2)}(i,j,k,l;\sigma ,\sigma ^{})`$
$`+`$ $`{\displaystyle }d^3r\{{\displaystyle }d^3r^{}{\displaystyle \frac{e^2}{|𝐫𝐫^{}|}}n(𝐫^{})+{\displaystyle \frac{\delta E_{\mathrm{rxc}}}{\delta n(𝐫)}}+v_{\mathrm{ext}}(𝐫)`$
$`v_{\mathrm{eff}}(𝐫)\}\delta n^{(o)}(𝐫).`$
### 4.2 Variation of the Energy Functional with respect to the Density Matrices

The variation $`\delta n^{(d)}(𝐫)`$ is given by variation of the density matrices, $`\delta \rho (i,j,\sigma )`$ and $`\delta \rho ^{(2)}(i,j,k,l;\sigma ,\sigma ^{})`$. Variation of energy functionals, $`\delta T_0^{(d)}[n]`$ and $`\delta V_{\mathrm{ee}}^{(d)}[n]`$, due to $`\delta n^{(d)}(𝐫)`$ are given as,

$`T_0^{(d)}[n]`$ (4.11)
$`=`$ $`{\displaystyle \underset{i,\sigma }{}}\epsilon _i\delta \rho (i,i,\sigma ){\displaystyle d^3rv_{\mathrm{eff}}(𝐫)\delta n^{(d)}(𝐫)},`$
$`V_{\mathrm{ee}}^{(d)}[n]`$ (4.12)
$`=`$ $`{\displaystyle \underset{i,j,k,l,\sigma ,\sigma ^{}}{\overset{(2)}{}}}{\displaystyle \frac{1}{2}}(\varphi _{i,\sigma }\varphi _{j,\sigma ^{}}|{\displaystyle \frac{e^2}{r}}|\varphi _{k,\sigma ^{}}\varphi _{l,\sigma })\delta \overline{\rho }^{(2)}(i,j,k,l;\sigma ,\sigma ^{})`$
$`+{\displaystyle d^3rd^3r^{}\frac{e^2}{|𝐫𝐫^{}|}n(𝐫^{})\delta n^{(d)}(𝐫)}.`$
We obtain the next expression for $`\delta E^{(d)}`$.

$`\delta E^{(d)}[n]`$ (4.13)
$`=`$ $`{\displaystyle \underset{i,\sigma }{}}\epsilon _i\delta \rho (i,i,\sigma )`$
$`+`$ $`{\displaystyle \underset{i,j,k,l,\sigma ,\sigma ^{}}{\overset{(2)}{}}}{\displaystyle \frac{1}{2}}(\varphi _{i,\sigma }\varphi _{j,\sigma ^{}}|{\displaystyle \frac{e^2}{r}}|\varphi _{k,\sigma ^{}}\varphi _{l,\sigma })\delta \overline{\rho }^{(2)}(i,j,k,l;\sigma ,\sigma ^{})`$
$`+`$ $`{\displaystyle }d^3r\{{\displaystyle }d^3r^{}{\displaystyle \frac{e^2}{|𝐫𝐫^{}|}}n(𝐫^{})+{\displaystyle \frac{\delta E_{\mathrm{rxc}}}{\delta n(𝐫)}}+v_{\mathrm{ext}}(𝐫)`$
$`v_{\mathrm{eff}}(𝐫)\}\delta n^{(d)}(𝐫).`$
Thus the total variation of the energy functional is given as follows.

$`\delta E[n]`$ (4.14)
$`=`$ $`{\displaystyle \underset{i,\sigma }{}}\epsilon _i\delta \rho (i,i,\sigma )`$
$`+`$ $`{\displaystyle \underset{ij,\sigma }{}}(\epsilon _i\epsilon _j){\displaystyle d^3r\varphi _{i,\sigma }^{}(𝐫)\delta \varphi _{j,\sigma }(𝐫)\rho (i,j,\sigma )}`$
$`+`$ $`{\displaystyle \underset{i,j,k,l,\sigma ,\sigma ^{}}{\overset{(2)}{}}}{\displaystyle \frac{1}{2}}\{\delta (\varphi _{i,\sigma }\varphi _{j,\sigma ^{}}|{\displaystyle \frac{e^2}{r}}|\varphi _{k,\sigma ^{}}\varphi _{l,\sigma })\overline{\rho }^{(2)}(i,j,k,l;\sigma ,\sigma ^{})`$
$`+`$ $`(\varphi _{i,\sigma }\varphi _{j,\sigma ^{}}|{\displaystyle \frac{e^2}{r}}|\varphi _{k,\sigma ^{}}\varphi _{l,\sigma })\delta \overline{\rho }^{(2)}(i,j,k,l;\sigma ,\sigma ^{})\}`$
$`+`$ $`{\displaystyle }d^3r\{{\displaystyle }d^3r^{}{\displaystyle \frac{e^2}{|𝐫𝐫^{}|}}n(𝐫^{})+{\displaystyle \frac{\delta E_{\mathrm{rxc}}}{\delta n(𝐫)}}+v_{\mathrm{ext}}(𝐫)`$
$`v_{\mathrm{eff}}(𝐫)\}\delta n(𝐫).`$
This equation tells us that the energy functional becomes stationary with respect to the density variation, if a next equation is satisfied for the solution of the auxiliary problem.

$`{\displaystyle \underset{i,\sigma }{}}\epsilon _i\delta \rho (i,i,\sigma )`$ (4.15)
$`+`$ $`{\displaystyle \underset{ij,\sigma }{}}(\epsilon _i\epsilon _j){\displaystyle d^3r\varphi _{i,\sigma }^{}(𝐫)\delta \varphi _{j,\sigma }(𝐫)\rho (i,j,\sigma )}`$
$`+`$ $`{\displaystyle \underset{i,j,k,l,\sigma ,\sigma ^{}}{\overset{(2)}{}}}{\displaystyle \frac{1}{2}}\{\delta (\varphi _{i,\sigma }\varphi _{j,\sigma ^{}}|{\displaystyle \frac{e^2}{r}}|\varphi _{k,\sigma ^{}}\varphi _{l,\sigma })\overline{\rho }^{(2)}(i,j,k,l;\sigma ,\sigma ^{})`$
$`+`$ $`(\varphi _{i,\sigma }\varphi _{j,\sigma ^{}}|{\displaystyle \frac{e^2}{r}}|\varphi _{k,\sigma ^{}}\varphi _{l,\sigma })\delta \overline{\rho }^{(2)}(i,j,k,l;\sigma ,\sigma ^{})\}=0`$
As shown in the next subsections, this condition for the density matrices is satisfied, if $`|\mathrm{\Psi }`$ is a solution of a many-body problem given by eq. (2.8) with properly chosen effective interactions. In other words, the present stationary condition on $`E[n]`$ is guaranteed provided that we can determine effective interactions such that the solution of eq. (2.8) satisfies eq. (4.15). In the following subsections, we show existence of required $`v_{\mathrm{eff}}(𝐫)`$, $`V_{\mathrm{eff}}^{(1)}(i,j,\sigma )`$, $`V_{\mathrm{eff}}^{(2)}(i,j,k,l;\sigma ,\sigma ^{})`$ and $`\widehat{H}_{\mathrm{rxc}}`$.

### 4.3 Variation of the Second Auxiliary Equation

Now we investigate variation of eq. (2.8). We note that not only the solution $`|\mathrm{\Psi }`$ but also all parameters, i.e. the orbital energy, interaction parameters and $`\widehat{H}_{\mathrm{rxc}}`$, are varied.

$`\{{\displaystyle \underset{i,\sigma }{}}(\epsilon _i+\delta \epsilon )c_{i,\sigma }^{}c_{i,\sigma }`$
$`+{\displaystyle \underset{i,j,\sigma }{}}(V_{\mathrm{eff}}^{(1)}(i,j,\sigma )+\delta V_{\mathrm{eff}}^{(1)}(i,j,\sigma ))c_{i,\sigma }^{}c_{j,\sigma }`$
$`+{\displaystyle \underset{i,j,k,l,\sigma ,\sigma ^{}}{\overset{(2)}{}}}(V_{\mathrm{eff}}^{(2)}(i,j,k,l;\sigma ,\sigma ^{})+\delta V_{\mathrm{eff}}^{(2)}(i,j,k,l;\sigma ,\sigma ^{}))`$
$`\times c_{i,\sigma }^{}c_{j,\sigma ^{}}^{}c_{k,\sigma ^{}}c_{l,\sigma }`$
$`+(\widehat{H}_{\mathrm{rxc}}+\delta \widehat{H}_{\mathrm{rxc}})\}|\mathrm{\Psi }+\delta \mathrm{\Psi }`$
$`=(E+\delta E)|\mathrm{\Psi }+\delta \mathrm{\Psi }.`$ (4.16)
As a result, we have the next equation which holds in the first order of variation.

$`\delta E`$ (4.17)
$`=`$ $`{\displaystyle \underset{i,\sigma }{}}(\delta \epsilon _i\rho (i,i,\sigma )+\epsilon _i\delta \rho (i,i,\sigma ))`$
$`+`$ $`{\displaystyle \underset{i,j,\sigma }{}}\left\{\delta V_{\mathrm{eff}}^{(1)}(i,j,\sigma )\rho (i,j,\sigma )+V_{\mathrm{eff}}^{(1)}(i,j,\sigma )\delta \rho (i,j,\sigma )\right\}`$
$`+`$ $`{\displaystyle \underset{i,j,k,l,\sigma ,\sigma ^{}}{\overset{(2)}{}}}\{\delta V_{\mathrm{eff}}^{(2)}(i,j,k,l;\sigma ,\sigma ^{})\rho ^{(2)}(i,j,k,l;\sigma ,\sigma ^{})`$
$`+V_{\mathrm{eff}}^{(2)}(i,j,k,l;\sigma ,\sigma ^{})\delta \rho ^{(2)}(i,j,k,l;\sigma ,\sigma ^{})\}`$
$`+`$ $`\delta \mathrm{\Psi }|\widehat{H}_{\mathrm{rxc}}|\mathrm{\Psi }`$
$`=`$ $`{\displaystyle \underset{i,\sigma }{}}\epsilon _i\delta \rho (i,i,\sigma )`$
$`+`$ $`{\displaystyle \underset{ij,\sigma }{}}(\epsilon _i\epsilon _j){\displaystyle d^3r\varphi _{i,\sigma }^{}(𝐫)\delta \varphi _{j,\sigma }(𝐫)\rho (i,j,\sigma )}`$
$`+`$ $`{\displaystyle \underset{i,j,k,l,\sigma ,\sigma ^{}}{\overset{(2)}{}}}\{\delta V_{\mathrm{eff}}^{(2)}(i,j,k,l;\sigma ,\sigma ^{})\overline{\rho }^{(2)}(i,j,k,l;\sigma ,\sigma ^{})`$
$`+`$ $`V_{\mathrm{eff}}^{(2)}(i,j,k,l;\sigma ,\sigma ^{})\delta \overline{\rho }^{(2)}(i,j,k,l;\sigma ,\sigma ^{})\}`$
$`+`$ $`{\displaystyle d^3r\delta v_{\mathrm{eff}}(𝐫)n(𝐫)}`$
$`+`$ $`{\displaystyle \underset{i,j,k,l,\sigma ,\sigma ^{}}{\overset{(2)}{}}}\{\delta V_{\mathrm{eff}}^{(2)}(i,j,k,l;\sigma ,\sigma ^{})\rho (i,l,\sigma )\rho (j,k,\sigma ^{})`$
$`+V_{\mathrm{eff}}^{(2)}(i,j,k,l;\sigma ,\sigma ^{})(\delta \rho (i,l,\sigma )\rho (j,k,\sigma ^{})`$
$`+\rho (i,l,\sigma )\delta \rho (j,k,\sigma ^{}))\}`$
$`+`$ $`{\displaystyle \underset{i,j,\sigma }{}}\left\{\delta V_{\mathrm{eff}}^{(1)}(i,j,\sigma )\rho (i,j,\sigma )+V_{\mathrm{eff}}^{(1)}(i,j,\sigma )\delta \rho (i,j,\sigma )\right\}`$
$`+`$ $`\delta \mathrm{\Psi }|\widehat{H}_{\mathrm{rxc}}|\mathrm{\Psi }.`$
Here we have utilized eq. (4.5) and eq. (4.6). For the stationary solution of eq. (2.8), $`\delta E=0`$ and we have another useful relation.

### 4.4 Determination of Effective Interactions

The expression of eq. (4.14) tells that natural expression for $`v_{\mathrm{eff}}(𝐫)`$ is similar to the KS effective potential and reads as,

$$v_{\mathrm{eff}}(𝐫)=d^3r\frac{e^2}{|𝐫𝐫^{}|}n(𝐫^{})+\frac{\delta E_{\mathrm{rxc}}}{\delta n(𝐫)}+v_{\mathrm{ext}}(𝐫).$$
(4.18)
Then, we have an expression of $`\delta v_{\mathrm{eff}}(𝐫)`$ as,

$$\delta v_{\mathrm{eff}}(𝐫)=d^3r\frac{e^2}{|𝐫𝐫^{}|}\delta n(𝐫^{})+\frac{\delta ^2E_{\mathrm{rxc}}}{\delta n(𝐫^{})\delta n(𝐫)}\delta n(𝐫^{}).$$
(4.19)
Besides, comparing eq. (4.14) and eq. (4.17), we notice that $`V_{\mathrm{eff}}^{(2)}(i,j,k,l;\sigma ,\sigma ^{})`$ should be chosen as,

$$V_{\mathrm{eff}}^{(2)}(i,j,k,l;\sigma ,\sigma ^{})=\frac{1}{2}(\varphi _{i,\sigma }\varphi _{j,\sigma ^{}}|\frac{e^2}{r}|\varphi _{k,\sigma ^{}}\varphi _{l,\sigma }).$$
(4.20)
Substitute eqs. (4.20) into eq. (4.17) and let $`\delta E=0`$, then we find that the required stationary condition, eq. (4.15), is satisfied, if the next condition holds.

$`{\displaystyle \underset{i,j,\sigma }{}}\left\{\delta V_{\mathrm{eff}}^{(1)}(i,j,\sigma )\rho (i,j,\sigma )+V_{\mathrm{eff}}^{(1)}(i,j,\sigma )\delta \rho (i,j,\sigma )\right\}`$ (4.21)
$`+\delta \mathrm{\Psi }|\widehat{H}_{\mathrm{rxc}}|\mathrm{\Psi }`$
$`=`$ $`{\displaystyle d^3r\delta v_{\mathrm{eff}}(𝐫)n(𝐫)}`$
$`{\displaystyle \underset{i,j,k,l,\sigma ,\sigma ^{}}{\overset{(2)}{}}}\{\delta V_{\mathrm{eff}}^{(2)}(i,j,k,l;\sigma ,\sigma ^{})\rho (i,l,\sigma )\rho (j,k,\sigma ^{})`$
$`+V_{\mathrm{eff}}^{(2)}(i,j,k,l;\sigma ,\sigma ^{})(\delta \rho (i,l,\sigma )\rho (j,k,\sigma ^{})`$
$`+\rho (i,l,\sigma )\delta \rho (j,k,\sigma ^{}))\}`$
This new condition is useful for determination of $`V_{\mathrm{eff}}^{(1)}(i,j,𝐫)`$ and $`\widehat{H}_{\mathrm{rxc}}`$. Substitute eq. (4.19) and eq. (4.20) into eq. (4.21), we obtain the next condition for $`V_{\mathrm{eff}}^{(1)}(i,j,\sigma )`$ and $`\widehat{H}_{\mathrm{rxc}}`$.

$`{\displaystyle \underset{i,j,\sigma }{}}\left\{\delta V_{\mathrm{eff}}^{(1)}(i,j,\sigma )\rho (i,j,\sigma )+V_{\mathrm{eff}}^{(1)}(i,j,\sigma )\delta \rho (i,j,\sigma )\right\}`$ (4.22)
$`+\delta \mathrm{\Psi }|\widehat{H}_{\mathrm{rxc}}|\mathrm{\Psi }`$
$`=`$ $`{\displaystyle \underset{i,j,\sigma }{}}{\displaystyle d^3rd^3r^{}\frac{e^2}{|𝐫𝐫^{}|}n(𝐫)}`$
$`\times (\delta \varphi _{i,\sigma }^{}(𝐫^{})\varphi _{j,\sigma }(𝐫^{})+\varphi _{i,\sigma }^{}(𝐫^{})\delta \varphi _{j,\sigma }(𝐫^{}))\rho (i,j,\sigma )`$
$``$ $`{\displaystyle \underset{i,j,\sigma }{}}{\displaystyle d^3rd^3r^{}\frac{e^2}{|𝐫𝐫^{}|}n(𝐫)\varphi _{i,\sigma }^{}(𝐫^{})\varphi _{j,\sigma }(𝐫^{})\delta \rho (i,j,\sigma )}`$
$``$ $`{\displaystyle \underset{i,j,\sigma }{\overset{\overline{(2)}}{}}}\{{\displaystyle \underset{k,l,\sigma ^{}}{\overset{\overline{(2)}}{}}}{\displaystyle \frac{1}{2}}\delta (\varphi _{i,\sigma }\varphi _{k,\sigma ^{}}|{\displaystyle \frac{e^2}{r}}|\varphi _{l,\sigma ^{}}\varphi _{j,\sigma })\rho (k,l,\sigma ^{})\rho (i,j,\sigma )`$
$`+`$ $`{\displaystyle \underset{k,l,\sigma ^{}}{\overset{\overline{(2)}}{}}}(\varphi _{i,\sigma }\varphi _{k,\sigma ^{}}|{\displaystyle \frac{e^2}{r}}|\varphi _{l,\sigma ^{}}\varphi _{j,\sigma })\rho (k,l,\sigma ^{})\delta \rho (i,j,\sigma )\}`$
$``$ $`{\displaystyle d^3rd^3r^{}\frac{\delta ^2E_{\mathrm{rxc}}}{\delta n(𝐫^{})\delta n(𝐫)}n(𝐫)\delta n(𝐫^{})}.`$
Here, a summation, $`^{\overline{(2)}}`$, represents a restricted summation satisfying $`_{i,j,\sigma }^{\overline{(2)}}_{k,l,\sigma ^{}}^{\overline{(2)}}=_{i,k,l,j,\sigma ,\sigma ^{}}^{(2)}`$. The condition, eq. (4.22), is satisfied, if we choose $`V_{\mathrm{eff}}^{(1)}(i,j,𝐫)`$ and $`\widehat{H}_{\mathrm{rxc}}`$ as follows.

$`V_{\mathrm{eff}}^{(1)}(i,j,\sigma )`$ (4.23)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle d^3rd^3r^{}\frac{e^2}{|𝐫𝐫^{}|}n(𝐫)\varphi _{i,\sigma }^{}(𝐫^{})\varphi _{j,\sigma }(𝐫^{})}`$
$``$ $`{\displaystyle \underset{i_0,j_0,\sigma _0}{\overset{\overline{(2)}}{}}}\delta _{i,i_0}\delta _{j,j_0}\delta _{\sigma ,\sigma _0}`$
$`\times `$ $`{\displaystyle \underset{k,l,\sigma ^{}}{\overset{\overline{(2)}}{}}}{\displaystyle \frac{1}{2}}(\varphi _{i,\sigma }\varphi _{k,\sigma ^{}}|{\displaystyle \frac{e^2}{r}}|\varphi _{l,\sigma ^{}}\varphi _{j,\sigma })\rho (k,l,\sigma ^{}),`$
$`\mathrm{\Psi }|\widehat{H}_{\mathrm{rxc}}|\mathrm{\Psi }`$ $`=`$ $`{\displaystyle d^3r\frac{\delta E_{\mathrm{rxc}}}{\delta n(𝐫)}n(𝐫)}+E_{\mathrm{rxc}}[n],`$ (4.24)
$`\delta \mathrm{\Psi }|\widehat{H}_{\mathrm{rxc}}|\mathrm{\Psi }`$ $`=`$ $`{\displaystyle d^3rd^3r^{}\frac{\delta ^2E_{\mathrm{rxc}}}{\delta n(𝐫^{})\delta n(𝐫)}n(𝐫)\delta n(𝐫^{})}.`$
Here we note that required conditions for $`\widehat{H}_{\mathrm{rxc}}`$ are only eq. (4.24) and eq. (LABEL:exp-dhrxc). These equations are satisfied, if we choose an explicit expression of $`\widehat{H}_{\mathrm{rxc}}`$ as,

$`\widehat{H}_{\mathrm{rxc}}`$ $`=`$ $`\epsilon _{\mathrm{rxc}}[n]\widehat{I},`$ (4.26)
$`\epsilon _{\mathrm{rxc}}[n]`$ $`=`$ $`{\displaystyle d^3r\frac{\delta E_{\mathrm{rxc}}}{\delta n(𝐫)}n(𝐫)}+E_{\mathrm{rxc}}[n].`$ (4.27)
The operator, $`\widehat{I}`$, represents an identity operator in the phase space of the $`N`$-particle system under discussion. Finding of these expressions for effective interactions is an important result of the present theory.

## 5 The System of Auxiliary Equations

We have obtained a new set of equations, which gives a generalization of KSE. The system is summarized as follows.

$$\left\{\frac{\mathrm{}^2}{2m}\mathrm{\Delta }_𝐫+v_{\mathrm{eff}}(𝐫)\right\}\varphi _{i,\sigma }(𝐫)=\epsilon _i\varphi _{i,\sigma }(𝐫),$$
(5.1)
$`\{{\displaystyle \underset{i,\sigma }{}}\epsilon _ic_{i,\sigma }^{}c_{i,\sigma }+{\displaystyle \underset{i,j,\sigma }{}}V_{\mathrm{eff}}^{(1)}(i,j,\sigma )c_{i,\sigma }^{}c_{j,\sigma }`$ (5.2)
$`+`$ $`{\displaystyle \underset{i,j,k,l,\sigma ,\sigma ^{}}{\overset{(2)}{}}}V_{\mathrm{eff}}^{(2)}(i,j,k,l;\sigma ,\sigma ^{})c_{i,\sigma }^{}c_{j,\sigma ^{}}^{}c_{k,\sigma ^{}}c_{l,\sigma }`$
$`+`$ $`\epsilon _{\mathrm{rxc}}\widehat{I}\}|\mathrm{\Psi }=E|\mathrm{\Psi },`$
$$v_{\mathrm{eff}}(𝐫)=d^3r\frac{e^2}{|𝐫𝐫^{}|}n(𝐫^{})+\frac{\delta E_{\mathrm{rxc}}}{\delta n(𝐫)}+v_{\mathrm{ext}}(𝐫),$$
(5.3)
$$V_{\mathrm{eff}}^{(2)}(i,j,k,l;\sigma ,\sigma ^{})=\frac{1}{2}(\varphi _{i,\sigma }\varphi _{j,\sigma ^{}}|\frac{e^2}{r}|\varphi _{k,\sigma ^{}}\varphi _{l,\sigma }),$$
(5.4)
$`V_{\mathrm{eff}}^{(1)}(i,j,\sigma )`$ (5.5)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle d^3rd^3r^{}\frac{e^2}{|𝐫𝐫^{}|}n(𝐫)\varphi _{i,\sigma }^{}(𝐫^{})\varphi _{j,\sigma }(𝐫^{})}`$
$``$ $`{\displaystyle \underset{i_0,j_0,\sigma _0}{\overset{\overline{(2)}}{}}}\delta _{i,i_0}\delta _{j,j_0}\delta \sigma ,\sigma _0`$
$`\times `$ $`{\displaystyle \underset{k,l,\sigma ^{}}{\overset{\overline{(2)}}{}}}{\displaystyle \frac{1}{2}}(\varphi _{i,\sigma }\varphi _{k,\sigma ^{}}|{\displaystyle \frac{e^2}{r}}|\varphi _{l,\sigma ^{}}\varphi _{j,\sigma })\rho (k,l,\sigma ^{}),`$
$`\epsilon _{\mathrm{rxc}}[n]`$ $`=`$ $`E_{\mathrm{rxc}}[n]{\displaystyle d^3r\frac{\delta E_{\mathrm{rxc}}}{\delta n(𝐫)}n(𝐫)}.`$ (5.6)
In the following subsections, we discuss interpretation of the new equations and possible calculational strategies.

### 5.1 Interpretation of the Equations

Each equation in the new set of auxiliary equations has its own meaning. We review important properties.

Eq. (5.1) is a direct generalization of KSE. The effective single-particle potential, eq. (5.3), is comprised of a Hartree term, a new exchange-correlation potential and the external potential. However, we have to comment that $`n(𝐫)`$ is given by a multi Slater determinant, $`|\mathrm{\Psi }_{\mathrm{GS}}`$, and definition of $`E_{\mathrm{rxc}}`$ is different from that utilized in the original KS scheme.

Eq. (5.2) gives an auxiliary many-body problem. The summation, $`^{(2)}`$, is arbitrary as far as $`\widehat{H}_{\mathrm{eff}}`$ is Hermite. If $`^{(2)}`$ is taken for all possible combinations, $`E_{\mathrm{rxc}}0`$ and the problem becomes identical to the original Coulomb-potential problem for a non-uniform electron gas. If the summation is not taken for any combination, our system of equations reduces to KSE. Thus, our new equations is a natural extension of KSE which makes a connection between the usual KS scheme and the Coulomb problem.

We have obtained eq. (5.4) from observation of eq. (4.14) and eq. (4.17). A natural question may be “Why does any screened interaction appear?” This point is discussed in §6, too.

Eq. (5.5) is naturally interpreted as a counter term required by definitions of eq. (5.3) and eq. (5.4). This term represents a Hartree-like mean-field. However, we have to be careful with appearance of the second term with a summation $`^{\overline{(2)}}`$.

Eq. (5.6) is a correction coming from the reduced exchange correlation functional. This expression indicates that $`\epsilon _{\mathrm{rxc}}`$ is a functional of $`n(𝐫)`$.

### 5.2 Calculational Strategy

If we assume that we can solve eq. (5.1) and eq. (5.2), a process to obtain the ground state will be the following. Here, we also assume that we know a value of $`E_{\mathrm{rxc}}[n]`$ and its variational derivative for any $`n(𝐫)`$.

1. Start from an initial single-particle density $`n(𝐫)`$.
2. Once $`n(𝐫)`$ is given, $`v_{\mathrm{eff}}(𝐫)`$ is given.
3. Solve eq. (5.1) to obtain $`\varphi _{i,\sigma }(𝐫)`$. Here $`n(𝐫)`$ is assumed to be fixed.
4. Solve eq. (5.2) by a proper method. Since $`V_{\mathrm{eff}}^{(1)}(i,j,\sigma )`$ and $`\widehat{H}_{\mathrm{rxc}}`$ are given by the density matrices and $`n(𝐫)`$ of the solution, the determining equation has to be solved in a self-consistent manner.
5. Recalculate $`n(𝐫)`$ and goto the step 2 until convergence in $`n(𝐫)`$ is obtained.

The above optimization process is done in a phase space of $`n(𝐫)`$ given by $`|\mathrm{\Psi }`$.

## 6 Interpretation of New Energy Functionals

As discussed in §3, we considered that functionals, $`T_0[n]`$, $`V_{\mathrm{ee}}[n]`$ and $`E_{\mathrm{rxc}}[n]`$, are density functionals. We discuss this important point in more details. We hope that the explanation below will answer remaining questions for our formalism of DFT.

Let us review the HT theory at first. To formulate rigorous scheme of KSE, they considered a functional of a single Slater determinant $`|\varphi `$. Here we use a notation, $`n_\varphi (𝐫)=\varphi |\widehat{n}(𝐫)|\varphi `$.

$`G_H[\varphi ]`$ $`=`$ $`\varphi |\widehat{T}|\varphi +\mathrm{\Delta }T(\varphi )+{\displaystyle d^3rv_{\mathrm{ext}}(𝐫)n_\varphi (𝐫)}`$ (6.1)
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle d^3rd^3r^{}\frac{e^2}{|𝐫𝐫^{}|}n_\varphi (𝐫)n_\varphi (𝐫^{})}`$
$`+`$ $`E_{\mathrm{xc}}(\varphi ).`$
To define $`\mathrm{\Delta }T(\varphi )`$, we are required to determine a multi Slater determinant $`|\psi `$ and a single Slater determinant $`|\varphi ^{}`$ via two minimization processes given by next two functionals.

$`G_{T+V}[\varphi ]`$ $`=`$ $`\underset{\psi n_\varphi }{\mathrm{min}}\psi |\widehat{T}+\widehat{V}\psi ,`$ (6.2)
$`G_T[\varphi ]`$ $`=`$ $`\underset{\varphi ^{}n_\varphi }{\mathrm{min}}\varphi ^{}|\widehat{T}\varphi ^{}.`$ (6.3)
Using $`|\psi `$ and $`|\varphi ^{}`$, $`\mathrm{\Delta }T`$ is defined as,

$$\mathrm{\Delta }T[\varphi ]=\psi |\widehat{T}|\psi \varphi ^{}|\widehat{T}|\varphi ^{}.$$
(6.4)
Another functional $`E_{\mathrm{xc}}(\varphi )`$ is given as,

$$E_{\mathrm{xc}}[\varphi ]=\psi |\widehat{V}|\psi \frac{1}{2}d^3rd^3r^{}\frac{e^2}{|𝐫𝐫^{}|}n_\varphi (𝐫)n_\varphi (𝐫^{}).$$
(6.5)
Since $`G_H`$ is rewritten as,

$`G_H[\varphi ]`$ $`=`$ $`\varphi |\widehat{T}|\varphi \varphi ^{}|\widehat{T}|\varphi ^{}`$ (6.6)
$`+`$ $`\psi |\widehat{T}+\widehat{V}|\psi +{\displaystyle d^3rv_{\mathrm{ext}}(𝐫)n_\varphi (𝐫)},`$
it is easily shown that the minimization of $`G_H[\varphi ]`$ with respect to $`|\varphi `$ is equivalent to obtain the true ground-state energy, $`E_0`$. Namely,

$$\underset{\varphi }{\mathrm{min}}G_H[\varphi ]=\underset{n}{\mathrm{min}}\left\{F[n]+d^3rv_{\mathrm{ext}}(𝐫)n(𝐫)\right\}=E_0,$$
(6.7)
where $`F[n]`$ is Levy’s DFT. The minimization of $`G_H`$ with respect to $`|\varphi `$ is shown to be identical to solve KSE.$`^{\text{?)}}`$

The strategy given in the present paper is to utilize a multi Slater determinant $`|\mathrm{\Psi }`$ in place of $`|\varphi `$ and replace $`\widehat{T}`$ with $`\widehat{T}+\widehat{V}_{\mathrm{reduced}}`$ defined below. Namely, our extension of KSE is rigorously given by the following functional of $`|\mathrm{\Psi }`$. We use a notation, $`n_\mathrm{\Psi }(𝐫)=\mathrm{\Psi }|\widehat{n}|\mathrm{\Psi }`$, too.

$`\overline{G}_H[\mathrm{\Psi }]`$ (6.8)
$`=`$ $`\mathrm{\Psi }|\widehat{T}+\widehat{V}_{\mathrm{red}}|\mathrm{\Psi }\underset{\mathrm{\Psi }^{}n_\mathrm{\Psi }}{\mathrm{min}}\mathrm{\Psi }^{}|\widehat{T}+\widehat{V}_{\mathrm{red}}\mathrm{\Psi }^{}`$
$`+`$ $`F[n_\mathrm{\Psi }]+{\displaystyle d^3rv_{\mathrm{ext}}(𝐫)n_\mathrm{\Psi }(𝐫)}.`$
Here a reduced interaction operator, $`\widehat{V}_{\mathrm{red}}`$, is defined as a following operation.

$`\mathrm{\Psi }|\widehat{V}_{\mathrm{red}}|\mathrm{\Psi }`$ (6.9)
$`=`$ $`{\displaystyle \underset{i,j,k,l,\sigma ,\sigma ^{}}{\overset{(2)}{}}}{\displaystyle \frac{1}{2}}{\displaystyle d^3rd^3r^{}\frac{e^2}{|𝐫𝐫^{}|}}`$
$`\times \varphi _{i,\sigma }^{}(𝐫)\varphi _{j,\sigma ^{}}^{}(𝐫^{})\varphi _{k,\sigma ^{}}(𝐫^{})\varphi _{l,\sigma }(𝐫)`$
$`\times \{\mathrm{\Psi }|c_{i,\sigma }^{}c_{j,\sigma ^{}}^{}c_{k,\sigma ^{}}c_{l,\sigma ^{}}|\mathrm{\Psi }`$
$`\mathrm{\Psi }|c_{i,\sigma }^{}c_{l,\sigma ^{}}|\mathrm{\Psi }\mathrm{\Psi }|c_{j,\sigma ^{}}^{}c_{k,\sigma ^{}}|\mathrm{\Psi }\}.`$
Then, we can rewrite $`\overline{G}_H`$ as,

$`\overline{G}_H[\mathrm{\Psi }]`$ (6.10)
$`=`$ $`\mathrm{\Psi }|\widehat{T}+\widehat{V}_{\mathrm{red}}|\mathrm{\Psi }+{\displaystyle d^3rv_{\mathrm{ext}}(𝐫)n_\mathrm{\Psi }(𝐫)}`$
$`+{\displaystyle \frac{1}{2}}{\displaystyle d^3rd^3r^{}\frac{e^2}{|𝐫𝐫^{}|}n_\mathrm{\Psi }(𝐫)n_\mathrm{\Psi }(𝐫^{})}`$
$`+E_{\mathrm{rxc}}[n_\mathrm{\Psi }].`$
Here $`E_{\mathrm{rxc}}[n_\mathrm{\Psi }]`$ is defined as,

$`E_{\mathrm{rxc}}[n_\mathrm{\Psi }]`$ (6.11)
$`=`$ $`F[n_\mathrm{\Psi }]\underset{\mathrm{\Psi }^{}n_\mathrm{\Psi }}{\mathrm{min}}\mathrm{\Psi }^{}|\widehat{T}+\widehat{V}_{\mathrm{red}}\mathrm{\Psi }^{}`$
$``$ $`{\displaystyle \frac{1}{2}}{\displaystyle d^3rd^3r^{}\frac{e^2}{|𝐫𝐫^{}|}n_\mathrm{\Psi }(𝐫)n_\mathrm{\Psi }(𝐫^{})}.`$
This definition of $`E_{\mathrm{rxc}}`$ is equivalent to eq. (3.3), if we identify $`n_\mathrm{\Psi }(𝐫)`$ with $`n(𝐫)`$. It is straight-forward to derive our set of auxiliary equations by a variational principle on functionals $`\overline{G}_H`$.

These expressions tell a lot about possibility of our extensions. Eq. (6.9) requires to determine not only range of $`^{(2)}`$ but also KS orbitals, $`\varphi _{i,\sigma }(𝐫)`$. However, we may utilize any localized orbitals in the definition of $`\widehat{V}_{\mathrm{red}}`$. Since eq. (6.9) takes a form of a correction to the mean-field approximation, we should include a key process, e.g. interactions between electrons in the same localized orbital, to reduce amount of $`E_{\mathrm{rxc}}`$. Besides, we might be allowed to utilize a screened interaction in eq. (6.9), if $`E_{\mathrm{rxc}}`$ became a small correction.

Note that definition of $`\widehat{V}_{\mathrm{red}}`$ requires to fix interaction terms defined by $`\varphi _{i,\sigma }(𝐫)`$. In the previous derivation, we implicitly assumed an optimization process of orbitals, $`\varphi _{i,\sigma }(𝐫)`$, too. Once $`\varphi _{i,\sigma }(𝐫)`$ is varied, values of $`E_{\mathrm{rxc}}`$ change in principle. In several expected situation like the valence fluctuation, however, the Hartree-term would be dominant to determine the shape of local orbitals. Thus, to approximate $`E_{\mathrm{rxc}}`$ independent of a small change of $`\varphi _{i,\sigma }(𝐫)`$ might be practical.

## 7 Possible Approximation

For a practical calculation, we have to decide an approximate evaluation of $`E_{\mathrm{rxc}}`$ as the original KS scheme. However, we have to consider that $`E_{\mathrm{rxc}}`$ depends on the summation, $`^{(2)}`$, in eq. (5.2). As a case study, let us consider a localized f-orbital embedded in conduction electrons. In this discussion, we omit the relativistic effects for simplicity.

Using a proper unitary transformation, a Hubbard-type local interaction between f-electrons can be described explicitly by 2-body interaction terms in eq. (5.2). Let us assume that other 2-body interactions are excluded in $`^{(2)}`$. This effective Hamiltonian corresponds to an Anderson impurity Hamiltonian,$`^{\text{?}\text{)}}`$ where both the f-orbital and conduction states are given by solutions of eq. (5.1). The residual exchange-correlation term, $`E_{\mathrm{rxc}}`$, represents not only all of the exchange-correlation effects among conduction electrons but also a part of exchange-correlation effects between an f-electron and conduction electrons. The latter exchange-correlation includes an effect by charge fluctuation on the f-level, which introduce effective expansion or contraction of f-charge affecting the effective potential, $`v_{\mathrm{eff}}(𝐫)`$, too. However, we can expect that most of the effects coming from this charge redistribution can be represented by the Hartree term in $`v_{\mathrm{eff}}(𝐫)`$.

Thus, we can expect that $`E_{\mathrm{rxc}}`$ is not large in amount, as far as we include the most important exchange-correlation processes in our auxiliary equations through the local interaction, $`V_{\mathrm{eff}}^{(2)}(i,j,k,l;\sigma ,\sigma ^{})`$. The simplest practical approximation for $`E_{\mathrm{rxc}}`$ is a kind of the local density approximation (LDA). Many strategy of LDA for $`E_{\mathrm{rxc}}`$ would be possible. However, to determine this function, we have to decide inputs. If we adopt representation by localized orbitals,$`^{\text{?)}}`$ we can define a single particle density of f-electrons, $`n_\mathrm{f}(𝐫)`$, and that of conduction electrons, $`n_\mathrm{c}(𝐫)`$. Then, some versions of LDA would be classified by the inputs as, 1) $`n(𝐫)`$ at $`𝐫`$, 2) $`n_\mathrm{f}(𝐫)`$ and $`n_\mathrm{c}(𝐫)`$ at $`𝐫`$, 3) $`n_\mathrm{f}(𝐫)`$, $`n_\mathrm{c}(𝐫)`$, $`\mathrm{\Delta }n_\mathrm{f}(𝐫)`$ and $`\mathrm{\Delta }n_\mathrm{c}(𝐫)`$ at $`𝐫`$. The category 3) corresponds to GGA or more closely to meta-GGA.$`^{\text{?}\text{)}}`$ Determination of a practical method is a remaining important problem.

Reliable $`E_{\mathrm{rxc}}`$ should be determined by accurate numerical calculations. As has been done for determination of LDA using the diffusion Monte-Carlo (DMC) data$`^{\text{?}\text{}\text{?}\text{)}}`$ for the uniform electron gas,$`^{\text{?)}}`$ parameterization might be possible if there were data for an impurity problem by DMC or other accurate methods.

Another possibility is to adopt the optimized effective potential method with exact evaluation of the exchange interaction (EXX).$`^{\text{?}\text{}\text{?}\text{}\text{?}\text{)}}`$ For example, if evaluation of $`E_{\mathrm{rxc}}`$ by EXX with RPA for correlation energy (EXX+RPA)$`^{\text{?}\text{)}}`$ is tractable for evaluation of $`E_{\mathrm{rxc}}`$, and if the short-range dynamical correlation described in $`\widehat{H}_{\mathrm{eff}}`$ is well-described by DMF, a hybrid method of DMF+EXX+RPA might be possible.

## 8 Discussion

Starting from Levy’s DFT, we have formulate an extension of the Kohn-Sham equations in order to simulate correlation effects coming from localized nature of electrons in materials. The system of equations consists of 1) determination equations of single-particle orbitals and 2) a many-body model problem. These auxiliary equations should be determined self-consistently. In §6, we have given a rigorous derivation of the extended KSE.

As commented in §5.1, our new scheme provides a connection of the original Coulomb problem and KSE. This property is very useful for application. Conversely, a serious problem is how to determine $`^{(2)}`$ most efficiently. Without ruining qualitative accuracy, we have to control this summation to optimize simplicity for actual calculations. The problem, we think, may be solved by knowledge of information science.

From a practical point of view, our theory may lead to a simple approximation method including non-local correlation effects, which cannot be treated even by GGA.$`^{\text{?}\text{)}}`$ Comparison between our formalism and meta-GGA$`^{\text{?)}}`$ might be an interesting problem. Since our strategy supplies a kind of the so-called first-principles calculation, we hope that the present scheme may be a step toward a unified theory for the interacting electron systems.

Eq. (5.2) can be interpreted as a kind of the extended Hubbard model. Although 2-body interactions appearing in eq. (5.2) is bare Coulomb (and exchange) interactions, the model system can describe every magnetic systems as well as non-magnetic materials. Namely, if we properly solve the auxiliary problem, spontaneous symmetry breaking can be discussed in a realistic electronic-state calculation. Besides, if a superconducting state was realized only by electron-electron repulsion, our strategy could be utilized to explore the superconductivity by a DFT-based calculation.

As well-known, the relativistic DFT (RDFT),$`^{\text{?}\text{}\text{?}\text{}\text{?}\text{}\text{?}\text{}\text{?}\text{)}}`$ the current density functional theory (CDFT)$`^{\text{?}\text{}\text{?}\text{)}}`$ and the unified version, i.e. the relativistic current density functional theory (RCDFT)$`^{\text{?}\text{}\text{?}\text{)}}`$ have been established. There should hopefully appear an extended version of the present theory based on RCDFT. Besides, there are many remaining important directions to be investigated, e.g. an extension of KSE for correlated systems at finite-temperatures$`^{\text{?)}}`$ and DFT-description of time-dependent phenomena.$`^{\text{?, }\text{?}\text{}\text{?}\text{}\text{?}\text{)}}`$

## Acknowledgements

The author would like to thank fruitful discussions and encouragements given by Dr. T. Kotani, Dr. E. Tamura and Dr. Y. Morikawa in the early stage of the project. He is grateful for important comments and discussions given by Prof. S. Tsuneyuki and Prof. M. Imada. The present work has been supported by Grant-in-Aid from Ministry of Education, Culture, Sports, Science and Technology Japan.