Hartree-Fock method for Helium excited state

There are two excited wavefunctions, singlet and triplet state. In the singlet state, the total spin is zero, while in the triplet state, the total spin is 1.

In any case, the Slater determinant is

$\displaystyle \Psi(x,y) = \frac{1}{\sqrt{2}} \begin{vmatrix} \phi_1(x) & \phi_2(x) \\ \phi_1(y) & \phi_2(y) \end{vmatrix}$

For spin-singlet state, we set,

$\phi_1(q) = \uparrow \phi_{1s}(r), \phi_2(q) = \downarrow \phi_{2s}(r)$

For triplet state,

$\phi_1(q) = \uparrow \phi_{1s}(r), \phi_2(q) = \uparrow \phi_{2s}(r)$

Since the Hamitonian is spin independent, we only take care of the spatial part.

For the triplet state, if we put in the wavefucntion in the Hamiltonian. Using the basis set approximation, The Fock matrix is

$\displaystyle F_{ij} = H_{ij} + \sum_{hk}^M a_{\nu h}^*a_{\nu k} (G_{ij}^{hk} - \delta(\sigma_h, \sigma_k)G_{ik}^{hj})$

The factor $latex \delta(\sigma_h, \sigma_k)$ is due to the spin component.

From the beginning, the Fock matrix is

$\displaystyle \sum_{j}^M F_{ij} a_{\mu j} = \sum_{j}^M H_{ij} a_{\mu j} + \sum_{j}^M a_{\mu j} \sum_{hk}^M a_{\nu h}^*a_{\nu k} G_{ij}^{hk} - \sum_{j}^M a_{\nu j} \sum_{hk}^M a_{\nu h}^*a_{\mu k} G_{ij}^{hk}$

We can see the exchange term is difference. We can interchange $k \leftrightarrow j$

$\displaystyle \sum_{j}^M a_{\nu j} \sum_{hk}^M a_{\nu h}^*a_{\mu k} G_{ij}^{hk} = \sum_{k}^M a_{\nu k} \sum_{hj}^M a_{\nu h}^*a_{\mu j} G_{ik}^{hj}$

and pull the $a_{\mu j}$ out.

The calculation is straight forward. First, we calculate the 1s state due to the mean field created by the 2s orbital, then calculate the 2s state due to the mean field by the 1s orbital. And then repeat until converged.

Using the Hydrogen 1s, 2s, 3s, and 4s orbtial, we calculated

1s KE+PE = -2 Hartree = -54.42 eV
2s KE+PE = -0.4349 Hartree = -11.83 eV
The direct term = 0.2800 Hartree = 7.62 eV
The exchange term = 0.016 Hartree = 0.44 eV

The total energy = -2.171 Hartree = -59.07 eV.

Compare to the calculated ground state energy ( -77.11 eV), the excitation energy is 18.04 eV. And the ionization energy is 4.65 eV.

The experimental value for the triplet state excitation energy is 19.77 eV, and the ionization energy is 4.70 eV.

Consider the ground state energy is different from the experimental value by 2 eV. the calculation is fairly good.

For the singlet state, because the spin state of the two electron are opposite, the Fock matrix becomes Hartree matrix.

$\displaystyle F_{ij} = H_{ij} + \sum_{hk}^M a_{\nu h}^*a_{\nu k} G_{ij}^{hk}$

The calculated result is,

1s KE+PE = -2.000 Hartree = -54.42 eV
2s KE+PE = -0.400 Hartree = -10.89 eV
The direct term = 0.2531 Hartree = 6.89 eV
The exchange term = 0.040 Hartree = 1.09 eV

The total energy = -2.147 Hartree = -58.43 eV.

Compare to the calculated ground state energy ( -77.11 eV), the excitation energy is 18.68 eV. And the ionization energy is 4.00 eV.

The experimental value for the triplet state excitation energy is 20.55 eV, and the ionization energy is 3.92 eV.

It is interesting to compare with other formalism. Many other textbook use singlet wavefunction,

$\displaystyle \Psi_S(x,y) = \frac{1}{\sqrt{2}}( \phi_1(x) \phi_2(y) + \phi_2(x) \phi_1(y)) \frac{1}{\sqrt{2}}(\uparrow \downarrow - \downarrow \uparrow)$

The triplet state is the same,

$\displaystyle \Psi_T(x,y) = \frac{1}{\sqrt{2}}( \phi_1(x) \phi_2(y) - \phi_2(x) \phi_1(y)) \frac{1}{\sqrt{2}}(\uparrow \uparrow)$

Since the Hamiltonian is spin independent, we can see, when we calculate the total energy of the singlet state and triplet state are

