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.
In Hartree-Fock theory one works with large basis sets. This poses a problem for large systems. An alternative to the HF methods is DFT. DFT takes into account electron correlations but is less demanding computationally than full scale diagonalization, Coupled Cluster theory or say Monte Carlo methods.
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{_auto1} \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.
As defined above, the one-body density in coordinate space is defined as
$$ \rho({\bf r}) = \sum_{i=1}^N \delta({\bf r}-{\bf r}_i). $$In second quantization this becomes
$$ \hat{\rho}({\bf r}) = \sum_{\alpha\beta}^{\infty}\rho_{\alpha,\beta}({\bf r}) a^{\dagger}_{\alpha}a_{\beta}. $$where
$$ \rho_{\alpha,\beta}({\bf r})= \psi^*_{\alpha}({\bf r})\psi_{\beta}({\bf r}). $$The number of particles is \( N \) and the integral of the expectation value of the one-body density operator should give you \( N \) particles. With an appropriate similarity transformation we can make this operator diagonal in the single-particle basis \( \psi_{\alpha} \). That is
$$ \hat{\rho}({\bf r}) = \sum_{\alpha=1}^{\infty} |\psi_{\alpha}({\bf r})|^2a^{\dagger}_{\alpha}a_{\alpha}= \hat{\rho}({\bf r}) = \sum_{\alpha=1}^{\infty}\rho_{\alpha\alpha}({\bf r}) a^{\dagger}_{\alpha}a_{\alpha}. $$The ground state wave function \( |\Psi_0\rangle \) is a linear combination of all \( D \) possible Slater determinants \( |\Phi_i\rangle \)
$$ |\Psi_0\rangle=\sum_{i=1}^DC_{0i}|\Phi_i\rangle, $$where the coefficients \( C_{0i} \) could arise from an FCI calculation using a given Slater determinant basis
$$ |\Phi_i\rangle = a^{\dagger}_1a^{\dagger}_2\dots a^{\dagger}_N |0\rangle. $$The ground state expectation value of the one-body density operator is
$$ \langle\hat{\rho}({\bf r})\rangle = \langle \Psi_0|\sum_{\alpha=1}^{\infty} \rho_{\alpha,\alpha}({\bf r})a^{\dagger}_{\alpha}a_{\alpha}|\Psi_0\rangle, $$which translates into
$$ \langle\hat{\rho}({\bf r})\rangle = \sum_{ij=1}^DC^*_{0i}C_{0j}\langle \Phi_i|\sum_{\alpha=1}^{\infty} \rho_{\alpha,\alpha}({\bf r})a^{\dagger}_{\alpha}a_{\alpha}|\Phi_j\rangle. $$Integrating
$$ \int \langle\hat{\rho}({\bf r})\rangle d {\bf r}, $$gives us \( N \), the number of particles!
The two-body densities is a simple extension of the one-body density
$$ \rho({\bf r}_1, {\bf r}_2) = \sum_{ij=1}^N \delta({\bf r}_1-{\bf r}_i)\delta({\bf r}_2-{\bf r}_j), $$which in second quantization becomes
$$ \hat{\rho}({\bf r}_1, {\bf r}_2) = \sum_{\alpha\beta\gamma\delta}^{\infty}\rho_{\alpha,\gamma}({\bf r}_1)\rho_{\beta,\delta}({\bf r}_2)a^{\dagger}_{\alpha}a^{\dagger}_{\beta}a_{\delta}a_{\gamma}, $$meaning that the ground-state two-body density is
$$ \langle\hat{\rho}({\bf r}_1, {\bf r}_2) \rangle = \sum_{ij=1}^DC^*_{0i}C_{0j}\langle \Phi_i|\sum_{\alpha\beta\gamma\delta}^{\infty}\rho_{\alpha,\gamma}({\bf r}_1)\rho_{\beta,\delta}({\bf r}_2)a^{\dagger}_{\alpha}a^{\dagger}_{\beta}a_{\delta}a_{\gamma} |\Phi_j\rangle. $$The Hartree-Fock algorithm can be broken down as follows. We recall that our Hartree-Fock matrix is
$$ \hat{h}_{\alpha\beta}^{HF}=\langle \alpha \vert\hat{h}_0 \vert \beta \rangle+ \sum_{j=1}^N\sum_{\gamma\delta} C^*_{j\gamma}C_{j\delta}\langle \alpha\gamma|V|\beta\delta\rangle_{AS}. $$Normally we assume that the single-particle basis \( \vert\beta\rangle \) forms an eigenbasis for the operator \( \hat{h}_0 \), meaning that the Hartree-Fock matrix becomes
$$ \hat{h}_{\alpha\beta}^{HF}=\epsilon_{\alpha}\delta_{\alpha,\beta}+ \sum_{j=1}^N\sum_{\gamma\delta} C^*_{j\gamma}C_{j\delta}\langle \alpha\gamma|V|\beta\delta\rangle_{AS}. $$The Hartree-Fock eigenvalue problem
$$ \sum_{\beta}\hat{h}_{\alpha\beta}^{HF}C_{i\beta}=\epsilon_i^{\mathrm{HF}}C_{i\alpha}, $$can be written out in a more compact form as
$$ \hat{h}^{HF}\hat{C}=\epsilon^{\mathrm{HF}}\hat{C}. $$The equations are often rewritten in terms of a so-called density matrix, which is defined as
$$ \rho_{\gamma\delta}=\sum_{i=1}^{N}\langle\gamma|i\rangle\langle i|\delta\rangle = \sum_{i=1}^{N}C_{i\gamma}C^*_{i\delta}. $$It means that we can rewrite the Hartree-Fock Hamiltonian as
$$ \hat{h}_{\alpha\beta}^{HF}=\epsilon_{\alpha}\delta_{\alpha,\beta}+ \sum_{\gamma\delta} \rho_{\gamma\delta}\langle \alpha\gamma|V|\beta\delta\rangle_{AS}. $$It is convenient to use the density matrix since we can precalculate in every iteration the product of two eigenvector components \( C \).
Note that \( \langle \alpha\vert\hat{h}_0\vert\beta \rangle \) denotes the matrix elements of the one-body part of the starting hamiltonian.
The electronic energy \( E \) is said to be a \emph{functional} of the electronic density, \( E[\rho] \), in the sense that for a given function \( \rho(r) \), there is a single corresponding energy. The \emph{Hohenberg-Kohn theorem} confirms that such a functional exists, but does not tell us the form of the functional. As shown by Kohn and Sham, the exact ground-state energy \( E \) of an \( N \)-electron system can be written as
$$ E[\rho] = -\frac{1}{2} \sum_{i=1}^N\int \Psi_i^*(\mathbf{r_1})\nabla_1^2 \Psi_i(\mathbf{r_1}) d\mathbf{r_1} - \int \frac{Z}{r_1} \rho(\mathbf{r_1}) d\mathbf{r_1} + \frac{1}{2} \int\frac{\rho(\mathbf{r_1})\rho(\mathbf{r_2})}{r_{12}} d\mathbf{r_1}d\mathbf{r_2} + E_{XC}[\rho] $$with \( \Psi_i \) the \emph{Kohn-Sham} (KS) \emph{orbitals}. Note that we have limited ourselves to atomic physics here.
How do we arrive at the above equation?
The ground-state charge density is given by
$$ \rho(\mathbf{r}) = \sum_{i=1}^N|\Psi_i(\mathbf{r})|^2, $$where the sum is over the occupied Kohn-Sham orbitals. The last term, \( E_{XC}[\rho] \), is the \emph{exchange-correlation energy} which in theory takes into account all non-classical electron-electron interaction. However, we do not know how to obtain this term exactly, and are forced to approximate it. The KS orbitals are found by solving the \emph{Kohn-Sham equations}, which can be found by applying a variational principle to the electronic energy \( E[\rho] \). This approach is similar to the one used for obtaining the HF equation.
The KS equations reads
$$ \left\{ -\frac{1}{2}\nabla_1^2 - \frac{Z}{r_1} + \int \frac{\rho(\mathbf{r_2})}{r_{12}} d\mathbf{r_2} + V_{XC}(\mathbf{r_1}) \right\} \Psi_i(\mathbf{r_1}) = \epsilon_i \Psi_i(\mathbf{r_1}) $$where \( \epsilon_i \) are the KS orbital energies, and where the \emph{exchange-correlation potential} is given by
$$ V_{XC}[\rho] = \frac{\delta E_{XC}[\rho]}{\delta \rho}. $$The KS equations are solved in a self-consistent fashion. A variety of basis set functions can be used, and the experience gained in HF calculations are often useful. The computational time needed for a DFT calculation formally scales as the third power of the number of basis functions.
The main source of error in DFT usually arises from the approximate nature of \( E_{XC} \). In the \emph{local density approximation} (LDA) it is approximated as
$$ E_{XC} = \int \rho(\mathbf{r})\epsilon_{XC}[\rho(\mathbf{r})] d\mathbf{r}, $$where \( \epsilon_{XC}[\rho(\mathbf{r})] \) is the exchange-correlation energy per electron in a homogeneous electron gas of constant density. The LDA approach is clearly an approximation as the charge is not continuously distributed. To account for the inhomogeneity of the electron density, a nonlocal correction involving the gradient of \( \rho \) is often added to the exchange-correlation energy.
We assume that there is a \( \mathcal{V}_{\mathrm{ext}} = \) set of external single-particle \alert{potentials} \( v \) so that
$$ \begin{equation*} \hat{H}\vert \phi\rangle = \left(\hat{T}+\hat{V}_{\mathrm{ext}}+\hat{V}\right)=E\vert\phi\rangle,\qquad \hat{V}_{\mathrm{ext}}\in \mathcal{V}_{\mathrm{ext}},\nonumber \end{equation*} $$gives a \alert{non-degenerate} N-particle ground state \( \vert\Psi \rangle_0 \). For any system of interacting particles in an external potential \( \mathcal{V}_{\mathrm{ext}} \), the potential \( \mathcal{V}_{\mathrm{ext}} \) is uniquely determined (by a near constant) by the ground state density \( \rho_0 \). There is a corollary to this statement which states that since \( \hat{H} \) is determined, the many-body functions for all states are also determined. All properties of the system are determined via \( \rho_0 \).
The density (assuming normalized state vectors)
$$ \begin{equation*} \rho(\mathbf{r})=\sum_{i}\int dx_{2}\dots \int dx_{N}\vert \Psi(\mathbf{r},x_{2},\dots ,x_{N})\vert^{2} \nonumber \end{equation*} $$Theorem II states that a universal functional for the energy \( E[\rho] \) (function of \( \rho \)) can be defined for every external potential \( \mathcal{U}_{\mathrm{ext}} \). For a given external potential, the exact ground state energy of the system is a global minimum of this functional. The density which minimizes this functional is \( \rho_0 \).
Proof I Let us prove \( C:\mathcal{V}(C)\longrightarrow \Psi \) injective:
$$ \begin{equation*} \hat{V}\neq \hat{V}'+\text{constant} \qquad \stackrel{{\LARGE ?}}{\Longrightarrow } \qquad \vert \Psi\rangle\neq \vert \Psi'\rangle, \nonumber \end{equation*} $$where \( \hat{V}, \hat{V}' \in \mathcal{V} \)
\textbf{\emph{Reductio ad absurdum}}: Assume \( \vert \Psi\rangle=\vert \Psi'\rangle \) for some \( \hat{V}\neq \hat{V}'+\text{const} \), \( \hat{V}, \hat{V}' \in \mathcal{V} \)
\( \hat{T}\neq \hat{T}[V] \), $\hat{W}\neq \hat{W}[V] \quad \Longrightarrow $\footnote{Unique continuation theorem: \( \vert \Psi\rangle\neq 0 \) on a set of positive measure}
$$ \begin{equation*} \left(\hat{V}-\hat{V}'\right)\vert \Psi \rangle=\left(E_{gs}-E_{gs}'\right)\vert \Psi\rangle.\nonumber \end{equation*} $$ $$ \begin{align} \Longrightarrow \qquad & \hat{V}-\hat{V}'=E_{gs}-E_{gs}' \nonumber \\ \Longrightarrow \qquad & \hat{V}=\hat{V}'+\text{constant} \qquad \alert{\text{Contradiction!}}\nonumber \end{align} $$Proof II: Let us prove \( D:\Psi \longrightarrow \mathcal{N} \) injective:
$$ \begin{equation*} \vert \Psi\rangle\neq \vert \Psi'\rangle \qquad \stackrel{{\LARGE ?}}{\Longrightarrow } \qquad \rho(\mathbf{r})\neq n'(\mathbf{r}) \nonumber \end{equation*} $$\textbf{\emph{Reductio ad absurdum}}: Assume \( \rho(\mathbf{r})=n'(\mathbf{r}) \) for some \( \vert \Psi\rangle\neq \vert \Psi'\rangle \)
Ritz principle \( \quad \Longrightarrow \)
$$ \begin{equation} E_{gs}=\langle \Psi\vert \hat{H}\vert \Psi\rangle \langle\Psi'\vert \hat{H}\vert \Psi'\rangle. \label{_auto2} \end{equation} $$ $$ \begin{equation*} \langle\Psi'\vert \hat{H}\vert \Psi'\rangle\langle\Psi'\vert\hat{H}'+\hat{V}-\hat{V}'\vert \Psi'\rangle=E_{gs}'+\int n'(\mathbf{r})[v(\mathbf{r})-v'(\mathbf{r})]d^{3}r \nonumber \end{equation*} $$ $$ \begin{equation*} \Longrightarrow \qquad E_{gs}' < E_{gs}+\int n'(\mathbf{r})[v(\mathbf{r})-v'(\mathbf{r})]d^{3}r %\nonumber \end{equation*} $$By symmetry
$$ \begin{equation*} \Longrightarrow \qquad E_{gs}< E_{gs}'+\int n'(\mathbf{r})[v'(\mathbf{r})-v(\mathbf{r})]d^{3}r %\nonumber \end{equation*} $$ $$ \begin{equation*} E_{gs}+E_{gs}' < E_{gs}+E_{gs}' \qquad \text{\alert{Contradiction!}} \nonumber \end{equation*} $$Define
$$ \begin{equation} E_{v_{0}}[\rho]:=\langle \Psi[\rho]\vert \hat{T}+\hat{W}+\hat{V_{0}}\vert \Psi[\rho]\rangle \label{_auto3} \end{equation} $$\( \hat{V_{0}} = \) external potential, \( n_{0}(\mathbf{r}) = \) corresponding GS density, \( E_{0} = \) GS energy
Rayleigh-Ritz principle \( \quad \Longrightarrow \quad \) \alert{second statement of H-K theorem}:
$$ \begin{equation*} E_{0}=\min_{n\in \mathcal{N}}E_{v_{0}}[\rho] \nonumber \end{equation*} $$Last satement of H-K theorem:
$$ \begin{equation} F_{HK}[\rho]\equiv \langle\Psi[\rho]\vert\hat{T}+\hat{W}\vert \Psi[\rho]\rangle \label{_auto4} \end{equation} $$is \emph{universal} (\( F_{HK}\neq F_{HK}[\hat{V_{0}}] \))
The classic Kohn-Sham scheme:
$$ \left(-\frac{\hbar^{2}}{2m}\nabla^{2}+v_{s,0}(\mathbf{r})\right)\phi_{i,0}(\mathbf{r})=\varepsilon_{i}\phi_{i,0}(\mathbf{r}), \qquad \varepsilon_{1}\geq \varepsilon_{2} \geq \dots, $$where
$$ v_{s,0}(\mathbf{r})=v_{0}(\mathbf{r})+\int d^{3}r' w(\mathbf{r},\mathbf{r}')\rho_{0}(\mathbf{r}')+v_{\mathrm{XC}}([\rho_{0}];\mathbf{r}) $$The density is calculated as
$$ \rho_{0}(\mathbf{r})=\sum_{i=1}^{N}\vert \phi_{i,0}(\mathbf{r})\vert^{2}. $$The equation is solved selfconsistently. The total energy
$$ E=\sum_{i=1}^{N}\varepsilon_{i}-\frac{1}{2}\int d^{3}r d^{3}r'\rho(\mathbf{r})w(\mathbf{r},\mathbf{r}')\rho(\mathbf{r}')+E_{\mathrm{XC}}[\rho]-\int d^{3}r v_{\mathrm{XC}}([\rho];\mathbf{r})\rho(\mathbf{r}) $$We assume here that we are only interested in the ground state of the system and expand the exact wave function in term of a series of Slater determinants
$$ \vert \Psi_0\rangle = \vert \Phi_0\rangle + \sum_{m=1}^{\infty}C_m\vert \Phi_m\rangle, $$where we have assumed that the true ground state is dominated by the solution of the unperturbed problem, that is
$$ \hat{H}_0\vert \Phi_0\rangle= W_0\vert \Phi_0\rangle. $$The state \( \vert \Psi_0\rangle \) is not normalized, rather we have used an intermediate normalization \( \langle \Phi_0 \vert \Psi_0\rangle=1 \) since we have \( \langle \Phi_0\vert \Phi_0\rangle=1 \).
The Schroedinger equation is
$$ \hat{H}\vert \Psi_0\rangle = E\vert \Psi_0\rangle, $$and multiplying the latter from the left with \( \langle \Phi_0\vert \) gives
$$ \langle \Phi_0\vert \hat{H}\vert \Psi_0\rangle = E\langle \Phi_0\vert \Psi_0\rangle=E, $$and subtracting from this equation
$$ \langle \Psi_0\vert \hat{H}_0\vert \Phi_0\rangle= W_0\langle \Psi_0\vert \Phi_0\rangle=W_0, $$and using the fact that the both operators \( \hat{H} \) and \( \hat{H}_0 \) are hermitian results in
$$ \Delta E=E-W_0=\langle \Phi_0\vert \hat{H}_I\vert \Psi_0\rangle, $$which is an exact result. We call this quantity the correlation energy.
This equation forms the starting point for all perturbative derivations. However, as it stands it represents nothing but a mere formal rewriting of Schroedinger's equation and is not of much practical use. The exact wave function \( \vert \Psi_0\rangle \) is unknown. In order to obtain a perturbative expansion, we need to expand the exact wave function in terms of the interaction \( \hat{H}_I \).
Here we have assumed that our model space defined by the operator \( \hat{P} \) is one-dimensional, meaning that
$$ \hat{P}= \vert \Phi_0\rangle \langle \Phi_0\vert , $$and
$$ \hat{Q}=\sum_{m=1}^{\infty}\vert \Phi_m\rangle \langle \Phi_m\vert . $$We can thus rewrite the exact wave function as
$$ \vert \Psi_0\rangle= (\hat{P}+\hat{Q})\vert \Psi_0\rangle=\vert \Phi_0\rangle+\hat{Q}\vert \Psi_0\rangle. $$Going back to the Schr\"odinger equation, we can rewrite it as, adding and a subtracting a term \( \omega \vert \Psi_0\rangle \) as
$$ \left(\omega-\hat{H}_0\right)\vert \Psi_0\rangle=\left(\omega-E+\hat{H}_I\right)\vert \Psi_0\rangle, $$where \( \omega \) is an energy variable to be specified later.
We assume also that the resolvent of \( \left(\omega-\hat{H}_0\right) \) exits, that is it has an inverse which defined the unperturbed Green's function as
$$ \left(\omega-\hat{H}_0\right)^{-1}=\frac{1}{\left(\omega-\hat{H}_0\right)}. $$We can rewrite Schroedinger's equation as
$$ \vert \Psi_0\rangle=\frac{1}{\omega-\hat{H}_0}\left(\omega-E+\hat{H}_I\right)\vert \Psi_0\rangle, $$and multiplying from the left with \( \hat{Q} \) results in
$$ \hat{Q}\vert \Psi_0\rangle=\frac{\hat{Q}}{\omega-\hat{H}_0}\left(\omega-E+\hat{H}_I\right)\vert \Psi_0\rangle, $$which is possible since we have defined the operator \( \hat{Q} \) in terms of the eigenfunctions of \( \hat{H} \).
These operators commute meaning that
$$ \hat{Q}\frac{1}{\left(\omega-\hat{H}_0\right)}\hat{Q}=\hat{Q}\frac{1}{\left(\omega-\hat{H}_0\right)}=\frac{\hat{Q}}{\left(\omega-\hat{H}_0\right)}. $$With these definitions we can in turn define the wave function as
$$ \vert \Psi_0\rangle=\vert \Phi_0\rangle+\frac{\hat{Q}}{\omega-\hat{H}_0}\left(\omega-E+\hat{H}_I\right)\vert \Psi_0\rangle. $$This equation is again nothing but a formal rewrite of Schr\"odinger's equation and does not represent a practical calculational scheme. It is a non-linear equation in two unknown quantities, the energy \( E \) and the exact wave function \( \vert \Psi_0\rangle \). We can however start with a guess for \( \vert \Psi_0\rangle \) on the right hand side of the last equation.
The most common choice is to start with the function which is expected to exhibit the largest overlap with the wave function we are searching after, namely \( \vert \Phi_0\rangle \). This can again be inserted in the solution for \( \vert \Psi_0\rangle \) in an iterative fashion and if we continue along these lines we end up with
$$ \vert \Psi_0\rangle=\sum_{i=0}^{\infty}\left\{\frac{\hat{Q}}{\omega-\hat{H}_0}\left(\omega-E+\hat{H}_I\right)\right\}^i\vert \Phi_0\rangle, $$for the wave function and
$$ \Delta E=\sum_{i=0}^{\infty}\langle \Phi_0\vert \hat{H}_I\left\{\frac{\hat{Q}}{\omega-\hat{H}_0}\left(\omega-E+\hat{H}_I\right)\right\}^i\vert \Phi_0\rangle, $$which is now a perturbative expansion of the exact energy in terms of the interaction \( \hat{H}_I \) and the unperturbed wave function \( \vert \Psi_0\rangle \).
In our equations for \( \vert \Psi_0\rangle \) and \( \Delta E \) in terms of the unperturbed solutions \( \vert \Phi_i\rangle \) we have still an undetermined parameter \( \omega \) and a dependecy on the exact energy \( E \). Not much has been gained thus from a practical computational point of view.
In Brilluoin-Wigner perturbation theory it is customary to set \( \omega=E \). This results in the following perturbative expansion for the energy \( \Delta E \)
$$ \Delta E=\sum_{i=0}^{\infty}\langle \Phi_0\vert \hat{H}_I\left\{\frac{\hat{Q}}{\omega-\hat{H}_0}\left(\omega-E+\hat{H}_I\right)\right\}^i\vert \Phi_0\rangle= $$ $$ \langle \Phi_0\vert \left(\hat{H}_I+\hat{H}_I\frac{\hat{Q}}{E-\hat{H}_0}\hat{H}_I+ \hat{H}_I\frac{\hat{Q}}{E-\hat{H}_0}\hat{H}_I\frac{\hat{Q}}{E-\hat{H}_0}\hat{H}_I+\dots\right)\vert \Phi_0\rangle. $$This expression depends however on the exact energy \( E \) and is again not very convenient from a practical point of view. It can obviously be solved iteratively, by starting with a guess for \( E \) and then solve till some kind of self-consistency criterion has been reached.
Actually, the above expression is nothing but a rewrite again of the full Schr\"odinger equation.
Defining \( e=E-\hat{H}_0 \) and recalling that \( \hat{H}_0 \) commutes with \( \hat{Q} \) by construction and that \( \hat{Q} \) is an idempotent operator \( \hat{Q}^2=\hat{Q} \). Using this equation in the above expansion for \( \Delta E \) we can write the denominator
$$ \hat{Q}\frac{1}{\hat{e}-\hat{Q}\hat{H}_I\hat{Q}}= $$ $$ \hat{Q}\left[\frac{1}{\hat{e}}+\frac{1}{\hat{e}}\hat{Q}\hat{H}_I\hat{Q} \frac{1}{\hat{e}}+\frac{1}{\hat{e}}\hat{Q}\hat{H}_I\hat{Q} \frac{1}{\hat{e}}\hat{Q}\hat{H}_I\hat{Q}\frac{1}{\hat{e}}+\dots\right]\hat{Q}. $$Inserted in the expression for \( \Delta E \) leads to
$$ \Delta E= \langle \Phi_0\vert \hat{H}_I+\hat{H}_I\hat{Q}\frac{1}{E-\hat{H}_0-\hat{Q}\hat{H}_I\hat{Q}}\hat{Q}\hat{H}_I\vert \Phi_0\rangle. $$In RS perturbation theory we set \( \omega = W_0 \) and obtain the following expression for the energy difference
$$ \Delta E=\sum_{i=0}^{\infty}\langle \Phi_0\vert \hat{H}_I\left\{\frac{\hat{Q}}{W_0-\hat{H}_0}\left(\hat{H}_I-\Delta E\right)\right\}^i\vert \Phi_0\rangle= $$ $$ \langle \Phi_0\vert \left(\hat{H}_I+\hat{H}_I\frac{\hat{Q}}{W_0-\hat{H}_0}(\hat{H}_I-\Delta E)+ \hat{H}_I\frac{\hat{Q}}{W_0-\hat{H}_0}(\hat{H}_I-\Delta E)\frac{\hat{Q}}{W_0-\hat{H}_0}(\hat{H}_I-\Delta E)+\dots\right)\vert \Phi_0\rangle. $$Recalling that \( \hat{Q} \) commutes with \( \hat{H_0} \) and since \( \Delta E \) is a constant we obtain that
$$ \hat{Q}\Delta E\vert \Phi_0\rangle = \hat{Q}\Delta E\vert \hat{Q}\Phi_0\rangle = 0. $$Inserting this results in the expression for the energy results in
$$ \Delta E=\langle \Phi_0\vert \left(\hat{H}_I+\hat{H}_I\frac{\hat{Q}}{W_0-\hat{H}_0}\hat{H}_I+ \hat{H}_I\frac{\hat{Q}}{W_0-\hat{H}_0}(\hat{H}_I-\Delta E)\frac{\hat{Q}}{W_0-\hat{H}_0}\hat{H}_I+\dots\right)\vert \Phi_0\rangle. $$We can now this expression in terms of a perturbative expression in terms of \( \hat{H}_I \) where we iterate the last expression in terms of \( \Delta E \)
$$ \Delta E=\sum_{i=1}^{\infty}\Delta E^{(i)}. $$We get the following expression for \( \Delta E^{(i)} \)
$$ \Delta E^{(1)}=\langle \Phi_0\vert \hat{H}_I\vert \Phi_0\rangle, $$which is just the contribution to first order in perturbation theory,
$$ \Delta E^{(2)}=\langle\Phi_0\vert \hat{H}_I\frac{\hat{Q}}{W_0-\hat{H}_0}\hat{H}_I\vert \Phi_0\rangle, $$which is the contribution to second order.
$$ \Delta E^{(3)}=\langle \Phi_0\vert \hat{H}_I\frac{\hat{Q}}{W_0-\hat{H}_0}\hat{H}_I\frac{\hat{Q}}{W_0-\hat{H}_0}\hat{H}_I\Phi_0\rangle- \langle\Phi_0\vert \hat{H}_I\frac{\hat{Q}}{W_0-\hat{H}_0}\langle \Phi_0\vert \hat{H}_I\vert \Phi_0\rangle\frac{\hat{Q}}{W_0-\hat{H}_0}\hat{H}_I\vert \Phi_0\rangle, $$being the third-order contribution.
In the shell-model lectures we showed that we could rewrite the exact state function for say the ground state, as a linear expansion in terms of all possible Slater determinants. That is, we define the ansatz for the ground state as
$$ |\Phi_0\rangle = \left(\prod_{i\le F}\hat{a}_{i}^{\dagger}\right)|0\rangle, $$where the index \( i \) defines different single-particle states up to the Fermi level. We have assumed that we have \( N \) fermions.
A given one-particle-one-hole (\( 1p1h \)) state can be written as
$$ |\Phi_i^a\rangle = \hat{a}_{a}^{\dagger}\hat{a}_i|\Phi_0\rangle, $$while a \( 2p2h \) state can be written as
$$ |\Phi_{ij}^{ab}\rangle = \hat{a}_{a}^{\dagger}\hat{a}_{b}^{\dagger}\hat{a}_j\hat{a}_i|\Phi_0\rangle, $$and a general \( ApAh \) state as
$$ |\Phi_{ijk\dots}^{abc\dots}\rangle = \hat{a}_{a}^{\dagger}\hat{a}_{b}^{\dagger}\hat{a}_{c}^{\dagger}\dots\hat{a}_k\hat{a}_j\hat{a}_i|\Phi_0\rangle. $$We use letters \( ijkl\dots \) for states below the Fermi level and \( abcd\dots \) for states above the Fermi level. A general single-particle state is given by letters \( pqrs\dots \).
We can then expand our exact state function for the ground state as
$$ |\Psi_0\rangle=C_0|\Phi_0\rangle+\sum_{ai}C_i^a|\Phi_i^a\rangle+\sum_{abij}C_{ij}^{ab}|\Phi_{ij}^{ab}\rangle+\dots =(C_0+\hat{C})|\Phi_0\rangle, $$where we have introduced the so-called correlation operator
$$ \hat{C}=\sum_{ai}C_i^a\hat{a}_{a}^{\dagger}\hat{a}_i +\sum_{abij}C_{ij}^{ab}\hat{a}_{a}^{\dagger}\hat{a}_{b}^{\dagger}\hat{a}_j\hat{a}_i+\dots $$Since the normalization of \( \Psi_0 \) is at our disposal and since \( C_0 \) is by hypothesis non-zero, we may arbitrarily set \( C_0=1 \) with corresponding proportional changes in all other coefficients. Using this so-called intermediate normalization we have
$$ \langle \Psi_0 | \Phi_0 \rangle = \langle \Phi_0 | \Phi_0 \rangle = 1, $$resulting in
$$ |\Psi_0\rangle=(1+\hat{C})|\Phi_0\rangle. $$In an FCI calculation, the unknown coefficients in \( \hat{C} \) are the eigenvectors which result from the diagonalization of the Hamiltonian matrix.
How can we use perturbation theory to determine the same coefficients? Let us study the contributions to second order in the interaction, namely
$$ \Delta E^{(2)}=\langle\Phi_0\vert \hat{H}_I\frac{\hat{Q}}{W_0-\hat{H}_0}\hat{H}_I\vert \Phi_0\rangle. $$The intermediate states given by \( \hat{Q} \) can at most be of a \( 2p-2h \) nature if we have a two-body Hamiltonian. This means that second order in the perturbation theory can have \( 1p-1h \) and \( 2p-2h \) at most as intermediate states. When we diagonalize, these contributions are included to infinite order. This means that higher-orders in perturbation theory bring in more complicated correlations.
If we limit the attention to a Hartree-Fock basis, then we have that \( \langle\Phi_0\vert \hat{H}_I \vert 2p-2h\rangle \) is the only contribution and the contribution to the energy reduces to
$$ \Delta E^{(2)}=\frac{1}{4}\sum_{abij}\langle ij\vert \hat{v}\vert ab\rangle \frac{\langle ab\vert \hat{v}\vert ij\rangle}{\epsilon_i+\epsilon_j-\epsilon_a-\epsilon_b}. $$If we compare this to the correlation energy obtained from full configuration interaction theory with a Hartree-Fock basis, we found that
$$ E-E_0 =\Delta E= \sum_{abij}\langle ij | \hat{v}| ab \rangle C_{ij}^{ab}, $$where the energy \( E_0 \) is the reference energy and \( \Delta E \) defines the so-called correlation energy.
We see that if we set
$$ C_{ij}^{ab} =\frac{1}{4}\frac{\langle ab \vert \hat{v} \vert ij \rangle}{\epsilon_i+\epsilon_j-\epsilon_a-\epsilon_b}, $$we have a perfect agreement between FCI and MBPT. However, FCI includes such \( 2p-2h \) correlations to infinite order. In order to make a meaningful comparison we would at least need to sum such correlations to infinite order in perturbation theory.
Summing up, we can see that