Another possibility is to expand the single-particle functions in a known basis and vary the coefficients, that is, the new single-particle wave function is written as a linear expansion in terms of a fixed chosen orthogonal basis (for example the well-known harmonic oscillator functions or the hydrogen-like functions etc). We define our new Hartree-Fock single-particle basis by performing a unitary transformation on our previous basis (labelled with greek indices) as
$$ \begin{equation} \psi_p^{HF} = \sum_{\lambda} C_{p\lambda}\phi_{\lambda}. \label{eq:newbasis} \end{equation} $$In this case we vary the coefficients \( C_{p\lambda} \). If the basis has infinitely many solutions, we need to truncate the above sum. We assume that the basis \( \phi_{\lambda} \) is orthogonal.
It is normal to choose a single-particle basis defined as the eigenfunctions of parts of the full Hamiltonian. The typical situation consists of the solutions of the one-body part of the Hamiltonian, that is we have
$$ \hat{h}_0\phi_{\lambda}=\epsilon_{\lambda}\phi_{\lambda}. $$The single-particle wave functions \( \phi_{\lambda}(\mathbf{r}) \), defined by the quantum numbers \( \lambda \) and \( \mathbf{r} \) are defined as the overlap
$$ \phi_{\lambda}(\mathbf{r}) = \langle \mathbf{r} | \lambda \rangle . $$In deriving the Hartree-Fock equations, we will expand the single-particle functions in a known basis and vary the coefficients, that is, the new single-particle wave function is written as a linear expansion in terms of a fixed chosen orthogonal basis (for example the well-known harmonic oscillator functions or the hydrogen-like functions etc).
We stated that a unitary transformation keeps the orthogonality. To see this consider first a basis of vectors \( \mathbf{v}_i \),
$$ \mathbf{v}_i = \begin{bmatrix} v_{i1} \\ \dots \\ \dots \\v_{in} \end{bmatrix} $$We assume that the basis is orthogonal, that is
$$ \mathbf{v}_j^T\mathbf{v}_i = \delta_{ij}. $$An orthogonal or unitary transformation
$$ \mathbf{w}_i=\mathbf{U}\mathbf{v}_i, $$preserves the dot product and orthogonality since
$$ \mathbf{w}_j^T\mathbf{w}_i=(\mathbf{U}\mathbf{v}_j)^T\mathbf{U}\mathbf{v}_i=\mathbf{v}_j^T\mathbf{U}^T\mathbf{U}\mathbf{v}_i= \mathbf{v}_j^T\mathbf{v}_i = \delta_{ij}. $$This means that if the coefficients \( C_{p\lambda} \) belong to a unitary or orthogonal trasformation (using the Dirac bra-ket notation)
$$ \vert p\rangle = \sum_{\lambda} C_{p\lambda}\vert\lambda\rangle, $$orthogonality is preserved, that is \( \langle \alpha \vert \beta\rangle = \delta_{\alpha\beta} \) and \( \langle p \vert q\rangle = \delta_{pq} \).
This propertry is extremely useful when we build up a basis of many-body Stater determinant based states.
Note also that although a basis \( \vert \alpha\rangle \) contains an infinity of states, for practical calculations we have always to make some truncations.Before we develop the Hartree-Fock equations, there is another very useful property of determinants that we will use both in connection with Hartree-Fock calculations. This applies also to our previous discussion on full configuration interaction theory.
Consider the following determinant
$$ \left| \begin{array}{cc} \alpha_1b_{11}+\alpha_2sb_{12}& a_{12}\\ \alpha_1b_{21}+\alpha_2b_{22}&a_{22}\end{array} \right|=\alpha_1\left|\begin{array}{cc} b_{11}& a_{12}\\ b_{21}&a_{22}\end{array} \right|+\alpha_2\left| \begin{array}{cc} b_{12}& a_{12}\\b_{22}&a_{22}\end{array} \right| $$We can generalize this to an \( n\times n \) matrix and have
$$ \left| \begin{array}{cccccc} a_{11}& a_{12} & \dots & \sum_{k=1}^n c_k b_{1k} &\dots & a_{1n}\\ a_{21}& a_{22} & \dots & \sum_{k=1}^n c_k b_{2k} &\dots & a_{2n}\\ \dots & \dots & \dots & \dots & \dots & \dots \\ \dots & \dots & \dots & \dots & \dots & \dots \\ a_{n1}& a_{n2} & \dots & \sum_{k=1}^n c_k b_{nk} &\dots & a_{nn}\end{array} \right|= \sum_{k=1}^n c_k\left| \begin{array}{cccccc} a_{11}& a_{12} & \dots & b_{1k} &\dots & a_{1n}\\ a_{21}& a_{22} & \dots & b_{2k} &\dots & a_{2n}\\ \dots & \dots & \dots & \dots & \dots & \dots\\ \dots & \dots & \dots & \dots & \dots & \dots\\ a_{n1}& a_{n2} & \dots & b_{nk} &\dots & a_{nn}\end{array} \right| . $$This is a property we will use in our Hartree-Fock discussions.
We can generalize the previous results, now with all elements \( a_{ij} \) being given as functions of linear combinations of various coefficients \( c \) and elements \( b_{ij} \),
$$ \left| \begin{array}{cccccc} \sum_{k=1}^n b_{1k}c_{k1}& \sum_{k=1}^n b_{1k}c_{k2} & \dots & \sum_{k=1}^n b_{1k}c_{kj} &\dots & \sum_{k=1}^n b_{1k}c_{kn}\\ \sum_{k=1}^n b_{2k}c_{k1}& \sum_{k=1}^n b_{2k}c_{k2} & \dots & \sum_{k=1}^n b_{2k}c_{kj} &\dots & \sum_{k=1}^n b_{2k}c_{kn}\\ \dots & \dots & \dots & \dots & \dots & \dots \\ \dots & \dots & \dots & \dots & \dots &\dots \\ \sum_{k=1}^n b_{nk}c_{k1}& \sum_{k=1}^n b_{nk}c_{k2} & \dots & \sum_{k=1}^n b_{nk}c_{kj} &\dots & \sum_{k=1}^n b_{nk}c_{kn}\end{array} \right|=det(\mathbf{C})det(\mathbf{B}), $$where \( det(\mathbf{C}) \) and \( det(\mathbf{B}) \) are the determinants of \( n\times n \) matrices with elements \( c_{ij} \) and \( b_{ij} \) respectively. This is a property we will use in our Hartree-Fock discussions. Convince yourself about the correctness of the above expression by setting \( n=2 \).
With our definition of the new basis in terms of an orthogonal basis we have
$$ \psi_p(x) = \sum_{\lambda} C_{p\lambda}\phi_{\lambda}(x). $$If the coefficients \( C_{p\lambda} \) belong to an orthogonal or unitary matrix, the new basis is also orthogonal. Our Slater determinant in the new basis \( \psi_p(x) \) is written as
$$ \frac{1}{\sqrt{N!}} \left| \begin{array}{ccccc} \psi_{p}(x_1)& \psi_{p}(x_2)& \dots & \dots & \psi_{p}(x_N)\\ \psi_{q}(x_1)&\psi_{q}(x_2)& \dots & \dots & \psi_{q}(x_N)\\ \dots & \dots & \dots & \dots & \dots \\ \dots & \dots & \dots & \dots & \dots \\ \psi_{t}(x_1)&\psi_{t}(x_2)& \dots & \dots & \psi_{t}(x_N)\end{array} \right|=\frac{1}{\sqrt{N!}} \left| \begin{array}{ccccc} \sum_{\lambda} C_{p\lambda}\phi_{\lambda}(x_1)& \sum_{\lambda} C_{p\lambda}\phi_{\lambda}(x_2)& \dots & \dots & \sum_{\lambda} C_{p\lambda}\phi_{\lambda}(x_N)\\ \sum_{\lambda} C_{q\lambda}\phi_{\lambda}(x_1)&\sum_{\lambda} C_{q\lambda}\phi_{\lambda}(x_2)& \dots & \dots & \sum_{\lambda} C_{q\lambda}\phi_{\lambda}(x_N)\\ \dots & \dots & \dots & \dots & \dots \\ \dots & \dots & \dots & \dots & \dots \\ \sum_{\lambda} C_{t\lambda}\phi_{\lambda}(x_1)&\sum_{\lambda} C_{t\lambda}\phi_{\lambda}(x_2)& \dots & \dots & \sum_{\lambda} C_{t\lambda}\phi_{\lambda}(x_N)\end{array} \right|, $$which is nothing but \( det(\mathbf{C})det(\Phi) \), with \( det(\Phi) \) being the determinant given by the basis functions \( \phi_{\lambda}(x) \).
In our discussions hereafter we will use our definitions of single-particle states above and below the Fermi (\( F \)) level given by the labels \( ijkl\dots \le F \) for so-called single-hole states and \( abcd\dots > F \) for so-called particle states. For general single-particle states we employ the labels \( pqrs\dots \).
We have
$$ E[\Phi] = \sum_{\mu=1}^N \langle \mu | h | \mu \rangle + \frac{1}{2}\sum_{{\mu}=1}^N\sum_{{\nu}=1}^N \langle \mu\nu|\hat{v}|\mu\nu\rangle_{AS}, $$we found the expression for the energy functional in terms of the basis function \( \phi_{\lambda}(\mathbf{r}) \). We then varied the above energy functional with respect to the basis functions \( |\mu \rangle \).
Now we are interested in defining a new basis defined in terms of a chosen basis as defined in Eq. \eqref{eq:newbasis}. We can then rewrite the energy functional as
$$ \begin{equation} E[\Phi^{HF}] = \sum_{i=1}^N \langle i | h | i \rangle + \frac{1}{2}\sum_{ij=1}^N\langle ij|\hat{v}|ij\rangle_{AS}, \label{FunctionalEPhi2} \end{equation} $$where \( \Phi^{HF} \) is the new Slater determinant defined by the new basis of Eq. \eqref{eq:newbasis}.
Using Eq. \eqref{eq:newbasis} we can rewrite Eq. \eqref{FunctionalEPhi2} as
$$ \begin{equation} E[\Psi] = \sum_{i=1}^N \sum_{\alpha\beta} C^*_{i\alpha}C_{i\beta}\langle \alpha | h | \beta \rangle + \frac{1}{2}\sum_{ij=1}^N\sum_{{\alpha\beta\gamma\delta}} C^*_{i\alpha}C^*_{j\beta}C_{i\gamma}C_{j\delta}\langle \alpha\beta|\hat{v}|\gamma\delta\rangle_{AS}. \label{FunctionalEPhi3} \end{equation} $$We wish now to minimize the above functional. We introduce again a set of Lagrange multipliers, noting that since \( \langle i | j \rangle = \delta_{i,j} \) and \( \langle \alpha | \beta \rangle = \delta_{\alpha,\beta} \), the coefficients \( C_{i\gamma} \) obey the relation
$$ \langle i | j \rangle=\delta_{i,j}=\sum_{\alpha\beta} C^*_{i\alpha}C_{i\beta}\langle \alpha | \beta \rangle= \sum_{\alpha} C^*_{i\alpha}C_{i\alpha}, $$which allows us to define a functional to be minimized that reads
$$ \begin{equation} F[\Phi^{HF}]=E[\Phi^{HF}] - \sum_{i=1}^N\epsilon_i\sum_{\alpha} C^*_{i\alpha}C_{i\alpha}. \label{_auto1} \end{equation} $$Minimizing with respect to \( C^*_{i\alpha} \), remembering that the equations for \( C^*_{i\alpha} \) and \( C_{i\alpha} \) can be written as two independent equations, we obtain
$$ \frac{d}{dC^*_{i\alpha}}\left[ E[\Phi^{HF}] - \sum_{j}\epsilon_j\sum_{\alpha} C^*_{j\alpha}C_{j\alpha}\right]=0, $$which yields for every single-particle state \( i \) and index \( \alpha \) (recalling that the coefficients \( C_{i\alpha} \) are matrix elements of a unitary (or orthogonal for a real symmetric matrix) matrix) the following Hartree-Fock equations
$$ \sum_{\beta} C_{i\beta}\langle \alpha | h | \beta \rangle+ \sum_{j=1}^N\sum_{\beta\gamma\delta} C^*_{j\beta}C_{j\delta}C_{i\gamma}\langle \alpha\beta|\hat{v}|\gamma\delta\rangle_{AS}=\epsilon_i^{HF}C_{i\alpha}. $$We can rewrite this equation as (changing dummy variables)
$$ \sum_{\beta} \left\{\langle \alpha | h | \beta \rangle+ \sum_{j}^N\sum_{\gamma\delta} C^*_{j\gamma}C_{j\delta}\langle \alpha\gamma|\hat{v}|\beta\delta\rangle_{AS}\right\}C_{i\beta}=\epsilon_i^{HF}C_{i\alpha}. $$Note that the sums over greek indices run over the number of basis set functions (in principle an infinite number).
Defining
$$ h_{\alpha\beta}^{HF}=\langle \alpha | h | \beta \rangle+ \sum_{j=1}^N\sum_{\gamma\delta} C^*_{j\gamma}C_{j\delta}\langle \alpha\gamma|\hat{v}|\beta\delta\rangle_{AS}, $$we can rewrite the new equations as
$$ \begin{equation} \sum_{\beta}h_{\alpha\beta}^{HF}C_{i\beta}=\epsilon_i^{HF}C_{i\alpha}. \label{eq:newhf} \end{equation} $$The latter is nothing but a standard eigenvalue problem.
We see that we do not need to compute any integrals in an iterative procedure for solving the equations. It suffices to tabulate the matrix elements \( \langle \alpha | h | \beta \rangle \) and \( \langle \alpha\gamma|\hat{v}|\beta\delta\rangle_{AS} \) once and for all. Successive iterations require thus only a look-up in tables over one-body and two-body matrix elements. These details will be discussed below when we solve the Hartree-Fock equations numerical.
Our Hartree-Fock matrix is thus
$$ \hat{h}_{\alpha\beta}^{HF}=\langle \alpha | \hat{h}_0 | \beta \rangle+ \sum_{j=1}^N\sum_{\gamma\delta} C^*_{j\gamma}C_{j\delta}\langle \alpha\gamma|\hat{v}|\beta\delta\rangle_{AS}. $$The Hartree-Fock equations are solved in an iterative waym starting with a guess for the coefficients \( C_{j\gamma}=\delta_{j,\gamma} \) and solving the equations by diagonalization till the new single-particle energies \( \epsilon_i^{\mathrm{HF}} \) do not change anymore by a prefixed quantity.
Normally we assume that the single-particle basis \( |\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|\hat{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 Hartree-Fock equations are, in their simplest form, solved in an iterative way, starting with a guess for the coefficients \( C_{i\alpha} \). We label the coefficients as \( C_{i\alpha}^{(n)} \), where the subscript \( n \) stands for iteration \( n \). To set up the algorithm we can proceed as follows:
Solving the Hartree-Fock eigenvalue problem yields then new eigenvectors \( C_{i\alpha}^{(1)} \) and eigenvalues \( \epsilon_i^{HF(1)} \).
The diagonalization with the new Hartree-Fock potential yields new eigenvectors and eigenvalues. This process is continued till for example
$$ \frac{\sum_{p} |\epsilon_i^{(n)}-\epsilon_i^{(n-1)}|}{m} \le \lambda, $$where \( \lambda \) is a user prefixed quantity (\( \lambda \sim 10^{-8} \) or smaller) and \( p \) runs over all calculated single-particle energies and \( m \) is the number of single-particle states.
The equations are often rewritten in terms of a so-called density matrix, which is defined as
$$ \begin{equation} \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}. \label{_auto2} \end{equation} $$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 \).
import numpy as np
class HartreeFock:
def __init__(self, num_electrons, num_orbitals):
self.num_electrons = num_electrons
self.num_orbitals = num_orbitals
# You would need to replace these by proper integrals
self.h = np.random.rand(num_orbitals, num_orbitals) # One-electron integrals
self.coulomb = np.random.rand(num_orbitals, num_orbitals, num_orbitals, num_orbitals) # Two-electron integrals
def build_fock_matrix(self, density_matrix):
fock_matrix = self.h.copy()
for i in range(self.num_orbitals):
for j in range(self.num_orbitals):
fock_matrix[i, j] += np.sum(density_matrix * self.coulomb[i, j])
return fock_matrix
def build_density_matrix(self, coefficients):
density_matrix = np.zeros((self.num_orbitals, self.num_orbitals))
for i in range(self.num_electrons):
density_matrix += np.outer(coefficients[:, i], coefficients[:, i])
return density_matrix
def diagonalize(self, fock_matrix):
energy, coefficients = np.linalg.eigh(fock_matrix)
return energy, coefficients
def run(self, max_iter=100, tol=1e-6):
coeffs = np.zeros((self.num_orbitals, self.num_electrons))
density_matrix = np.zeros((self.num_orbitals, self.num_orbitals))
for iteration in range(max_iter):
fock_matrix = self.build_fock_matrix(density_matrix)
energies, coeffs = self.diagonalize(fock_matrix)
new_density_matrix = self.build_density_matrix(coeffs)
if np.linalg.norm(new_density_matrix - density_matrix) < tol:
print(f"Converged in {iteration} iterations")
break
density_matrix = new_density_matrix
return energies, coeffs
# Example usage
if __name__ == "__main__":
num_electrons = 2
num_orbitals = 4
hf = HartreeFock(num_electrons, num_orbitals)
energies, coefficients = hf.run()
print("Final energies:", energies)
print("Final coefficients:", coefficients)
We wish now to derive the Hartree-Fock equations using our second-quantized formalism and study the stability of the equations. Our ansatz for the ground state of the system is approximated as (this is our representation of a Slater determinant in second quantization)
$$ |\Phi_0\rangle = |c\rangle = a^{\dagger}_i a^{\dagger}_j \dots a^{\dagger}_l|0\rangle. $$We wish to determine \( \hat{u}^{HF} \) so that \( E_0^{HF}= \langle c|\hat{H}| c\rangle \) becomes a local minimum.
In our analysis here we will need Thouless' theorem, which states that an arbitrary Slater determinant \( |c'\rangle \) which is not orthogonal to a determinant \( | c\rangle ={\displaystyle\prod_{i=1}^{n}} a_{\alpha_{i}}^{\dagger}|0\rangle \), can be written as
$$ |c'\rangle=exp\left\{\sum_{a>F}\sum_{i\le F}C_{ai}a_{a}^{\dagger}a_{i}\right\}| c\rangle $$Let us give a simple proof of Thouless' theorem. The theorem states that we can make a linear combination av particle-hole excitations with respect to a given reference state \( \vert c\rangle \). With this linear combination, we can make a new Slater determinant \( \vert c'\rangle \) which is not orthogonal to \( \vert c\rangle \), that is
$$ \langle c|c'\rangle \ne 0. $$To show this we need some intermediate steps. The exponential product of two operators \( \exp{\hat{A}}\times\exp{\hat{B}} \) is equal to \( \exp{(\hat{A}+\hat{B})} \) only if the two operators commute, that is
$$ [\hat{A},\hat{B}] = 0. $$If the operators do not commute, we need to resort to the Baker-Campbell-Hauersdorf. This relation states that
$$ \exp{\hat{C}}=\exp{\hat{A}}\exp{\hat{B}}, $$with
$$ \hat{C}=\hat{A}+\hat{B}+\frac{1}{2}[\hat{A},\hat{B}]+\frac{1}{12}[[\hat{A},\hat{B}],\hat{B}]-\frac{1}{12}[[\hat{A},\hat{B}],\hat{A}]+\dots $$From these relations, we note that in our expression for \( |c'\rangle \) we have commutators of the type
$$ [a_{a}^{\dagger}a_{i},a_{b}^{\dagger}a_{j}], $$and it is easy to convince oneself that these commutators, or higher powers thereof, are all zero. This means that we can write out our new representation of a Slater determinant as
$$ |c'\rangle=exp\left\{\sum_{a>F}\sum_{i\le F}C_{ai}a_{a}^{\dagger}a_{i}\right\}| c\rangle=\prod_{i}\left\{1+\sum_{a>F}C_{ai}a_{a}^{\dagger}a_{i}+\left(\sum_{a>F}C_{ai}a_{a}^{\dagger}a_{i}\right)^2+\dots\right\}| c\rangle $$We note that
$$ \prod_{i}\sum_{a>F}C_{ai}a_{a}^{\dagger}a_{i}\sum_{b>F}C_{bi}a_{b}^{\dagger}a_{i}| c\rangle =0, $$and all higher-order powers of these combinations of creation and annihilation operators disappear due to the fact that \( (a_i)^n| c\rangle =0 \) when \( n > 1 \). This allows us to rewrite the expression for \( |c'\rangle \) as
$$ |c'\rangle=\prod_{i}\left\{1+\sum_{a>F}C_{ai}a_{a}^{\dagger}a_{i}\right\}| c\rangle, $$which we can rewrite as
$$ |c'\rangle=\prod_{i}\left\{1+\sum_{a>F}C_{ai}a_{a}^{\dagger}a_{i}\right\}| a^{\dagger}_{i_1} a^{\dagger}_{i_2} \dots a^{\dagger}_{i_n}|0\rangle. $$The last equation can be written as
$$ \begin{align} |c'\rangle&=\prod_{i}\left\{1+\sum_{a>F}C_{ai}a_{a}^{\dagger}a_{i}\right\}| a^{\dagger}_{i_1} a^{\dagger}_{i_2} \dots a^{\dagger}_{i_n}|0\rangle=\left(1+\sum_{a>F}C_{ai_1}a_{a}^{\dagger}a_{i_1}\right)a^{\dagger}_{i_1} \label{_auto3}\\ & \times\left(1+\sum_{a>F}C_{ai_2}a_{a}^{\dagger}a_{i_2}\right)a^{\dagger}_{i_2} \dots |0\rangle=\prod_{i}\left(a^{\dagger}_{i}+\sum_{a>F}C_{ai}a_{a}^{\dagger}\right)|0\rangle. \label{_auto4} \end{align} $$If we define a new creation operator
$$ \begin{equation} b^{\dagger}_{i}=a^{\dagger}_{i}+\sum_{a>F}C_{ai}a_{a}^{\dagger}, \label{eq:newb} \end{equation} $$we have
$$ |c'\rangle=\prod_{i}b^{\dagger}_{i}|0\rangle=\prod_{i}\left(a^{\dagger}_{i}+\sum_{a>F}C_{ai}a_{a}^{\dagger}\right)|0\rangle, $$meaning that the new representation of the Slater determinant in second quantization, \( |c'\rangle \), looks like our previous ones. However, this representation is not general enough since we have a restriction on the sum over single-particle states in Eq. \eqref{eq:newb}. The single-particle states have all to be above the Fermi level.
The question then is whether we can construct a general representation of a Slater determinant with a creation operator
$$ \tilde{b}^{\dagger}_{i}=\sum_{p}f_{ip}a_{p}^{\dagger}, $$where \( f_{ip} \) is a matrix element of a unitary matrix which transforms our creation and annihilation operators \( a^{\dagger} \) and \( a \) to \( \tilde{b}^{\dagger} \) and \( \tilde{b} \). These new operators define a new representation of a Slater determinant as
$$ |\tilde{c}\rangle=\prod_{i}\tilde{b}^{\dagger}_{i}|0\rangle. $$We need to show that \( |\tilde{c}\rangle= |c'\rangle \). We need also to assume that the new state is not orthogonal to \( |c\rangle \), that is \( \langle c| \tilde{c}\rangle \ne 0 \). From this it follows that
$$ \langle c| \tilde{c}\rangle=\langle 0| a_{i_n}\dots a_{i_1}\left(\sum_{p=i_1}^{i_n}f_{i_1p}a_{p}^{\dagger} \right)\left(\sum_{q=i_1}^{i_n}f_{i_2q}a_{q}^{\dagger} \right)\dots \left(\sum_{t=i_1}^{i_n}f_{i_nt}a_{t}^{\dagger} \right)|0\rangle, $$which is nothing but the determinant \( det(f_{ip}) \) which we can, using the intermediate normalization condition, normalize to one, that is
$$ det(f_{ip})=1, $$meaning that \( f \) has an inverse defined as (since we are dealing with orthogonal, and in our case unitary as well, transformations)
$$ \sum_{k} f_{ik}f^{-1}_{kj} = \delta_{ij}, $$and
$$ \sum_{j} f^{-1}_{ij}f_{jk} = \delta_{ik}. $$Using these relations we can then define the linear combination of creation (and annihilation as well) operators as
$$ \sum_{i}f^{-1}_{ki}\tilde{b}^{\dagger}_{i}=\sum_{i}f^{-1}_{ki}\sum_{p=i_1}^{\infty}f_{ip}a_{p}^{\dagger}=a_{k}^{\dagger}+\sum_{i}\sum_{p=i_{n+1}}^{\infty}f^{-1}_{ki}f_{ip}a_{p}^{\dagger}. $$Defining
$$ c_{kp}=\sum_{i \le F}f^{-1}_{ki}f_{ip}, $$we can redefine
$$ a_{k}^{\dagger}+\sum_{i}\sum_{p=i_{n+1}}^{\infty}f^{-1}_{ki}f_{ip}a_{p}^{\dagger}=a_{k}^{\dagger}+\sum_{p=i_{n+1}}^{\infty}c_{kp}a_{p}^{\dagger}=b_k^{\dagger}, $$our starting point.
We have shown that our general representation of a Slater determinant
$$ |\tilde{c}\rangle=\prod_{i}\tilde{b}^{\dagger}_{i}|0\rangle=|c'\rangle=\prod_{i}b^{\dagger}_{i}|0\rangle, $$with
$$ b_k^{\dagger}=a_{k}^{\dagger}+\sum_{p=i_{n+1}}^{\infty}c_{kp}a_{p}^{\dagger}. $$This means that we can actually write an ansatz for the ground state of the system as a linear combination of terms which contain the ansatz itself \( |c\rangle \) with an admixture from an infinity of one-particle-one-hole states. The latter has important consequences when we wish to interpret the Hartree-Fock equations and their stability. We can rewrite the new representation as
$$ |c'\rangle = |c\rangle+|\delta c\rangle, $$where \( |\delta c\rangle \) can now be interpreted as a small variation. If we approximate this term with contributions from one-particle-one-hole (1p-1h) states only, we arrive at
$$ |c'\rangle = \left(1+\sum_{ai}\delta C_{ai}a_{a}^{\dagger}a_i\right)|c\rangle. $$In our derivation of the Hartree-Fock equations we have shown that
$$ \langle \delta c| \hat{H} | c\rangle =0, $$which means that we have to satisfy
$$ \langle c|\sum_{ai}\delta C_{ai}\left\{a_{a}^{\dagger}a_i\right\} \hat{H} | c\rangle =0. $$With this as a background, we are now ready to study the stability of the Hartree-Fock equations. This is the topic for week 40.
The variational condition for deriving the Hartree-Fock equations guarantees only that the expectation value \( \langle c | \hat{H} | c \rangle \) has an extreme value, not necessarily a minimum. To figure out whether the extreme value we have found is a minimum, we can use second quantization to analyze our results and find a criterion for the above expectation value to a local minimum. We will use Thouless' theorem and show that
$$ \frac{\langle c' |\hat{H} | c'\rangle}{\langle c' |c'\rangle} \ge \langle c |\hat{H} | c\rangle= E_0, $$with
$$ {|c'\rangle} = {|c\rangle + |\delta c\rangle}. $$Using Thouless' theorem we can write out \( |c'\rangle \) as
$$ \begin{align} {|c'\rangle}&=\exp\left\{\sum_{a > F}\sum_{i \le F}\delta C_{ai}a_{a}^{\dagger}a_{i}\right\}| c\rangle \label{_auto5}\\ &=\left\{1+\sum_{a > F}\sum_{i \le F}\delta C_{ai}a_{a}^{\dagger} a_{i}+\frac{1}{2!}\sum_{ab > F}\sum_{ij \le F}\delta C_{ai}\delta C_{bj}a_{a}^{\dagger}a_{i}a_{b}^{\dagger}a_{j}+\dots\right\} \label{_auto6} \end{align} $$where the amplitudes \( \delta C \) are small.
The norm of \( |c'\rangle \) is given by (using the intermediate normalization condition \( \langle c' |c\rangle=1 \))
$$ \langle c' | c'\rangle = 1+\sum_{a>F} \sum_{i\le F}|\delta C_{ai}|^2+O(\delta C_{ai}^3). $$The expectation value for the energy is now given by (using the Hartree-Fock condition)
$$ \langle c' |\hat{H} | c'\rangle=\langle c |\hat{H} | c\rangle + \sum_{ab>F} \sum_{ij\le F}\delta C_{ai}^*\delta C_{bj}\langle c |a_{i}^{\dagger}a_{a}\hat{H}a_{b}^{\dagger}a_{j}|c\rangle+ $$ $$ \frac{1}{2!}\sum_{ab>F} \sum_{ij\le F}\delta C_{ai}\delta C_{bj}\langle c |\hat{H}a_{a}^{\dagger}a_{i}a_{b}^{\dagger}a_{j}|c\rangle+\frac{1}{2!}\sum_{ab>F} \sum_{ij\le F}\delta C_{ai}^*\delta C_{bj}^*\langle c|a_{j}^{\dagger}a_{b}a_{i}^{\dagger}a_{a}\hat{H}|c\rangle +\dots $$We have already calculated the second term on the right-hand side of the previous equation
$$ \begin{align} \langle c | \left(\{a^\dagger_i a_a\} \hat{H} \{a^\dagger_b a_j\} \right) | c\rangle&=\sum_{pq} \sum_{ijab}\delta C_{ai}^*\delta C_{bj} \langle p|\hat{h}_0 |q\rangle \langle c | \left(\{a^{\dagger}_i a_a\}\{a^{\dagger}_pa_q\} \{a^{\dagger}_b a_j\} \right)| c\rangle \label{_auto7}\\ & +\frac{1}{4} \sum_{pqrs} \sum_{ijab}\delta C_{ai}^*\delta C_{bj} \langle pq| \hat{v}|rs\rangle \langle c | \left(\{a^\dagger_i a_a\}\{a^{\dagger}_p a^{\dagger}_q a_s a_r\} \{a^{\dagger}_b a_j\} \right)| c\rangle , \label{_auto8} \end{align} $$resulting in
$$ E_0\sum_{ai}|\delta C_{ai}|^2+\sum_{ai}|\delta C_{ai}|^2(\varepsilon_a-\varepsilon_i)-\sum_{ijab} \langle aj|\hat{v}| bi\rangle \delta C_{ai}^*\delta C_{bj}. $$ $$ \frac{1}{2!}\langle c |\left(\{a^\dagger_j a_b\} \{a^\dagger_i a_a\} \hat{V}_N \right) | c\rangle = \frac{1}{2!}\langle c |\left( \hat{V}_N \{a^\dagger_a a_i\} \{a^\dagger_b a_j\} \right)^{\dagger} | c\rangle $$which is nothing but
$$ \frac{1}{2!}\langle c | \left( \hat{V}_N \{a^\dagger_a a_i\} \{a^\dagger_b a_j\} \right) | c\rangle^* =\frac{1}{2} \sum_{ijab} (\langle ij|\hat{v}|ab\rangle)^*\delta C_{ai}^*\delta C_{bj}^* $$or
$$ \frac{1}{2} \sum_{ijab} (\langle ab|\hat{v}|ij\rangle)\delta C_{ai}^*\delta C_{bj}^* $$where we have used the relation
$$ \langle a |\hat{A} | b\rangle = (\langle b |\hat{A}^{\dagger} | a\rangle)^* $$due to the hermiticity of \( \hat{H} \) and \( \hat{V} \).
We define two matrix elements
$$ A_{ai,bj}=-\langle aj|\hat{v} bi\rangle $$and
$$ B_{ai,bj}=\langle ab|\hat{v}|ij\rangle $$both being anti-symmetrized.
With these definitions we write out the energy as
$$ \begin{align} \langle c'|H|c'\rangle& = \left(1+\sum_{ai}|\delta C_{ai}|^2\right)\langle c |H|c\rangle+\sum_{ai}|\delta C_{ai}|^2(\varepsilon_a^{HF}-\varepsilon_i^{HF})+\sum_{ijab}A_{ai,bj}\delta C_{ai}^*\delta C_{bj}+ \label{_auto9}\\ &\frac{1}{2} \sum_{ijab} B_{ai,bj}^*\delta C_{ai}\delta C_{bj}+\frac{1}{2} \sum_{ijab} B_{ai,bj}\delta C_{ai}^*\delta C_{bj}^* +O(\delta C_{ai}^3), \label{_auto10} \end{align} $$which can be rewritten as
$$ \langle c'|H|c'\rangle = \left(1+\sum_{ai}|\delta C_{ai}|^2\right)\langle c |H|c\rangle+\Delta E+O(\delta C_{ai}^3), $$and skipping higher-order terms we arrived
$$ \frac{\langle c' |\hat{H} | c'\rangle}{\langle c' |c'\rangle} =E_0+\frac{\Delta E}{\left(1+\sum_{ai}|\delta C_{ai}|^2\right)}. $$We have defined
$$ \Delta E = \frac{1}{2} \langle \chi | \hat{M}| \chi \rangle $$with the vectors
$$ \chi = \left[ \delta C\hspace{0.2cm} \delta C^*\right]^T $$and the matrix
$$ \hat{M}=\left(\begin{array}{cc} \Delta + A & B \\ B^* & \Delta + A^*\end{array}\right), $$with \( \Delta_{ai,bj} = (\varepsilon_a-\varepsilon_i)\delta_{ab}\delta_{ij} \).
The condition
$$ \Delta E = \frac{1}{2} \langle \chi | \hat{M}| \chi \rangle \ge 0 $$for an arbitrary vector
$$ \chi = \left[ \delta C\hspace{0.2cm} \delta C^*\right]^T $$means that all eigenvalues of the matrix have to be larger than or equal zero. A necessary (but no sufficient) condition is that the matrix elements (for all \( ai \) )
$$ (\varepsilon_a-\varepsilon_i)\delta_{ab}\delta_{ij}+A_{ai,bj} \ge 0. $$This equation can be used as a first test of the stability of the Hartree-Fock equation.
The electron gas is perhaps the only realistic model of a system of many interacting particles that allows for an analytical solution of the Hartree-Fock equations. Furthermore, to first order in the interaction, one can also obtain an analytical expression for the total energy and several other properties of a many-particle systems. The model gives a very good approximation to the properties of valence electrons in metals. The assumptions are
The homogeneus electron gas is a system of electrons that is not influenced by external forces except by an attraction provided by a uniform background of ions. These ions give rise to a uniform background charge. The ions are stationary and the system as a whole is neutral. Irrespective of this simplicity, this system, in both two and three-dimensions, has eluded a proper description of correlations in terms of various first principle methods, except perhaps for quantum Monte Carlo methods. In particular, the diffusion Monte Carlo calculations of Ceperley and Ceperley and Tanatar are presently still considered as the best possible benchmarks for the two- and three-dimensional electron gas.
The electron gas, in two or three dimensions is thus interesting as a test-bed for electron-electron correlations. The three-dimensional electron gas is particularly important as a cornerstone of the local-density approximation in density-functional theory. In the physical world, systems similar to the three-dimensional electron gas can be found in, for example, alkali metals and doped semiconductors. Two-dimensional electron fluids are observed on metal and liquid-helium surfaces, as well as at metal-oxide-semiconductor interfaces. However, the Coulomb interaction has an infinite range, and therefore long-range correlations play an essential role in the electron gas.
At low densities, the electrons become localized and form a lattice. This so-called Wigner crystallization is a direct consequence of the long-ranged repulsive interaction. At higher densities, the electron gas is better described as a liquid. When using, for example, Monte Carlo methods the electron gas must be approximated by a finite system. The long-range Coulomb interaction in the electron gas causes additional finite-size effects that are not present in other infinite systems like nuclear matter or neutron star matter. This poses additional challenges to many-body methods when applied to the electron gas.
This is a homogeneous system and the one-particle wave functions are given by plane wave functions normalized to a volume \( \Omega \) 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{\Omega}}\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 electrons 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 can be written as
$$ \hat{T}=\sum_{\mathbf{k}\sigma}\frac{\hbar^2k^2}{2m}a_{\mathbf{k}\sigma}^{\dagger}a_{\mathbf{k}\sigma}. $$The Hamiltonian operator is given by
$$ \hat{H}=\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/\Omega \). 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.
In the first exercise below we show that the Hartree-Fock energy can be written as
$$ \varepsilon_{k}^{HF}=\frac{\hbar^{2}k^{2}}{2m_e}-\frac{e^{2}} {\Omega^{2}}\sum_{k'\leq k_{F}}\int d\mathbf{r}e^{i(\mathbf{k}'-\mathbf{k})\mathbf{r}}\int d\mathbf{r'}\frac{e^{i(\mathbf{k}-\mathbf{k}')\mathbf{r}'}} {\vert\mathbf{r}-\mathbf{r}'\vert} $$resulting in
$$ \varepsilon_{k}^{HF}=\frac{\hbar^{2}k^{2}}{2m_e}-\frac{e^{2} k_{F}}{2\pi} \left[ 2+\frac{k_{F}^{2}-k^{2}}{kk_{F}}ln\left\vert\frac{k+k_{F}} {k-k_{F}}\right\vert \right] $$The previous result can be rewritten in terms of the density
$$ n= \frac{k_F^3}{3\pi^2}=\frac{3}{4\pi r_s^3}, $$where \( n=N_e/\Omega \), \( N_e \) being the number of electrons, and \( r_s \) is the radius of a sphere which represents the volum per conducting electron. It can be convenient to use the Bohr radius \( a_0=\hbar^2/e^2m_e \). For most metals we have a relation \( r_s/a_0\sim 2-6 \). The quantity \( r_s \) is dimensionless.
In the second exercise below we find that the total energy \( E_0/N_e=\langle\Phi_{0}|\hat{H}|\Phi_{0}\rangle/N_e \) for for this system to first order in the interaction is given as
$$ E_0/N_e=\frac{e^2}{2a_0}\left[\frac{2.21}{r_s^2}-\frac{0.916}{r_s}\right]. $$The electron gas model allows closed form solutions for quantities like the single-particle Hartree-Fock energy. The latter quantity is given by the following expression
$$ \varepsilon_{k}^{HF}=\frac{\hbar^{2}k^{2}}{2m}-\frac{e^{2}} {V^{2}}\sum_{k'\leq k_{F}}\int d\mathbf{r}e^{i(\mathbf{k'}-\mathbf{k})\mathbf{r}}\int d\mathbf{r}'\frac{e^{i(\mathbf{k}-\mathbf{k'})\mathbf{r}'}} {\vert\mathbf{r}-\mathbf{r'}\vert} $$a) Show first that
$$ \varepsilon_{k}^{HF}=\frac{\hbar^{2}k^{2}}{2m}-\frac{e^{2} k_{F}}{2\pi} \left[ 2+\frac{k_{F}^{2}-k^{2}}{kk_{F}}ln\left\vert\frac{k+k_{F}} {k-k_{F}}\right\vert \right] $$
Hint.
Hint: Introduce the convergence factor \( e^{-\mu\vert\mathbf{r}-\mathbf{r}'\vert} \) in the potential and use \( \sum_{\mathbf{k}}\rightarrow \frac{V}{(2\pi)^{3}}\int d\mathbf{k} \)
Solution.
We want to show that, given the Hartree-Fock equation for the electron gas
$$ \varepsilon_{k}^{HF}=\frac{\hbar^{2}k^{2}}{2m}-\frac{e^{2}} {V^{2}}\sum_{p\leq k_{F}}\int d\mathbf{r}\exp{(i(\mathbf{p}-\mathbf{k})\mathbf{r})}\int d\mathbf{r}'\frac{\exp{(i(\mathbf{k}-\mathbf{p})\mathbf{r}'})} {\vert\mathbf{r}-\mathbf{r'}\vert} $$the single-particle energy can be written as
$$ \varepsilon_{k}^{HF}=\frac{\hbar^{2}k^{2}}{2m}-\frac{e^{2} k_{F}}{2\pi} \left[ 2+\frac{k_{F}^{2}-k^{2}}{kk_{F}}ln\left\vert\frac{k+k_{F}} {k-k_{F}}\right\vert \right]. $$We introduce the convergence factor \( e^{-\mu\vert\mathbf{r}-\mathbf{r}'\vert} \) in the potential and use \( \sum_{\mathbf{k}}\rightarrow \frac{V}{(2\pi)^{3}}\int d\mathbf{k} \). We can then rewrite the integral as
$$ \begin{align} \frac{e^{2}} {V^{2}}\sum_{k'\leq k_{F}}\int d\mathbf{r}\exp{(i(\mathbf{k'}-\mathbf{k})\mathbf{r})}\int d\mathbf{r}'\frac{\exp{(i(\mathbf{k}-\mathbf{p})\mathbf{r}'})} {\vert\mathbf{r}-\mathbf{r'}\vert}= & \label{_auto11}\\ \frac{e^{2}}{V (2\pi)^3} \int d\mathbf{r}\int \frac{d\mathbf{r}'}{\vert\mathbf{r}-\mathbf{r'}\vert}\exp{(-i\mathbf{k}(\mathbf{r}-\mathbf{r}'))}\int d\mathbf{p}\exp{(i\mathbf{p}(\mathbf{r}-\mathbf{r}'))}, \label{_auto12} \end{align} $$and introducing the abovementioned convergence factor we have
$$ \begin{align} \lim_{\mu \to 0}\frac{e^{2}}{V (2\pi)^3} \int d\mathbf{r}\int d\mathbf{r}'\frac{\exp{(-\mu\vert\mathbf{r}-\mathbf{r}'\vert})}{\vert\mathbf{r}-\mathbf{r'}\vert}\int d\mathbf{p}\exp{(i(\mathbf{p}-\mathbf{k})(\mathbf{r}-\mathbf{r}'))}. \label{_auto13} \end{align} $$With a change variables to \( \mathbf{x} = \mathbf{r}-\mathbf{r}' \) and \( \mathbf{y}=\mathbf{r}' \) we rewrite the last integral as
$$ \lim_{\mu \to 0}\frac{e^{2}}{V (2\pi)^3} \int d\mathbf{p}\int d\mathbf{y}\int d\mathbf{x}\exp{(i(\mathbf{p}-\mathbf{k})\mathbf{x})}\frac{\exp{(-\mu\vert\mathbf{x}\vert})}{\vert\mathbf{x}\vert}. $$The integration over \( \mathbf{x} \) can be performed using spherical coordinates, resulting in (with \( x=\vert \mathbf{x}\vert \))
$$ \int d\mathbf{x}\exp{(i(\mathbf{p}-\mathbf{k})\mathbf{x})}\frac{\exp{(-\mu\vert\mathbf{x}\vert})}{\vert\mathbf{x}\vert}=\int x^2 dx d\phi d\cos{(\theta)}\exp{(i(\mathbf{p}-\mathbf{k})x\cos{(\theta))}}\frac{\exp{(-\mu x)}}{x}. $$We obtain
$$ \begin{align} 4\pi \int dx \frac{ \sin{(\vert \mathbf{p}-\mathbf{k}\vert)x} }{\vert \mathbf{p}-\mathbf{k}\vert}{\exp{(-\mu x)}}= \frac{4\pi}{\mu^2+\vert \mathbf{p}-\mathbf{k}\vert^2}. \label{_auto14} \end{align} $$This results gives us
$$ \begin{align} \lim_{\mu \to 0}\frac{e^{2}}{V (2\pi)^3} \int d\mathbf{p}\int d\mathbf{y}\frac{4\pi}{\mu^2+\vert \mathbf{p}-\mathbf{k}\vert^2}=\lim_{\mu \to 0}\frac{e^{2}}{ 2\pi^2} \int d\mathbf{p}\frac{1}{\mu^2+\vert \mathbf{p}-\mathbf{k}\vert^2}, \label{_auto15} \end{align} $$where we have used that the integrand on the left-hand side does not depend on \( \mathbf{y} \) and that \( \int d\mathbf{y}=V \).
Introducing spherical coordinates we can rewrite the integral as
$$ \begin{align} \lim_{\mu \to 0}\frac{e^{2}}{ 2\pi^2} \int d\mathbf{p}\frac{1}{\mu^2+\vert \mathbf{p}-\mathbf{k}\vert^2}=\frac{e^{2}}{ 2\pi^2} \int d\mathbf{p}\frac{1}{\vert \mathbf{p}-\mathbf{k}\vert^2}=& \label{_auto16}\\ \frac{e^{2}}{\pi} \int_0^{k_F} p^2dp\int_0^{\pi} d\theta\cos{(\theta)}\frac{1}{p^2+k^2-2pk\cos{(\theta)}}, \label{_auto17} \end{align} $$and with the change of variables \( \cos{(\theta)}=u \) we have
$$ \frac{e^{2}}{\pi} \int_0^{k_F} p^2dp\int_{0}^{\pi} d\theta\cos{(\theta)}\frac{1}{p^2+k^2-2pk\cos{(\theta)}}=\frac{e^{2}}{\pi} \int_0^{k_F} p^2dp\int_{-1}^{1} du\frac{1}{p^2+k^2-2pku}, $$which gives
$$ \frac{e^{2}}{k\pi} \int_0^{k_F} pdp\left\{ln(\vert p+k\vert)-ln(\vert p-k\vert)\right\}. $$Introducing new variables \( x=p+k \) and \( y=p-k \), we obtain after some straightforward reordering of the integral
$$ \frac{e^{2}}{k\pi}\left[ kk_F+\frac{k_{F}^{2}-k^{2}}{kk_{F}}ln\left\vert\frac{k+k_{F}} {k-k_{F}}\right\vert \right], $$which gives the abovementioned expression for the single-particle energy.
b) Rewrite the above result as a function of the density
$$ n= \frac{k_F^3}{3\pi^2}=\frac{3}{4\pi r_s^3}, $$where \( n=N/V \), \( N \) being the number of particles, and \( r_s \) is the radius of a sphere which represents the volum per conducting electron.
Solution.
Introducing the dimensionless quantity \( x=k/k_F \) and the function
$$ F(x) = \frac{1}{2}+\frac{1-x^2}{4x}\ln{\left\vert \frac{1+x}{1-x}\right\vert}, $$we can rewrite the single-particle Hartree-Fock energy as
$$ \varepsilon_{k}^{HF}=\frac{\hbar^{2}k^{2}}{2m}-\frac{2e^{2} k_{F}}{\pi}F(k/k_F), $$and dividing by the non-interacting contribution at the Fermi level,
$$ \varepsilon_{0}^{F}=\frac{\hbar^{2}k_F^{2}}{2m}, $$we have
$$ \frac{\varepsilon_{k}^{HF} }{\varepsilon_{0}^{F}}=x^2-\frac{e^2m}{\hbar^2 k_F\pi}F(x)=x^2-\frac{4}{\pi k_Fa_0}F(x), $$where \( a_0=0.0529 \) nm is the Bohr radius, setting thereby a natural length scale.
By introducing the radius \( r_s \) of a sphere whose volume is the volume occupied by each electron, we can rewrite the previous equation in terms of \( r_s \) using that the electron density \( n=N/V \)
$$ n=\frac{k_F^3}{3\pi^2} = \frac{3}{4\pi r_s^3}, $$we have (with \( k_F=1.92/r_s \),
$$ \frac{\varepsilon_{k}^{HF} }{\varepsilon_{0}^{F}}=x^2-\frac{e^2m}{\hbar^2 k_F\pi}F(x)=x^2-\frac{r_s}{a_0}0.663F(x), $$with \( r_s \sim 2-6 \) for most metals.
It can be convenient to use the Bohr radius \( a_0=\hbar^2/e^2m \). For most metals we have a relation \( r_s/a_0\sim 2-6 \).
c) Make a plot of the free electron energy and the Hartree-Fock energy and discuss the behavior around the Fermi surface. Extract also the Hartree-Fock band width \( \Delta\varepsilon^{HF} \) defined as
$$ \Delta\varepsilon^{HF}=\varepsilon_{k_{F}}^{HF}- \varepsilon_{0}^{HF}. $$Compare this results with the corresponding one for a free electron and comment your results. How large is the contribution due to the exchange term in the Hartree-Fock equation?
Solution.
We can now define the so-called band gap, that is the scatter between the maximal and the minimal value of the electrons in the conductance band of a metal (up to the Fermi level). For \( x=1 \) and \( r_s/a_0=4 \) we have
$$ \frac{\varepsilon_{k=k_F}^{HF} }{\varepsilon_{0}^{F}} = -0.326, $$and for \( x=0 \) we have
$$ \frac{\varepsilon_{k=0}^{HF} }{\varepsilon_{0}^{F}} = -2.652, $$which results in a gap at the Fermi level of
$$ \Delta \varepsilon^{HF} = \frac{\varepsilon_{k=k_F}^{HF} }{\varepsilon_{0}^{F}}-\frac{\varepsilon_{k=0}^{HF} }{\varepsilon_{0}^{F}} = 2.326. $$This quantity measures the deviation from the \( k=0 \) single-particle energy and the energy at the Fermi level. The general result is
$$ \Delta \varepsilon^{HF} = 1+\frac{r_s}{a_0}0.663. $$The following python code produces a plot of the electron energy for a free electron (only kinetic energy) and for the Hartree-Fock solution. We have chosen here a ratio \( r_s/a_0=4 \) and the equations are plotted as funtions of \( k/f_F \).
import numpy as np
from math import log
from matplotlib import pyplot as plt
from matplotlib import rc, rcParams
import matplotlib.units as units
import matplotlib.ticker as ticker
rc('text',usetex=True)
rc('font',**{'family':'serif','serif':['Hartree-Fock energy']})
font = {'family' : 'serif',
'color' : 'darkred',
'weight' : 'normal',
'size' : 16,
}
N = 100
x = np.linspace(0.0, 2.0,N)
F = 0.5+np.log(abs((1.0+x)/(1.0-x)))*(1.0-x*x)*0.25/x
y = x*x -4.0*0.663*F
plt.plot(x, y, 'b-')
plt.plot(x, x*x, 'r-')
plt.title(r'{\bf Hartree-Fock single-particle energy for electron gas}', fontsize=20)
plt.text(3, -40, r'Parameters: $r_s/a_0=4$', fontdict=font)
plt.xlabel(r'$k/k_F$',fontsize=20)
plt.ylabel(r'$\varepsilon_k^{HF}/\varepsilon_0^F$',fontsize=20)
# Tweak spacing to prevent clipping of ylabel
plt.subplots_adjust(left=0.15)
plt.savefig('hartreefockspelgas.pdf', format='pdf')
plt.show()
From the plot we notice that the exchange term increases considerably the band gap compared with the non-interacting gas of electrons.
We will now define a quantity called the effective mass. For \( \vert\mathbf{k}\vert \) near \( k_{F} \), we can Taylor expand the Hartree-Fock energy as
$$ \varepsilon_{k}^{HF}=\varepsilon_{k_{F}}^{HF}+ \left(\frac{\partial\varepsilon_{k}^{HF}}{\partial k}\right)_{k_{F}}(k-k_{F})+\dots $$If we compare the latter with the corresponding expressiyon for the non-interacting system
$$ \varepsilon_{k}^{(0)}=\frac{\hbar^{2}k^{2}_{F}}{2m}+ \frac{\hbar^{2}k_{F}}{m}\left(k-k_{F}\right)+\dots , $$we can define the so-called effective Hartree-Fock mass as
$$ m_{HF}^{*}\equiv\hbar^{2}k_{F}\left( \frac{\partial\varepsilon_{k}^{HF}} {\partial k}\right)_{k_{F}}^{-1} $$d) Compute \( m_{HF}^{*} \) and comment your results.
e) Show that the level density (the number of single-electron states per unit energy) can be written as
$$ n(\varepsilon)=\frac{Vk^{2}}{2\pi^{2}}\left( \frac{\partial\varepsilon}{\partial k}\right)^{-1} $$Calculate \( n(\varepsilon_{F}^{HF}) \) and comment the results.
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 \( \Omega \) 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{\Omega}}\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}. $$a) 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 \).
Solution.
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{_auto18} \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}{\Omega^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}{\Omega^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}{\Omega^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}{\Omega} \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}{2\Omega} \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}{2\Omega} \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} $$This equation will be useful for our nuclear matter calculations as well. 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 \).
b) Calculate thereafter the reference energy for the infinite electron gas in three dimensions using the above expressions for the kinetic energy and the potential energy.
Solution.
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{\Omega}{(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{\Omega}{(2\pi)^3} \left( 4\pi \int_0^{k_F} k^4 d\mathbf{k} \right) = \frac{4\pi\Omega}{(2\pi)^3} \frac{1}{5} k_F^5 = \frac{4\pi\Omega}{5(2\pi)^3} k_F^5 = \frac{\hbar^2 \Omega}{10\pi^2 m} k_F^5 . $$The density of states in momentum space is given by \( 2\Omega/(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{2\Omega}{(2\pi)^3} \frac{4}{3} \pi k_F^3 = \frac{\Omega}{3\pi^2} k_F^3 \quad \Rightarrow \quad k_F = \left( \frac{3\pi^2 N}{\Omega} \right)^{1/3}. $$This gives us
$$ \begin{equation} \langle \Phi_0\vert \hat{T} \vert \Phi_0\rangle = \frac{\hbar^2 \Omega}{10\pi^2 m} \left( \frac{3\pi^2 N}{\Omega} \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}{2\Omega} \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}{2\Omega} \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}{2\Omega} \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}{2\Omega} \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/\Omega \), with \( \Omega \) 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/\Omega \). 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}{\Omega}\frac{4\pi}{\mu^2}, $$and
$$ \hat{H}_{el-b}=-e^2\frac{N^2}{\Omega}\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.
c) Show thereafter 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}}{2\Omega}{\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}}. $$d) 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/\Omega \), \( 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.
e) Plot your results. Why is this system stable? Calculate thermodynamical quantities like the pressure, given by
$$ P=-\left(\frac{\partial E}{\partial \Omega}\right)_N, $$and the bulk modulus
$$ B=-\Omega\left(\frac{\partial P}{\partial \Omega}\right)_N, $$and comment your results.