$E_S = E_0 + J + K$

$E_T = E_0 + J - K$

Thus, the energy difference between singlet and triplet state is 2 exchange term.

In our calculation, the wave function of the singlet state is difference. We directly calculate from the Slater determinant. And the Fock matrix is lack of the exchange term due to the opposite spin. However, the triplet state is the same, the exchange term is 0.44 eV. If we do not do Hartree-Fock on the singlet state, but calculate the singlet state using the exchange term, the singlet state energy would be -58.19 eV. Comparing with the Hartree-Fock result (-58.43 eV), I would say they agree. The small difference is due to the difference of the mean fields.

Hartree method for Helium ground state

After long preparation, I am ready to do this problem.

The two electron in the helium ground state occupy same spacial orbital but difference spin. Thus, the total wavefunction is

$\displaystyle \Psi(x,y) = \frac{1}{\sqrt{2}}(\uparrow \downarrow - \downarrow \uparrow) \psi(x) \psi(y)$

Since the Coulomb potential is spin-independent, the Hartree-Fock method reduce to Hartree method. The Hartree operator is

$F(x) = H(x) + \langle \psi(y)|G(x,y) |\psi(y) \rangle$

where the single-particle Hamiltonian and mutual interaction are

$\displaystyle H(x) = -\frac{\hbar^2}{2m} \nabla^2 - \frac{Ze^2}{4\pi\epsilon_0 x} = -\frac{1}{2}\nabla^2 - \frac{Z}{x}$

$\displaystyle G(x,y) = \frac{e^2}{4\pi\epsilon_0|x-y|} = \frac{1}{|x-y|}$

In the last step, we use atomic unit, such that $\hbar = 1, m=1, e^2 = 4\pi\epsilon_0$. And the energy is in unit of Hartree, $1 \textrm{H} = 27.2114 \textrm{eV}$.

We are going to use Hydrogen-like orbital as a basis set.

$\displaystyle b_i(r) = R_{nl}(r)Y_{lm}(\Omega) \\= \sqrt{\frac{(n-l-1)!Z}{n^2(n+l)!}}e^{-\frac{Z}{n}r} \left( \frac{2Z}{n}r \right)^{l+1} L_{n-l-1}^{2l+1}\left( \frac{2Z}{n} r \right) \frac{1}{r} Y_{lm}(\Omega)$

I like the left the $1/r$, because in the integration $\int b^2 r^2 dr$, the $r^2$ can be cancelled. Also, the $i = nlm$ is a compact index of the orbital.

Using basis set expansion, we need to calculate the matrix elements of

$\displaystyle H_{ij}=\langle b_i(x) |H(x)|b_j(x)\rangle = -\delta \frac{Z^2}{2n^2}$

$\displaystyle G_{ij}^{hk} = \langle b_i(x) b_h(y) |G(x,y) |b_j(x) b_k(y) \rangle$

Now, we will concentrate on evaluate the mutual interaction integral.

Using the well-known expansion,

$\displaystyle G(x,y) = \frac{1}{|x-y|}=\frac{1}{r_{12}} = \sum_{l=0}^{\infty} \sum_{m=-l}^{l} \frac{4\pi}{2l+1} \frac{r_<^l}{r_>^{l+1}} Y_{lm}^{*}(\Omega_1)Y_{lm}(\Omega_2)$

The angular integral

$\displaystyle \langle Y_i(x) Y_h(y)| Y_{lm}^{*}(x) Y_{lm}(y) | Y_j(x) Y_k(y) \rangle \\ = \big( \int Y_i^{*}(x) Y_{lm}^{*}(x) Y_j(x) dx \big) \big( \int Y_h^{*}(y) Y_{lm}(y) Y_k(y) dy \big)$

where the integral $\int dx = \int_{0}^{\pi} \int_{0}^{2\pi} \sin(\theta_x) d\theta_x d\phi_x$.

From this post, the triplet integral of spherical harmonic is easy to compute.

$\displaystyle \int Y_h^{*}(y) Y_{lm}(y) Y_k(y) dy = \sqrt{\frac{(2l+1)(2l_k+1)}{4\pi (2l_h+1)}} C_{l0l_k0}^{l_h0} C_{lm l_km_k}^{l_hm_h}$

The Clebsch-Gordon coefficient imposed a restriction on $l,m$.

