We consider a system of electrons in infinite matter, the so-called electron gas. This is a homogeneous system and the one-particle states are given by plane wave function normalized to a volume \( V \) for a box with length \( L \) (the limit \( L\rightarrow \infty \) is to be taken after we have computed various expectation values)
$$ \psi_{\mathbf{k}\sigma}(\mathbf{r})= \frac{1}{\sqrt{V}}\exp{(i\mathbf{kr})}\xi_{\sigma} $$where \( \mathbf{k} \) is the wave number and \( \xi_{\sigma} \) is a spin function for either spin up or down
$$ \xi_{\sigma=+1/2}=\left(\begin{array}{c} 1 \\ 0 \end{array}\right) \hspace{0.5cm} \xi_{\sigma=-1/2}=\left(\begin{array}{c} 0 \\ 1 \end{array}\right). $$We assume that we have periodic boundary conditions which limit the allowed wave numbers to
$$ k_i=\frac{2\pi n_i}{L}\hspace{0.5cm} i=x,y,z \hspace{0.5cm} n_i=0,\pm 1,\pm 2, \dots $$We assume first that the particles interact via a central, symmetric and translationally invariant interaction \( V(r_{12}) \) with \( r_{12}=|\mathbf{r}_1-\mathbf{r}_2| \). The interaction is spin independent.
The total Hamiltonian consists then of kinetic and potential energy
$$ \hat{H} = \hat{T}+\hat{V}. $$The operator for the kinetic energy is given by
$$ \hat{T}=\sum_{\mathbf{k}\sigma}\frac{\hbar^2k^2}{2m}a_{\mathbf{k}\sigma}^{\dagger}a_{\mathbf{k}\sigma}. $$Find the expression for the interaction \( \hat{V} \) expressed with creation and annihilation operators.
The expression for the interaction has to be written in \( k \) space, even though \( V \) depends only on the relative distance. It means that you need to set up the Fourier transform \( \langle \mathbf{k}_i\mathbf{k}_j| V | \mathbf{k}_m\mathbf{k}_n\rangle \).
A general two-body interaction element is given by (not using anti-symmetrized matrix elements)
$$ \hat{V} = \frac{1}{2} \sum_{pqrs} \langle pq \hat{v} \vert rs\rangle a_p^\dagger a_q^\dagger a_s a_r , $$where \( \hat{v} \) is assumed to depend only on the relative distance between two interacting particles, that is \( \hat{v} = v(\vec r_1, \vec r_2) = v(|\vec r_1 - \vec r_2|) = v(r) \), with \( r = |\vec r_1 - \vec r_2| \)).
In our case we have, writing out explicitely the spin degrees of freedom as well
$$ \begin{equation} \hat{V} = \frac{1}{2} \sum_{\substack{\sigma_p \sigma_q \\ \sigma_r \sigma_s}} \sum_{\substack{\mathbf{k}_p \mathbf{k}_q \\ \mathbf{k}_r \mathbf{k}_s}} \langle \mathbf{k}_p \sigma_p, \mathbf{k}_q \sigma_2\vert v \vert \mathbf{k}_r \sigma_3, \mathbf{k}_s \sigma_s\rangle a_{\mathbf{k}_p \sigma_p}^\dagger a_{\mathbf{k}_q \sigma_q}^\dagger a_{\mathbf{k}_s \sigma_s} a_{\mathbf{k}_r \sigma_r} . \label{_auto1} \end{equation} $$Inserting plane waves as eigenstates we can rewrite the matrix element as
$$ \langle \mathbf{k}_p \sigma_p, \mathbf{k}_q \sigma_q\vert \hat{v} \vert \mathbf{k}_r \sigma_r, \mathbf{k}_s \sigma_s\rangle = \frac{1}{V^2} \delta_{\sigma_p \sigma_r} \delta_{\sigma_q \sigma_s} \int\int \exp{-i(\mathbf{k}_p \cdot \mathbf{r}_p)} \exp{-i( \mathbf{k}_q \cdot \mathbf{r}_q)} \hat{v}(r) \exp{i(\mathbf{k}_r \cdot \mathbf{r}_p)} \exp{i( \mathbf{k}_s \cdot \mathbf{r}_q)} d\mathbf{r}_p d\mathbf{r}_q , $$where we have used the orthogonality properties of the spin functions. We change now the variables of integration by defining \( \mathbf{r} = \mathbf{r}_p - \mathbf{r}_q \), which gives \( \mathbf{r}_p = \mathbf{r} + \mathbf{r}_q \) and \( d^3 \mathbf{r} = d^3 \mathbf{r}_p \).
The limits are not changed since they are from \( -\infty \) to \( \infty \) for all integrals. This results in
$$ \begin{align*} \langle \mathbf{k}_p \sigma_p, \mathbf{k}_q \sigma_q\vert \hat{v} \vert \mathbf{k}_r \sigma_r, \mathbf{k}_s \sigma_s\rangle &= \frac{1}{V^2} \delta_{\sigma_p \sigma_r} \delta_{\sigma_q \sigma_s} \int\exp{i (\mathbf{k}_s - \mathbf{k}_q) \cdot \mathbf{r}_q} \int v(r) \exp{i(\mathbf{k}_r - \mathbf{k}_p) \cdot ( \mathbf{r} + \mathbf{r}_q)} d\mathbf{r} d\mathbf{r}_q \\ &= \frac{1}{V^2} \delta_{\sigma_p \sigma_r} \delta_{\sigma_q \sigma_s} \int v(r) \exp{i\left[(\mathbf{k}_r - \mathbf{k}_p) \cdot \mathbf{r}\right]} \int \exp{i\left[(\mathbf{k}_s - \mathbf{k}_q + \mathbf{k}_r - \mathbf{k}_p) \cdot \mathbf{r}_q\right]} d\mathbf{r}_q d\mathbf{r} . \end{align*} $$We recognize the integral over \( \mathbf{r}_q \) as a \( \delta \)-function, resulting in
$$ \langle \mathbf{k}_p \sigma_p, \mathbf{k}_q \sigma_q\vert \hat{v} \vert \mathbf{k}_r \sigma_r, \mathbf{k}_s \sigma_s\rangle = \frac{1}{V} \delta_{\sigma_p \sigma_r} \delta_{\sigma_q \sigma_s} \delta_{(\mathbf{k}_p + \mathbf{k}_q),(\mathbf{k}_r + \mathbf{k}_s)} \int v(r) \exp{i\left[(\mathbf{k}_r - \mathbf{k}_p) \cdot \mathbf{r}\right]} d^3r . $$For this equation to be different from zero, we must have conservation of momenta, we need to satisfy \( \mathbf{k}_p + \mathbf{k}_q = \mathbf{k}_r + \mathbf{k}_s \).
We can use the conservation of momenta to remove one of the summation variables resulting in
$$ \hat{V} = \frac{1}{2V} \sum_{\sigma \sigma'} \sum_{\mathbf{k}_p \mathbf{k}_q \mathbf{k}_r} \left[ \int v(r) \exp{i\left[(\mathbf{k}_r - \mathbf{k}_p) \cdot \mathbf{r}\right]} d^3r \right] a_{\mathbf{k}_p \sigma}^\dagger a_{\mathbf{k}_q \sigma'}^\dagger a_{\mathbf{k}_p + \mathbf{k}_q - \mathbf{k}_r, \sigma'} a_{\mathbf{k}_r \sigma}, $$which can be rewritten as
$$ \begin{equation} \hat{V} = \frac{1}{2V} \sum_{\sigma \sigma'} \sum_{\mathbf{k} \mathbf{p} \mathbf{q}} \left[ \int v(r) \exp{-i( \mathbf{q} \cdot \mathbf{r})} d\mathbf{r} \right] a_{\mathbf{k} + \mathbf{q}, \sigma}^\dagger a_{\mathbf{p} - \mathbf{q}, \sigma'}^\dagger a_{\mathbf{p} \sigma'} a_{\mathbf{k} \sigma}, \label{eq:V} \end{equation} $$In the last equation we defined the quantities \( \mathbf{p} = \mathbf{k}_p + \mathbf{k}_q - \mathbf{k}_r \), \( \mathbf{k} = \mathbf{k}_r \) og \( \mathbf{q} = \mathbf{k}_p - \mathbf{k}_r \).
Let us now compute the expectation value of the reference energy using the expressions for the kinetic energy operator and the interaction. We need to compute \( \langle \Phi_0\vert \hat{H} \vert \Phi_0\rangle = \langle \Phi_0\vert \hat{T} \vert \Phi_0\rangle + \langle \Phi_0\vert \hat{V} \vert \Phi_0\rangle \), where \( \vert \Phi_0\rangle \) is our reference Slater determinant, constructed from filling all single-particle states up to the Fermi level. Let us start with the kinetic energy first
$$ \langle \Phi_0\vert \hat{T} \vert \Phi_0\rangle = \langle \Phi_0\vert \left( \sum_{\mathbf{p} \sigma} \frac{\hbar^2 p^2}{2m} a_{\mathbf{p} \sigma}^\dagger a_{\mathbf{p} \sigma} \right) \vert \Phi_0\rangle = \sum_{\mathbf{p} \sigma} \frac{\hbar^2 p^2}{2m} \langle \Phi_0\vert a_{\mathbf{p} \sigma}^\dagger a_{\mathbf{p} \sigma} \vert \Phi_0\rangle . $$From the possible contractions using Wick's theorem, it is straightforward to convince oneself that the expression for the kinetic energy becomes
$$ \langle \Phi_0\vert \hat{T} \vert \Phi_0\rangle = \sum_{\mathbf{i} \leq F} \frac{\hbar^2 k_i^2}{m} = \frac{V}{(2\pi)^3} \frac{\hbar^2}{m} \int_0^{k_F} k^2 d\mathbf{k}. $$The sum of the spin degrees of freedom results in a factor of two only if we deal with identical spin \( 1/2 \) fermions. Changing to spherical coordinates, the integral over the momenta \( k \) results in the final expression
$$ \langle \Phi_0\vert \hat{T} \vert \Phi_0\rangle = \frac{V}{(2\pi)^3} \left( 4\pi \int_0^{k_F} k^4 d\mathbf{k} \right) = \frac{4\pi V}{(2\pi)^3} \frac{1}{5} k_F^5 = \frac{4\pi V}{5(2\pi)^3} k_F^5 = \frac{\hbar^2 V}{10\pi^2 m} k_F^5 . $$The density of states in momentum space is given by \( 2V/(2\pi)^3 \), where we have included the degeneracy due to the spin degrees of freedom. The volume is given by \( 4\pi k_F^3/3 \), and the number of particles becomes
$$ N = \frac{2V}{(2\pi)^3} \frac{4}{3} \pi k_F^3 = \frac{V}{3\pi^2} k_F^3 \quad \Rightarrow \quad k_F = \left( \frac{3\pi^2 N}{V} \right)^{1/3}. $$This gives us
$$ \begin{equation} \langle \Phi_0\vert \hat{T} \vert \Phi_0\rangle = \frac{\hbar^2 V}{10\pi^2 m} \left( \frac{3\pi^2 N}{V} \right)^{5/3} = \frac{\hbar^2 (3\pi^2)^{5/3} N}{10\pi^2 m} \rho^{2/3} , \label{eq:T_forventning} \end{equation} $$We are now ready to calculate the expectation value of the potential energy
$$ \begin{align*} \langle \Phi_0\vert \hat{V} \vert \Phi_0\rangle &= \langle \Phi_0\vert \left( \frac{1}{2V} \sum_{\sigma \sigma'} \sum_{\mathbf{k} \mathbf{p} \mathbf{q} } \left[ \int v(r) \exp{-i (\mathbf{q} \cdot \mathbf{r})} d\mathbf{r} \right] a_{\mathbf{k} + \mathbf{q}, \sigma}^\dagger a_{\mathbf{p} - \mathbf{q}, \sigma'}^\dagger a_{\mathbf{p} \sigma'} a_{\mathbf{k} \sigma} \right) \vert \Phi_0\rangle \\ &= \frac{1}{2V} \sum_{\sigma \sigma'} \sum_{\mathbf{k} \mathbf{p} \mathbf{q}} \left[ \int v(r) \exp{-i (\mathbf{q} \cdot \mathbf{r})} d\mathbf{r} \right]\langle \Phi_0\vert a_{\mathbf{k} + \mathbf{q}, \sigma}^\dagger a_{\mathbf{p} - \mathbf{q}, \sigma'}^\dagger a_{\mathbf{p} \sigma'} a_{\mathbf{k} \sigma} \vert \Phi_0\rangle . \end{align*} $$The only contractions which result in non-zero results are those that involve states below the Fermi level, that is \( k \leq k_F \), \( p \leq k_F \), \( |\mathbf{p} - \mathbf{q}| < \mathbf{k}_F \) and \( |\mathbf{k} + \mathbf{q}| \leq k_F \). Due to momentum conservation we must also have \( \mathbf{k} + \mathbf{q} = \mathbf{p} \), \( \mathbf{p} - \mathbf{q} = \mathbf{k} \) and \( \sigma = \sigma' \) or \( \mathbf{k} + \mathbf{q} = \mathbf{k} \) and \( \mathbf{p} - \mathbf{q} = \mathbf{p} \). Summarizing, we must have
$$ \mathbf{k} + \mathbf{q} = \mathbf{p} \quad \text{and} \quad \sigma = \sigma', \qquad \text{or} \qquad \mathbf{q} = \mathbf{0} . $$We obtain then
$$ \langle \Phi_0\vert \hat{V} \vert \Phi_0\rangle = \frac{1}{2V} \left( \sum_{\sigma \sigma'} \sum_{\mathbf{q} \mathbf{p} \leq F} \left[ \int v(r) d\mathbf{r} \right] - \sum_{\sigma} \sum_{\mathbf{q} \mathbf{p} \leq F} \left[ \int v(r) \exp{-i (\mathbf{q} \cdot \mathbf{r})} d\mathbf{r} \right] \right). $$The first term is the so-called direct term while the second term is the exchange term.
We can rewrite this equation as (and this applies to any potential which depends only on the relative distance between particles)
$$ \begin{equation} \langle \Phi_0\vert \hat{V} \vert \Phi_0\rangle = \frac{1}{2V} \left( N^2 \left[ \int v(r) d\mathbf{r} \right] - N \sum_{\mathbf{q}} \left[ \int v(r) \exp{-i (\mathbf{q}\cdot \mathbf{r})} d\mathbf{r} \right] \right), \label{eq:V_b} \end{equation} $$where we have used the fact that a sum like \( \sum_{\sigma}\sum_{\mathbf{k}} \) equals the number of particles. Using the fact that the density is given by \( \rho = N/V \), with \( V \) being our volume, we can rewrite the last equation as
$$ \langle \Phi_0\vert \hat{V} \vert \Phi_0\rangle = \frac{1}{2} \left( \rho N \left[ \int v(r) d\mathbf{r} \right] - \rho\sum_{\mathbf{q}} \left[ \int v(r) \exp{-i (\mathbf{q}\cdot \mathbf{r})} d\mathbf{r} \right] \right). $$For the electron gas the interaction part of the Hamiltonian operator is given by
$$ \hat{H}_I=\hat{H}_{el}+\hat{H}_{b}+\hat{H}_{el-b}, $$with the electronic part
$$ \hat{H}_{el}=\sum_{i=1}^N\frac{p_i^2}{2m}+\frac{e^2}{2}\sum_{i\ne j}\frac{e^{-\mu |\mathbf{r}_i-\mathbf{r}_j|}}{|\mathbf{r}_i-\mathbf{r}_j|}, $$where we have introduced an explicit convergence factor (the limit \( \mu\rightarrow 0 \) is performed after having calculated the various integrals).
Correspondingly, we have
$$ \hat{H}_{b}=\frac{e^2}{2}\int\int d\mathbf{r}d\mathbf{r}'\frac{n(\mathbf{r})n(\mathbf{r}')e^{-\mu |\mathbf{r}-\mathbf{r}'|}}{|\mathbf{r}-\mathbf{r}'|}, $$which is the energy contribution from the positive background charge with density \( n(\mathbf{r})=N/V \). Finally,
$$ \hat{H}_{el-b}=-\frac{e^2}{2}\sum_{i=1}^N\int d\mathbf{r}\frac{n(\mathbf{r})e^{-\mu |\mathbf{r}-\mathbf{x}_i|}}{|\mathbf{r}-\mathbf{x}_i|}, $$is the interaction between the electrons and the positive background.
We can show that
$$ \hat{H}_{b}=\frac{e^2}{2}\frac{N^2}{V}\frac{4\pi}{\mu^2}, $$and
$$ \hat{H}_{el-b}=-e^2\frac{N^2}{V}\frac{4\pi}{\mu^2}. $$For the electron gas and a Coulomb interaction, these two terms are cancelled (in the thermodynamic limit) by the contribution from the direct term arising from the repulsive electron-electron interaction. What remains then when computing the reference energy is only the kinetic energy contribution and the contribution from the exchange term. For other interactions, like nuclear forces with a short range part and no infinite range, we need to compute both the direct term and the exchange term.
We can show that the final Hamiltonian can be written as
$$ H=H_{0}+H_{I}, $$with
$$ H_{0}={\displaystyle\sum_{\mathbf{k}\sigma}} \frac{\hbar^{2}k^{2}}{2m}a_{\mathbf{k}\sigma}^{\dagger} a_{\mathbf{k}\sigma}, $$and
$$ H_{I}=\frac{e^{2}}{2V}{\displaystyle\sum_{\sigma_{1}\sigma_{2}}}{\displaystyle\sum_{\mathbf{q}\neq 0,\mathbf{k},\mathbf{p}}}\frac{4\pi}{q^{2}} a_{\mathbf{k}+\mathbf{q},\sigma_{1}}^{\dagger} a_{\mathbf{p}-\mathbf{q},\sigma_{2}}^{\dagger} a_{\mathbf{p}\sigma_{2}}a_{\mathbf{k}\sigma_{1}}. $$Calculate \( E_0/N=\langle \Phi_{0}\vert H\vert \Phi_{0}\rangle/N \) for for this system to first order in the interaction. Show that, by using
$$ \rho= \frac{k_F^3}{3\pi^2}=\frac{3}{4\pi r_0^3}, $$with \( \rho=N/V \), \( r_0 \) being the radius of a sphere representing the volume an electron occupies and the Bohr radius \( a_0=\hbar^2/e^2m \), that the energy per electron can be written as
$$ E_0/N=\frac{e^2}{2a_0}\left[\frac{2.21}{r_s^2}-\frac{0.916}{r_s}\right]. $$Here we have defined \( r_s=r_0/a_0 \) to be a dimensionless quantity.
Plot the energy as function of \( r_s \). Why is this system stable? Calculate thermodynamical quantities like the pressure, given by
$$ P=-\left(\frac{\partial E}{\partial V}\right)_N, $$and the bulk modulus
$$ B=-V\left(\frac{\partial P}{\partial V}\right)_N, $$and comment your results.
Hohenberg and Kohn proved that the total energy of a system including that of the many-body effects of electrons (exchange and correlation) in the presence of static external potential (for example, the atomic nuclei) is a unique functional of the charge density. The minimum value of the total energy functional is the ground state energy of the system. The electronic charge density which yields this minimum is then the exact single particle ground state energy.
The electronic energy \( E \) is said to be a functional of the electronic density, \( E[\rho] \), in the sense that for a given function \( \rho(r) \), there is a single corresponding energy. The Hohenberg-Kohn theorems (two) confirms that such a functional exists, but does not tell us the form of the functional.
Density functional theory has turned out to be a resounding success. A large majority of crystalline materials, many molecules and molecular structures have been explained using DFT. Walter Kohn, involved in both the Hohenberg-Kohn theory and development of the Kohn-Sham equations, received the Nobel Prize in Chemistry for DFT and computational chemistry in 1998. DFT is used in a broad spectrum of disciplines, including quantum chemistry, materials science, condensed matter physics and low-energy nuclear physics
Any material on earth, whether in crystals, amorphous solids, molecules or yourself, consists of nothing else than a bunch of atoms, ions and electrons bound together by electric forces. All these possible forms of matter can be explained by virtue of one simple equation: the many-particle Schr\"{o}dinger equation,
$$ \begin{equation} i \hbar \frac{\partial}{\partial t} \Phi( \boldsymbol{r} ; t) = \left( - \sum_{i}^N \frac{\hbar^2}{2m_i} \frac{\partial^2}{\partial \boldsymbol{r}_i^2} + \sum_{i < j}^N \frac{e^2 Z_i Z_j }{ | \boldsymbol{r}_i - \boldsymbol{r}_j| } \right) \Phi( \boldsymbol{r} ; t). \label{_auto2} \end{equation} $$Here \( \Phi(\boldsymbol{r} ; t) \) is the many-body wavefunction for \( N \) particles, where each particle has its own mass \( m_i \), charge \( Z_i \) and position \( \boldsymbol{r}_i \). The only interaction is the Coulomb interaction \( e^2/r \).
In atomic and molecular physics, materials science and quantum chemsitry, the Born-Oppenheimer approximation arises from the physical problem we want to study: the ground state of a collection of interacting ions and electrons. Because even the lightest ion is more than a thousand times heavier than an electron, we approximate our Hamiltonian to not include the dynamics of the ions all-together. This is known as the Born-Oppenheimer approximation.
We then write the time-independent Schr\"{o}dinger equation for a collection of \( N \) electrons subject to the electric potential created by the fixed ions, labeled \( v_{\mathrm{ext}}(\boldsymbol{r}) \),
$$ \left( \sum_{i}^N \left( - \frac{\hbar^2}{2m} \frac{\partial^2}{\partial \boldsymbol{r}_i^2} + v_{\mathrm{ext}}(\boldsymbol{r}_i) \right) + \sum_{i < j}^N \frac{e^2}{|\boldsymbol{r}_i - \boldsymbol{r}_j|} \right) \Psi(\boldsymbol{r}) = E_0 \Psi( \boldsymbol{r} ), $$where \( \boldsymbol{r}_i \) are the positions of the electrons.
We can replace the explicit Coulomb repulsion with a generic two-body interaction
$$ \sum_{i}^N \left( - \frac{\hbar^2}{2m} \frac{\partial^2}{\partial \boldsymbol{r}_i^2} + v_{\mathrm{ext}}(\boldsymbol{r}_i) \right) + \sum_{i < j}^N v(r_{ij})( = E_0 \Psi( \boldsymbol{r} ), $$where \( r_{ij}=\vert \boldsymbol{r}_i-\boldsymbol{r}_j\vert \).
The potential \( v_{\mathrm{ext}}(\boldsymbol{r}_i) \) is created by the charged ions,
$$ v_{\mathrm{ext}}(\boldsymbol{r}_i) = - \sum_j \frac{e^2 Z_j}{| \boldsymbol{r}_i - \boldsymbol{R}_j|} $$where \( \boldsymbol{R} \) is the (static) positions of the ions and \( Z_j \) their charge.
The one-body reduced density matrix of a many-body system at zero temperature gives direct access to many observables, such as the charge density, kinetic energy and occupation numbers. In a coordinate representation it is defined as as the expectation value of the number operator
$$ \hat{n}(\boldsymbol{r})=\sum_{i=1}^N\delta(\boldsymbol{r}-\boldsymbol{r}_i), $$resulting in
$$ n(\boldsymbol{r})=\frac{\langle \Psi \vert \hat{n}(\boldsymbol{r}) \vert \Psi \rangle}{\langle \Psi\vert \Psi \rangle}. $$Here \( \Psi \) can be any state. In the equations that follow we assume it has been normalized properly.
The one-body or just electronic density here, is thus obtained by integrating out all electron degrees of freedom except one, that is (setting \( \boldsymbol{r}=\boldsymbol{r}_1 \))
$$ n(\boldsymbol{r})=n(\boldsymbol{r}_1) = \int d^3 \boldsymbol{r}_2 \cdots d^3 \boldsymbol{r}_N \left| \Psi (\boldsymbol{r}_1 \cdots \boldsymbol{r}_n) \right|^2. $$If \( \Psi \) is a single reference Slater determinant defined in terms of the single-particle functions \( \psi_i \), then we have
$$ n(\boldsymbol{r})=\sum_{i=1}^N\vert \psi_i(\boldsymbol{r})\vert^2, $$as derived earlier.
Text to come here.
Assume we found a solution for the Hamiltonian from the Born-Oppenheimer approximation, with ground state energy \( E_0 \) and a certain electronic density \( n(\boldsymbol{r}) \). The strength of the Coulomb interaction and the mass of an electron are constants of nature, so the only input that can possibly influence the electronic density \( n(\boldsymbol{r}) \) and the energy \( E_0 \) of our ground state is our choice of potential \( v_{\mathrm{ext}}(\boldsymbol{r}) \). In other words, the ground state energy is a {\em functional} of the input potential,
$$ E_0 [v_{\mathrm{ext}}(\boldsymbol{r})] = \mathcal{F}_E[ v_{\mathrm{ext}}(\boldsymbol{r}) ] $$A functional is nothing else than a function whose input is another function; in this case the functional \( \mathcal{F} \) takes as input the electric potential generated by the ions and outputs the ground state energy based on thr Born-Oppenheimer approximation.
At first this results seems counterintuitive. After all, the ground state energy clearly contains the kinetic energy \( T \), the interaction energy \( U \) and the potential energy. Only the latter term explicitly depends on the potential. We can thus write the ground state energy in terms of a separate functional for the kinetic and interaction energy, and the potential energy
Knowing this functional, for any given potential \( V(\boldsymbol{r}) \) we minimize the right hand side by checking all possible electronic density distributions.
There are only two minor problems. We don't know what this functional looks like. And even if we did, we don't know how to find the right electronic density.
The unknown functional \( \mathcal{F}[n(\boldsymbol{r})] \) should describe the kinetic and interaction energy of a system described by for example the Born-Oppenheimer approximation. Even though we cannot find its exact shape, we can look at its shape in some limiting cases that we can solve.
We know that a free homogeneous electron gas (HEG) with density \( n \) has a ground state energy of
$$ E_0 = \frac{3 \hbar^2 \left( 3 \pi^2 \right)^{2/3}}{10 m} n_0^{5/3}. $$For a slowly varying electronic density, we can approximate the kinetic energy contribution to the full functional \( \mathcal{F}[n(\boldsymbol{r})] \) as the energy evaluated at each point separately,
$$ \mathcal{T}_0 [n(\boldsymbol{r})] = \frac{3 \hbar^2 \left( 3 \pi^2 \right)^{2/3}}{10 m} \int d^3\boldsymbol{r} ( n(\boldsymbol{r}) )^{5/3}. $$We have also earlier derived that the simplest energy contribution from Coulomb interactions is given by the Hartree term
$$ \mathcal{U}_{\mathrm{H}} [n(\boldsymbol{r})] = \frac{e^2}{2} \int d^3\boldsymbol{r} d^3\boldsymbol{r}' \, \frac{n(\boldsymbol{r}) n(\boldsymbol{r}')}{|\boldsymbol{r} - \boldsymbol{r}'|}. $$It is natural to write out the full functional as containing the homogeneous electron gas term and the Hartree term. The remaining terms, though still unknown, should be small. This unknown part is conventionally called the {\em exchange-correlation potential} \( E_{xc}[n(\boldsymbol{r})] \). The full Hohenberg-Kohn functional, including the potential energy, is thus
$$ \mathcal{E}_{\mathrm{HK}}[n(\boldsymbol{r})] = \mathcal{T}_0 [n(\boldsymbol{r})] + \int d^3 \boldsymbol{r} v_{\mathrm{ext}}(r) n(\boldsymbol{r}) + \mathcal{U}_{\mathrm{H}} [n(\boldsymbol{r})] + E_{xc} [ n(\boldsymbol{r})]. $$We replaced an intractable problem with the task of minimizing an unknown functional \( \mathcal{F}[n(\boldsymbol{r})] \) over infinitely many possible electronic densities \( n(\boldsymbol{r}) \). In the previous section we already gave some first suggestions for the functional. But once we found it, how to find the right electronic density \( n(\boldsymbol{r}) \)?
Because the correct density minimizes the functional, we can find the functional by setting it's derivative to zero,
$$ \frac{\delta \mathcal{F}[n(\boldsymbol{r})] }{\delta n(\boldsymbol{r})} = 0. $$Using the functional Eq(insert), we write out
$$ \frac{\delta \mathcal{T}[n(\boldsymbol{r})] }{\delta n(\boldsymbol{r})} + v_{\mathrm{ext}}(\boldsymbol{r}) + \int \frac{n(\boldsymbol{r}')}{|\boldsymbol{r}-\boldsymbol{r}'|}d^3 \boldsymbol{r}' + \frac{\delta E_{xc} [n(\boldsymbol{r})] }{\delta n(\boldsymbol{r})} =0. $$The idea of Kohn and Sham was to treat this as if it is a single-particle problem. The first term represents the kinetic energy, and the remaining terms form the Kohn-Sham potential
$$ V_{\mathrm{KS}}(\boldsymbol{r}) = V(\boldsymbol{r}) + \int \frac{n(\boldsymbol{r}')}{|\boldsymbol{r}-\boldsymbol{r}'|}d^3 \boldsymbol{r}' + \frac{\delta E_{xc} [n(\boldsymbol{r})] }{\delta n(\boldsymbol{r})}. $$The Kohn-Sham equation is the single-particle Schr\"{o}dinger equation with the potential given by the last equation,
$$ \left( - \frac{\hbar^2}{2m} \frac{\partial^2}{\partial \boldsymbol{r}^2} + V_{\mathrm{KS}}(\boldsymbol{r}) \right) \psi_i(\boldsymbol{r}) = \epsilon_i \psi_i (\boldsymbol{r}). $$The Kohn-Sham equations have to be solved numerically. As with our discussions of Hartree-Fock theory, this is tractable because the equations are just given by a set of coupled differential equation. The electronic density is obtained by occupying the \( N \) solutions \( \psi_i(\boldsymbol{r}) \) with the lowest energy,
$$ n(\boldsymbol{r}) = \sum_{i=1}^N |\psi_i (\boldsymbol{r}) |^2. $$Now the electronic density obtained this way can be used to calculate a new Kohn-Sham potential. We continue this iterative procedure until we reach convergence.