In order to test my understanding on the Hartree method. I played a simple case. The potential is

$V(x) = \left\{ \begin{matrix} 0 , 0 \leq x \leq \pi \\ \infty , else \end{matrix} \right\}$

The mutual interaction between 2 particles is, for simple calculation,

$G(x_1, x_2) = g_0 \cos(x_1 - x_2)$

This mutual interaction is repulsive, we can expect the eigen wave function will be more spread out.

We use the eigen state of the $V(x)$ potential to be the basis,

$\displaystyle \phi_n(x) = \sqrt{\frac{2}{\pi}} \sin(nx)$

The Hartree method is basically solving the Hartree equation

$\displaystyle \left( -\frac{\hbar^2}{2m} \frac{d^2}{dx_1^2} + V(x_1) + \int G(x_1, x_2) (\psi_k(x_2))^2 dx_2 \right) \psi_n(x_1) = \epsilon_n \psi_n(x_1)$

We assume the two particles are in the same ground state, thus, $k = n = 1$. We assume,

$\psi_1(x) = \sum a_i \phi_i(x)$

The algorithm is that, for the 1 particles matrix

$\displaystyle F_{ij} = \int \psi_i(x) \left( -\frac{\hbar^2}{2m} \frac{d^2}{dx^2} + V(x) \right) \psi_j(x) dx$

Define

$\displaystyle G_0(x) = \int G(x, x_2) (\psi_1(x_2))^2 dx_2$

The 2-particle matrix

$\displaystyle G_{ij} = \int \psi_i(x) G_0(x) \psi_j(x) dx$

For the 1st trial, we set

$\psi_i^{(0)}(x) = \phi_{(2i+1)}(x)$

We only use “even” function with respect to $x = \pi/2$, because mixing parity will result non-reasonable result. ( why ? )

Then we get the matrix of the Hartree Hamiltonian $H_{ij} = F_{ij} + G_{ij}$, and solve for the eigen value and eigen vector,

$\vec{v}^{(0)} = (v_1^{(0)}, v_2^{(0)}, .... v_n^{(0)})$

The next trial is

$\displaystyle \psi_i^{(1)}(x) = \sum_i v_i^{(0)} \psi_i^{(0)}$

The process repeated until converge. The convergent can be exam using the Hartree Hamiltonian that the non-diagonal terms are vanished, i.e.

$F_{ij} + G_{ij} = 0, i\neq j$

I tried with 2 basis and 6 basis, I set $\hbar^2/2m = 1$ for simplicity, the Hartree ground state is 1.7125 for both basis. The eigen wavefunction for the 2-basis is

$\psi(x) = 0.797752 \sin(x) + 0.0145564 \sin(3x)$

with the mean field

$V_m(x) = 0.842569 \sin(x)$

For 6 basis, the wave function is

$\psi(x) = 0.797751 \sin(x) + 0.0145729 \sin(3 x) + 0.000789143 \sin(5 x) + 0.000125871 \sin(7 x) + 0.0000336832 \sin(9 x) + 0.0000119719 \sin(11 x)$

we can see, the contribution of the higher energy states are getting smaller and smaller. The higher energy states only contributes 0.03%. Following is the plot compare $\phi_1(x)$ (orange) and $\psi(x)$ (blue). We can see that the wave function becomes little spread out.

If we increase the strength of the mutual interaction by 10 times, the wave function becomes more spread out.

To cross check the result,  we can use the mean field and solve the Schrodinger equation using variational method. The method gives the ground state energy be 1.7125 and converge.

This is interesting that the variational method needs 4 or 5 basis to converge, while Hartree method only need 2.

At last, we substitute the wave functon back to the Hartree equation, it turns out the result is not satisfied 100%. The total wave function is

$\Psi(x_1, x_2) = \psi(x_1) \psi(x_2)$

The Schrodinger equation is

$\displaystyle H \Psi = -\frac{\hbar^2}{2m}\left(\frac{d^2}{dx_1^2} + \frac{d^2}{dx_2^2}\right) \Psi(x_1,x_2) + G(x_1,x_2)\Psi(x_1,x_2) = E \Psi(x_1,x_2)$

However, the left side and the right side are not equal, but the integration

$\langle \Psi |H|\Psi\rangle = \langle |\Psi|\Psi \rangle$

is equal.

I plot the $\Psi H\Psi$ (orange) and $E \Psi^2$ (blue), and there difference is the following

We can imagine, the integration of the bottom plot is zero. I don’t know the reason, but I guess it is not because of the number of basis, because the higher energy states only contributes for 0.03%. I guess, it is because of the energy is the quantity that is optimized rather then the wave function. And I am not sure that the reverse of the following is true.

$H|\Psi\rangle = E |\Psi\rangle \rightarrow \langle \Psi|H|\Psi \rangle = E \langle \Psi|\Psi \rangle$

For the case of many particles, the $G_0$ has to sum over all other particle. When all particles are in the same states, the sum is simple. For particles distribute among different states, in general, the sum is

$\displaystyle G_0 =\sum_{k\neq i} \int G(x_1,x_2) \psi_k(x_2)^2 dx$

Here, the sum is excluded the particle $x_1$, which is in state $i$. For example, if the system consists of 2 particles on state-1, and 2 particles on state-2. The mean field for the particle in state-1 is constructed using 1 particles on state-1 and 2 particles on state-2. Thus the energy for the state-2 may be not correct.

The mean field for the particle on state-2, the mean field should be constructed by 2 particles on state-1 and 1 particle on state-2.

For higher states that without particle, the mean field is built by 2 particles on state-1 and 2 particles on state-2.

In general, the mean field is constructed using other particles.

This fixed mean field method is simple but we have to define the mean field for each state and calculate again and again.