$\displaystyle \langle R_i(x) R_h(y)| \frac{r_<^l}{r_>^{l+1}} | R_j(x) R_k(y) \rangle \\ = \int_0^{\infty} \int_{0}^{\infty} R_i(x) R_h(y) \frac{r_<^l}{r_>^{l+1}} R_j(x) R_k(y) y^2 x^2 dy dx \\ = \int_0^{\infty} R_i(x) R_j(x) \\ \left( \int_{0}^{x} R_h(y) R_k(y) \frac{y^l}{x^{l+1}} y^2dy + \int_{x}^{\infty} R_h(x)R_k(x) \frac{x^l}{y^{l+1}} y^2 dy \right) x^2 dx$

The algebraic calculation of the integral is complicated, but after the restriction of $l$ from the Clebsch-Gordon coefficient, only few terms need to be calculated.

The general consideration is done. now, we use the first 2 even states as a basis set.

$\displaystyle b_{1s}(r) = R_{10}(r)Y_{00}(\Omega) = 2Z^{3/2}e^{-Zr}Y_{00}(\Omega)$

$\displaystyle b_{2s}(r) = R_{20}(r)Y_{00}(\Omega) = \frac{1}{\sqrt{8}}Z^{3/2}(2-Zr)e^{-Zr/2}Y_{00}(\Omega)$

These are both s-state orbital. Thus, the Clebsch-Gordon coefficient

$\displaystyle C_{lm l_k m_k}^{l_h m_h} = C_{lm00}^{00}$

The radial sum only has 1 term. And the mutual interaction becomes

$\displaystyle G(x,y) = \frac{1}{|x-y|}=\frac{1}{r_{12}} = 4\pi \frac{1}{r_>} Y_{00}^{*}(\Omega_1)Y_{00}(\Omega_2)$

The angular part

$\displaystyle \langle Y_i(x) Y_h(y)| Y_{lm}^{*}(x) Y_{lm}(y) | Y_j(x) Y_k(y) \rangle = \frac{1}{4\pi}$

Thus, the mutual interaction energy is

$G_{ij}^{hk} = \displaystyle \langle b_i(x) b_h(y) |G(x,y) |b_j(x) b_k(y) \rangle = \langle R_i(x) R_h(y)| \frac{1}{r_>} | R_j(x) R_k(y) \rangle$

$G_{ij}^{hk} = \displaystyle \langle R_i(x) R_h(y)| \frac{1}{r_>} | R_j(x) R_k(y) \rangle \\ \begin{pmatrix} G_{11}^{hk} & G_{12}^{hk} \\ G_{21}^{hk} & G_{22}^{hk} \end{pmatrix} = \begin{pmatrix} \begin{pmatrix} G_{11}^{11} & G_{11}^{12} \\ G_{11}^{21} & G_{11}^{22} \end{pmatrix} & \begin{pmatrix} G_{12}^{11} & G_{12}^{12} \\ G_{12}^{21} & G_{12}^{22} \end{pmatrix} \\ \begin{pmatrix} G_{21}^{11} & G_{21}^{12} \\ G_{21}^{21} & G_{21}^{22} \end{pmatrix} & \begin{pmatrix} G_{22}^{11} & G_{22}^{12} \\ G_{22}^{21} & G_{22}^{22} \end{pmatrix} \end{pmatrix} \\= \begin{pmatrix} \begin{pmatrix} 1.25 & 0.17871 \\ 0.17871 & 0.419753 \end{pmatrix} & \begin{pmatrix} 0.17871 & 0.0438957 \\ 0.0439857 & 0.0171633 \end{pmatrix} \\ \begin{pmatrix} 0.17871 & 0.0438957 \\ 0.0438957 & 0.0171633 \end{pmatrix} & \begin{pmatrix} 0.419753 & 0.0171633 \\ 0.0171633 & 0.300781 \end{pmatrix} \end{pmatrix}$

We can easy to see that $G_{ij}^{hk} = G_{ji}^{hk} = G_{ij}^{kh} = G_{hk}^{ij} = G_{ji}^{kh}$. Thus, if we flatten the matrix of matrix, it is Hermitian, or symmetric.

Now, we can start doing the Hartree method.

The general solution of the wave function is

$\psi(x) = a_1 b_{1s}(x) + a_2 b_{2s}(x)$

The Hartree matrix is

$F_{ij} = H_{ij} + \sum_{h,k} a_h a_k G_{ij}^{hk}$

The first trial wave function are the Hydrogen-like orbital,

$\psi^{(0)}(x) = b_{1s}(r)$

$F_{ij}^{(0)} = \begin{pmatrix} -2 & 0 \\ 0 & -0.5 \end{pmatrix} + \begin{pmatrix} 1.25 & 0.17871 \\ 0.17817 & 0.419753 \end{pmatrix}$

Solve for eigen system, we have the energy after 1st trial,

$\epsilon^{(1)} = -0.794702 , (a_1^{(1)}, a_2^{(1)}) = (-0.970112, 0.242659)$

After 13th trial,

$\epsilon^{(13)} = -0.880049 , (a_1^{(13)}, a_2^{(13)}) = (-0.981015, 0.193931)$

$F_{ij}^{(13)} = \begin{pmatrix} -2 & 0 \\ 0 & -0.5 \end{pmatrix} + \begin{pmatrix} 1.15078 & 0.155932 \\ 0.155932 & 0.408748 \end{pmatrix}$

Thus, the mixing of the 2s state is only 3.7%.

Since the eigen energy contains the 1-body energy and 2-body energy. So, the total energy for 2 electrons is

$E_2 = 2 * \epsilon^{(13)} - G = -2.82364 \textrm{H} = -76.835 \textrm{eV}$

In which ,

$G = \langle \psi(x) \psi(y) |G(x,y) |\psi(x) \psi(y) \rangle = 1.06354 \textrm{H} = 28.9403 \textrm{eV}$

So the energies for

From He to He++.  $E_2 = -2.82364 \textrm{H} = -76.835 \textrm{eV}$
From He+ to He++, $E_1^+ = -Z^2/2 = 2 \textrm{H} = -54.422 \textrm{eV}$.
From He to He+, is $E_1 = E_2 - E_1^+ = -0.823635 \textrm{H} = -22.4123 \textrm{eV}$

The experimental 1 electron ionization energy for Helium atom is

$E_1(exp) = -0.90357 \textrm{H} = -24.587 \textrm{eV}$
$E_1^+(exp) = -1.99982 \textrm{H} = -54.418 \textrm{eV}$
$E_2(exp) = -2.90339 \textrm{H} = -79.005 \textrm{eV}$

The difference with experimental value is 2.175 eV. The following plot shows the Coulomb potential, the screening due to the existence of the other electron, the resultant mean field, the energy, and $r \psi(x)$

Usually, the Hartree method will under estimate the energy, because it neglected the correlation, for example, pairing and spin dependence. In our calculation, the $E_2$ energy is under estimated.

From the $(a_1^{(13)}, a_2^{(13)}) = (-0.981015, 0.193931)$, we can see, the mutual interaction between 1s and 2s state is attractive. While the interaction between 1s-1s and 2s-2s states are repulsive. The repulsive can be easily understood. But I am not sure how to explain the attractive between 1s-2s state.

Since the mass correction and the fine structure correction is in order of $10^{-3} \textrm{eV}$, so the missing 0.2 eV should be due to something else, for example, the incomplete basis set.

If the basis set only contain the 1s orbit, the mutual interaction is 1.25 Hartree = 34.014 eV. Thus, the mixing reduce the interaction by 5.07 eV, just for 3.7% mixing

I included the 3s state,

$\epsilon^{(13)} = -0.888475 , (a_1^{(13)}, a_2^{(13)}, a_3^{(13)}) = (0.981096, -0.181995, -0.06579)$

The mutual energy is further reduced to 1.05415 Hartree = 28.6848 eV. The $E_2 = -77.038 \textrm{eV}$. If 4s orbital included, the $E_2 = -77.1058 \textrm{eV}$. We can expect, if more orbital in included, the $E_2$ will approach to $E_2(exp)$.

Hartree-Fock method for 1D infinite potential well

Following the previous post, I tested my understanding on the Hartree method. Now, I move to the Hartree-Fock method. The “mean field” energy of the Hartree-Fock is

$\displaystyle G_{\alpha \beta}= \sum_{j=i+1}^{N} \langle \psi_{\alpha}(i) \psi_{\nu} (j) | G(i,j) \left(| \phi_{\beta}(i) \phi_{\nu}(j) \rangle - |\phi_{\nu}(i) \phi_{\beta}(j) \rangle \right) \\ = \langle \alpha \nu | \beta \nu \rangle - \langle \alpha \nu | \nu \beta \rangle$

I also use the same method, in which, the trial wave function is replaced every iteration and the integration is calculated when needed.

Since the total wave function must be anti-symmetry under permutation, therefore one state can be occupied by only one particle. Thus, if we use the ground state in the mean field, the “meaningful” wave functions are the other states.

It is interesting that the mean field energy is zero when $\mu = \nu$, the consequence is no mean field for the same state. Suppose the mean field energy is constructed using the ground state, and we only use 3 states, the direct term is

$G_D = \begin{pmatrix} \langle 11|11 \rangle & \langle 11|21 \rangle & \langle 11|31 \rangle \\ \langle 21|11 \rangle & \langle 21|21 \rangle & \langle 21|31 \rangle \\ \langle 31|11 \rangle & \langle 31|21 \rangle & \langle 31|31 \rangle \end{pmatrix}$

The exchange term is

$G_E = \begin{pmatrix} \langle 11|11 \rangle & \langle 11|12 \rangle & \langle 11|13 \rangle \\ \langle 21|11 \rangle & \langle 21|12 \rangle & \langle 21|13 \rangle \\ \langle 31|11 \rangle & \langle 31|12 \rangle & \langle 31|13 \rangle \end{pmatrix}$

Due to the symmetry of the mutual interaction. We can see that some off-diagonal terms are cancelled. For example,

$\displaystyle \langle 1 1 | 3 1 \rangle = \int \psi_1^*(x) \psi_1^*(y) \cos(x-y) \psi_3(x) \psi_1(y) dy dx$

$\displaystyle \langle 1 1 | 1 3 \rangle = \int \psi_1^*(x) \psi_1^*(y) \cos(x-y) \psi_1(x) \psi_3(y) dy dx$

These two integrals are the same. In fact,

$latex \langle \alpha \nu | \beta \nu \rangle – \langle \alpha \nu | \nu \beta \rangle =$

whenever $\alpha = \nu$

$\displaystyle \langle \nu \nu | \beta \nu \rangle = \int \psi_\nu^*(x) \psi_\nu^*(y) \cos(x-y) \psi_\beta(x) \psi_\nu(y) dy dx$

$\displaystyle \langle \nu \nu | \nu \beta \rangle = \int \psi_\nu^*(x) \psi_\nu^*(y) \cos(x-y) \psi_\nu(x) \psi_\beta(y) dy dx$

We can see, when interchange $x \leftrightarrow y$, the direct term and the exchange term are identical, and then the mean field energy is zero. Also, when $\beta = \nu$ the mean field energy is also zero.

Due to the zero mean field, the off-diagonal terms of the Hamiltonian $H_{\alpha \beta}$ with $\alpha = \nu$ or $\beta = \nu$ are zero. Then, the eigen energy is the same as the diagonal term and the eigen vector is un-change.

Back to the case, the direct matrix at the 1st trial is,

$G_D = \begin{pmatrix} 0.720506 & 0 & -0.144101 \\ 0 & 0.576405 & 0 \\ -0.144101 & 0 & 0.555819 \end{pmatrix}$

The exchange matrix is

$G_E = \begin{pmatrix} 0.720506 & 0 & -0.144101 \\ 0 & 0.25 & 0 \\ -0.144101 & 0 & 0.0288202 \end{pmatrix}$

Thus, the Fock matrix is

$F = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 4.3264 & 0 \\ 0 & 0 & 9.527 \end{pmatrix}$

Therefore, the eigen states are the basis, unchanged

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

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

$\displaystyle \psi_3(x) = \sqrt{\frac{2}{\pi}} \sin(3x)$

Only the eigen energies are changed, as $\epsilon_1 = 1$$\epsilon_2 = 4.3264$$\epsilon_3 = 9.527$

The total wave function for 2 particles at state 1 and state-\mu is

$\Psi(x,y) = \frac{1}{\sqrt{2}} ( \psi_1(x) \psi_\mu(y) - \psi_\mu(x) \psi_1(y) )$

I found that the “mean field” function is not as trivial as in the Hartree case, because of the exchange term. In principle, the mean field for particle-i at state-$\mu$ is,

$G(i) = \int \phi_\nu^*(j) G(i,j) \phi_\nu(j) dj - \int \phi_\nu^*(j) G(i,j) \phi_\mu(j) dj$

However, the direct term is multiply with $\psi_\mu(i)$, but the exchange term is multiply with $\psi_\nu(i)$, which are two different functions. i.e., the “mean field” is affecting two functions, or the “mean field” is shared by two states.

Although the mean field for single state can be defined using exchange operator symbolically, but I don’t know how to really implement into a calculation. Thus, I don’t know how to cross-check the result.

Method on solving differential equation

The foundation of using Hartree-Fock method to get the self-consistence wave function and potential is solving the Hartree-Fock equation, which is kind of differential equation.

The idea of the Hartree-Fock method is

1. Assume the trial wave function $\psi^{(0)}_\mu(i)$
2. put the trial wave function into the Hartree-Fock equation $\displaystyle F(\psi_{\nu}^{(0)}(j))$
3. Solve the equation and obtain a new wave function $\psi^{(1)}_\mu(i)$
4. Go to step 2. until converge.

In this post, we will discuss at least 2 method on the step 3.

One method is solve the Hartree-Fock equation using numerical method, for example, Rungu-Kutta method, or Numerov’s method. Those method can also obtain the energy when impose some conditions, for example, the behaviour of the wave function at long range.

The 2nd method is basis expansion. We assume the wave function is a superposition of some basis,

$\displaystyle \psi_\mu(i) = \sum_{\alpha}^{N_\alpha} a_{\mu\alpha} \phi_\alpha(i)$

Substitute into the Hartree equation (Hartree-Fock equation is similar),

$\displaystyle F \left(\sum_{j=i+1}^{N}\psi_\nu(j) \right) \psi_\mu(i) = H(i) \psi_\mu(i) + \sum_{j=i+1}^{N}\langle \psi_\nu(j) | G(i,j) | \psi_\nu(j) \rangle \psi_\mu(i)$

$\displaystyle \sum_{\alpha} a_{\mu\alpha} F |\phi_{\alpha}(i) \rangle \\ = \sum_{\alpha} a_{\mu\alpha} \left( H(i) |\phi_{\alpha}(i) \rangle + \sum_{k,l} \sum_{j=i+1}^N a_{\nu k}^{*} a_{\nu l}\langle \phi_{\nu k}(j) | G(i,j) | \phi_{\nu l} (j) \rangle |\phi_{\mu \alpha}(i) \rangle \right) \\ = \sum_{\alpha} a_{\mu\alpha} \left( H(i) + G(i) \right) |\phi_{\alpha}(i) \rangle \\ = \sum_{\alpha} a_{\mu \alpha} E |\phi_{\alpha} \rangle$

Acting $\langle \phi_{\mu \beta}(i)|$

$\displaystyle \sum_{\alpha} a_{\mu \alpha} \langle \phi_{\mu \beta} (i) | F | \phi_{\mu \alpha} (i) \rangle = \sum_{\alpha} F_{\beta \alpha} a_{\mu \alpha} = E a_{\mu \beta}$

Thus, we can solve the coefficient $a_{\mu \alpha}$, which is the eigen vector. Then, new trial functions with the new coefficients are used.

In the 2nd method, we can see the integral

$H_{\beta \alpha} = \langle \phi_\beta(i) | H(i)| \phi_\alpha (i) \rangle$

$G_{\beta \alpha}^{k l} = \langle \phi_\beta(i) \phi_k(j) | G(i,j) | \phi_l(j) \phi_\alpha(i) \rangle$

can be pre-calculated. The calculation becomes iterating the coefficient $a_{\mu \alpha}$.

In method 2, the basis does not change, but only the coefficients. The advantage, also the disadvantage of this method is that, all possible integrals are pre calculated, even though some of them will never be used. Once the integrals are calculated, it can be used for solving many similar problems. But the calculation time for the integrals may be long and consume memory.

In previous post, I did not pre-calculate the integral, I simply put the new trial function and integrate again.

This method is very similar to the method 2. However, the basis is the wave function itself. After one calculation, a new basis, which is the wave function, will be used. And the old basis will be replaced. In method 2, only the coefficient is replaced.

This method may somehow faster then method 2, because only nesscary integrals are calculated. In Mathematica, method 2 is generally slower, especially for larger basis.

Most confusing thing for the method is that, weather the index is represent particle or state.

Look back the method. The total wave function is product of the wave function of each particle.

$\Psi = \phi_1(1) \phi_2(2) .. \phi_N(N)$

Here, we fixed our notation that the subscript is state, the parentheses is particle ID.

The single particle energy is,

$\displaystyle \langle \Psi |\sum_{i}^{N} H(i)|\Psi \rangle = \sum_{i}^{N} \langle \phi_\mu(i) | H(i)|\phi_\mu(i) \rangle$

Here, the $i$-particle in state $\mu$

The mutual interaction energy is

$\displaystyle \langle \Psi |\sum_{i}^{N-1}\sum_{j=i+1}^{N} G(i,j)|\Psi \rangle = \sum_{i}^{N-1}\sum_{j=i+1}^{N} \langle \phi_\mu(i) \phi_\nu(j) | G(i,j)|\phi_\mu(i) \phi_\nu(j) \rangle$

In both energy terms, we can factor out the sum of particle, we get the Hartree equation

$\displaystyle \left( H(i) + \sum_{j=i+1}^{N} \langle \phi_\nu(j) | G(i,j)| \phi_\nu(j) \rangle \right) |\phi_\mu(i) \rangle = \epsilon_\mu | \phi_\mu(i) \rangle$

Now, we can see, the sum is sum over the number of particle, not state.

When solving the Hartree equation using a given basis $b_k(x)$, we set the first trial of the solution to be

$| \phi_\mu(i) \rangle = | b_\mu(i) \rangle$

Then, in the first trial, the Hartree matrix is

$\displaystyle \hat{F}_{\alpha \beta} = \langle b_\alpha(i) | \left( H(i) + \sum_{j=i+1}^{N} \langle b_\nu(j) |G(i,j) | b_\nu(j) \rangle \right) | b_\beta(i) \rangle$

In here, we keep the particle ID. Inside the mean field integral, the trial state $b_\nu(j)$ is used for the particle-j.

After a trial, we will get the eigen vector for each state $v_\mu$, then we construct a new trial function and iterate until converge.

We can see, for multi-particle wave function, the sum is sum over the particle, even thought there are particles share the same state. In fact, since each particle should be at a state, the state should be a function of particle ID. So, the general Hartree matrix is

$\displaystyle \hat{F}_{\alpha \beta} = \langle \psi_\alpha(i)|\left( H(i) + \sum_{j=i+1}^{N} \langle \psi_{\nu(j)}(j) | G(i,j)| \psi_{\nu(j)}(j) \rangle \right) |\psi_\beta(i) \rangle$

where $\alpha, \beta = 1,2,3,..., k$, $k$ is number state to use.

In Hartree-Fock method, the sum is also over all other particle, not state.

Variational method on 1D potential well using plane wave basis

The potential is

$V(x) = \left \{ \begin{matrix} - V_0 & , -L/2 \leq x \leq L/2 \\ 0 & , else \end{matrix} \right \}$

The exact solution can be solved by matching boundary condition.

$\tan(\sqrt{\frac{2m}{\hbar^2}(E+V_0 ) \frac{L}{2}}) = \sqrt{\frac{|E|}{E+V_0}}$

The plane wave basis is

$b_n(x) =\frac{1}{\sqrt{a}} \exp\left( i \frac{2\pi}{a} n x \right)$

The matrix of the Hamiltonian is

$\displaystyle H_{ij} = \langle b_i | \frac{P^2}{2m} + V(x) | b_j \rangle$

$\displaystyle \langle b_i | \frac{P^2}{2m} | b_j \rangle = \delta_{ij} \frac{L}{a}\frac{\hbar^2}{2m} \left(\frac{2\pi}{a}\right)^2 i^2$

$\displaystyle \langle b_i | V(x) | b_j \rangle = -V_0 \frac{L}{a} \frac{\sin( \pi (i-j) L/a)}{\pi (i-j) L/a}$

Then, we solve for the eigen system.

I use mathematica for solving the eigen system.

For simplicity, I set $L = 2, V_0 = 1$. The exact solution is $E = -0.603898$.

The parameter $a$ controls the wavelength of the plane wave. By increasing the number of wave $n$, we effectively control the minimum wavelength $\lambda_0 = a/n$ Therefore, larger the $n$, the energy will approach to the actual value.

I generated $a$ from 10 to 200 with step 10, and $n$ from 5 to 60 with step 5. Here is the result. As we can see, the larger the $a$, we need more number of wave to go to the actual energy. When small $a$ and large $n$, the result quickly converge to the actual energy. For $n = 60, a = 40$, the energy is $E = -0.603792$, the difference is $\Delta E = 0.000106$ or 0.02% difference.

The following plot show $a = 40$

However, when the $a$ is too small, the plane wave cannot describe the actual wave function, then the converge fail. The following plot show $n = 60$

To investigate, I plot the case for $a = 3$ and $a = 10$.

We can see, since the well width is 2, if the maximum wavelength is only 3, then it cannot capture a longer wavelength component. And we can see the calculated wave functions are repeated.

Reason for n>l in atomic orbit

If we study atomic structure, we will find the principle quantum number must larger than azimuthal quantum number (the number of angular momentum), i.e.

$n > l$

But most text book give this result using mathematics that there is no solution for the radial equation when $n \leq l$. Some text book solves the radial equation using power series, and argue that the series has to be terminated at power of $n-1$ because the number of node is $n-1$. The exactly solution is related to Laguerre polynomial, and the polynomial is only defined for $n \geq l+1$.

Also, in nuclear physics, there is no restriction and $n$ and $l$ can take any integers. Why these two cases are different? of course, the key point is the potential difference. Coulomb potential is used in atomic orbit. Wood-Saxon potential (or finite square well ) is used in nuclear orbit.

The Schrodinger equation for potential $V(r)$ is

$\displaystyle \left( \frac{1}{r^2} \frac{\partial}{\partial r}\left( r^2 \frac{\partial}{\partial r}\right) - \frac{1}{r^2} L^2 - \frac{2m}{\hbar^2} V(r) \right) \psi(r, \Omega) = -\frac{2mE}{\hbar^2} \psi(r,\Omega)$

Separate the radial and spherical part $\psi(r, \Omega) = R(r) Y(\Omega)$

$\displaystyle L^2 Y(\Omega) = l(l+1) Y(\Omega)$

$\displaystyle \left( \frac{1}{r^2} \frac{d}{d r}\left( r^2 \frac{d}{d r}\right) - \frac{l(l+1)}{r^2} - \frac{2m}{\hbar^2} V(r) \right) R(r) = -\frac{2mE}{\hbar^2} R(r)$

The radial equation further reduce by using $u(r) = r R(r)$

$\displaystyle \frac{d^2 u}{dr^2} - \frac{2m}{\hbar^2} U(r) u(r) = \frac{2mE}{\hbar^2} u(r)$

we can see, the effective potential is

$\displaystyle U(r) = V(r) + \frac{\hbar^2}{2m} \frac{l(l+1)}{r^2}$

For Coulomb potential,

$\displaystyle U(r) = -\frac{e^2}{4\pi \epsilon_0 r} + \frac{\hbar^2}{2m} \frac{l(l+1)}{r^2}$

Introduce dimensionless quantities, or Express everything using the fine structure constant. $a_0$ is Bohr radius,

$\displaystyle \rho = \frac{2r}{a_0}, \lambda = \frac{e^2}{4\pi \epsilon_0} \sqrt{ - \frac{m}{2\hbar^2 E}} = n$

The equation becomes

$\displaystyle \frac{d^2u}{d\rho^2} + \left(\frac{\lambda}{\rho} - \frac{l(l+1)}{\rho^2}\right) u = \frac{1}{4} u$

The simplified effective potential,

$\displaystyle U(\rho) = \frac{n}{\rho} - \frac{l(l+1)}{\rho^2}$

Now, we related the principle quantum number and azimuthal quantum number in one simply formula.

Since the repulsive part due to the rotation is getting larger and larger, the minimum point of the effective potential is

$\displaystyle \rho_{min} = \frac{2l(l+1)}{n}$

When $l = n$, $\rho_{min} = 2(n+1) , r_{min} = (n+1)a_0$.

I plot the case for $n=2$

The bottom line is $l = 0$. We can see, when $l = 2$, the minimize point of the potential is at $r = 3 a_0$. The potential is vary shallow, and the bound state is very diffuses. Since the 2nd derivative of the wave function $u$ around the minimum point is

$\displaystyle \frac{d^2u}{d\rho^2} = \left( \frac{l(l+1)}{\rho^2} - \frac{n}{\rho} + \frac{1}{4} \right) u \approx \frac{1}{4} \frac{l(l+1) - n^2}{l(l+1)} u = k u$

For $l = n-1$ , $k$ is negative, that give oscillating wave function.

For $l = n$, $k$ is positive, that gives explosive wave function that cannot be normalized.

Therefore, $l < n$.

In fact, if we plot the effective potential,

$\displaystyle U(r) = -\frac{e^2}{4\pi \epsilon_0 r} + \frac{\hbar^2}{2m} \frac{l(l+1)}{r^2}$

we can find that when $l \geq 1$, the potential is also very shallow. For such a picture, we can see the wave function must be oscillating that the number of node must be larger then 1, so the principle quantum number should be larger than 1.

In nuclear orbit, the effective potential using Wood-Saxon potential is

$\displaystyle U(r) = -\frac{1}{1-\exp(\frac{r-R}{a})} + \frac{\hbar^2}{2m} \frac{l(l+1)}{r^2}$

where $R \approx 1.25 A^{1/3} [fm]$ and $a = 0.6 [fm]$. For simplicity, we use

$\displaystyle U(r) = -\frac{1}{1-\exp(\frac{r-R}{a})} + \frac{l(l+1)}{r^2}$

The position of minimum potential is no simple formula. I plot R=10,  a=1

We can see, for any $l < l_{u}$ there is a “well” that can bound a nucleon.

I haven’t solve the finite spherical well, but the infinite spherical well is solved. The solution is the spherical Bessel function of first kind. that function only has 1 parameter, which is the angular momentum. To determine the energy is simply find the node position, and the number of node is related to the principle quantum number. Therefore, the principle quantum number and angular momentum is unrelated. A similar result should be applied well in the case of finite spherical well like Wood-Saxon type, that, for a fixed angular momentum, there can be many solution with different number of nodes.

Wood-Saxon potential, the effective potential is always sloping, especially when outside of the well, that forces that wave function to be decay faster. However, the Coulomb potential is a long range force, the effective potential can be slowly changing for long distance, that allow the wave function to be oscillating for long range.