Hartree method for 1D infinite potential well

Leave a comment

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


\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.


Reason for n>l in atomic orbit

Leave a comment

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.

Express Hydrogen orbit using fine structure constant

Leave a comment

The fine structure constant is

\displaystyle \alpha = \frac{e^2}{4\pi \epsilon_0} \frac{1}{\hbar c}

The Bohr radius

\displaystyle a_0 = \frac{4\pi \epsilon_0}{e^2} \frac{\hbar^2}{m} = \frac{\hbar}{m c \alpha}

The energy

\displaystyle E_n  = -\frac{1}{2n^2} \left(\frac{e^2}{2\pi \epsilon_0}\right)^2 \frac{m}{\hbar^2} = - \frac{1}{2} m (c\alpha)^2 \frac{1}{n^2}

We can remember that the electron is traveling at \alpha c at ground state.

When solving the radial part of the schrodinger equation, we can use

\displaystyle \rho = \frac{2r}{a_0} = r\sqrt{- \frac{8mE}{\hbar^2}}

\displaystyle \lambda = n = \sqrt{ -\frac{mc^2 \alpha^2}{2 E}}

A little summery of what I am doing

Leave a comment

My current position is developing a computational program that can measure the system of polarized target automatically and repeatedly. The program needs to connect with the microwave generator, voltage supply, power meter, and oscilloscope. That is purely technical and nothing special.

Later, after the program started to collecting data, another program is needed to analysis the data. Although there are many program that can do analysis, but those program are not so easy to use, in both control and display. So, I make an analysis program. The program is simply read the 2-D time-domain signal, and then determine some parameters of a specific function that fit the signal. However, the signal could be very noisy, so FFTW and wavelet analysis were implemented. That’s why the wavelet analysis appeared in this blog.

After that, the sample of polarized target is made from various pentacene derivatives, that no known energy level exist. One way to find out the energy levels of the singlet excited states is to measure the absorption spectrum. However, the energy of the triplet state is difficult to measure. And the energy level of the triplet state is critical for Dynamic Nuclear Polarization to be happened. Thus, one way to find out is doing computation chemistry.

My last chemistry class was like 15 years ago, when I was junior high school. But the basic of computation chemistry is solving the Schrodinger equation. That is what theoretical nuclear physicists do! The variation method, the Hartree-fock method, I heard it and somehow know it, but never do it with my own hand and computer. That is why I revisit the Hartree-fock method, and found out my previous understanding is so naive.

In the course of studying Hartree-Fock method, and one problem is evaluating the overlap integral

\displaystyle \int \psi_a(r_1) \psi_b(r_2) \frac{1}{r_{12}} \psi_c(r_1) \psi_d(r_2) dr_1dr_2

For a 3-D system, the integral involves product of multiple spherical harmonics. That is really troublesome. Therefore, I move to study the spherical harmonics, and the related rotational invariant, Wigner D-matrix, Clebcsh-Gordon series and Fourier series.  The spherical harmonics arises from solving the Laplace equation in spherical coordinate. A general theory of the solution of Lapalace equation involves Legendre polynomial, which is a special case for Hypergeometric function.  And a very interesting connection is that the elliptical function of the 1st and 2nd kind are also two special cases of hypergeometric function, that, the solution of Laplace equation for elliptical boundary condition is elliptical functions. That connects spherical harmonics and elliptical function! Wow!

I am now a bit off-track, that I am very interesting on the function of all functions. We know that there are many elementary functions, such as sin, cos, and Log. And even more kind of special functions, such as

  • Hermite — solving 1-D harmonic oscillator
  • Laguerre — the radial function of hydrogen
  • Legendre — the solution of the “\theta” of  Laplace equation
  • Gamma — a continuation of factorial
  • Bessel — solution of 3-D infinite square well
  • Elliptic — magnetic field of a solenoid
  • Dirac delta
  • Gaussian or Error function

For discrete argument

  • Clebsch-Gordon
  • Factoral
  • Binomial

As far as I know, the Hypergemetric function is like “mother of functions”, although not all special functions can be expressed as it. I am driven by curiosity, so, I am not sure where I will go. For instance, the transformation of Hypergeometric function is very interesting.

So, for now, as an ending, I found one article is very interesting. The History and Future of Special Functions, by Stephen Wolfram.

Variation method & Eigenvalue Problem

Leave a comment

The relationship I haven’t stated explicitly, so, here.

In solving the Schrodinger equation (or a general Hamiltonian problem )

H |\psi \rangle = E|\psi \rangle

The variation method is a golden method to find the minimum energy. For the minimum energy solution \psi_0

\displaystyle \frac{\langle \psi_0|H|\psi_0\rangle}{\langle \psi_0|\psi_0 \rangle} = E_0

For other solution

\displaystyle \frac{\langle \psi|H|\psi\rangle}{\langle \psi|\psi \rangle} \geq E_0

This implies the diagonal elements of H are always larger or equal the minimum energy.

If we expand the solution in to an orthogonal complete basis \phi_i

\displaystyle |\psi\rangle = \sum_i a_i |\phi_i\rangle

The Schrodinger equation

\displaystyle \sum_i a_i H|\phi_i\rangle = \sum_i a_i E |\phi_i\rangle

since the basis is orthogonal, multiple

\displaystyle \sum_i a_i \langle \phi_j |H|\phi_i \rangle = \sum_i a_i E \langle \phi_j |\phi_i \rangle

\displaystyle \sum_i \langle \phi_j | H|\phi_i\rangle a_i = E a_j

\displaystyle \sum_i H_{ji} a_i = E a_j

In matrix form

 H \cdot \vec{a} = E \vec{a}

this is the eignevalue problem.

The eigenvalue problem does not solve the ground state, but also the excited states.

Exchange energy

Leave a comment

In QM, when treating a system with identical particles, the exchange term raised because the system must be the same after particle exchanged.

I am reluctant to use the term “exchange interaction” or “exchange force”, “exchange potential”, because there is no interaction, no force carrier, no potential at all.

Suppose the Hamiltonian for single fermion is

H\psi = \epsilon \psi

The 2-particle system, the Schrodinger equation is

H_T \Psi = H_1 + H_2 + H'= E \Psi

where H' is the interaction between the 2 particles.

Now, guessing the total wave function to be

\Psi = a \psi_1(r_1) \psi_2(r_2) + b \psi_1(r_2)\psi_2(r_1) = a \Psi_n + b\Psi_e

The first term is “normal”, the second term is “exchanged”. substitute to the equation

 a H_T \Psi_n + b H_T \Psi_e = a E \Psi_n + b E \Psi_e

multiply both side with \Psi_1 or \Psi_2, and integrate r_1, r_2, we have

\begin{pmatrix} J & K \\ K & J \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} = E \begin{pmatrix} a \\ b \end{pmatrix}


J = \int \psi_1^*(r_1) \psi_2^*(r_2) H_T \psi_1(r_1) \psi_2(r_2) dr_1 dr_2

K = \int \psi_1^*(r_2) \psi_2^*(r_1) H_T \psi_1(r_1) \psi_2(r_2) dr_1 dr_2

Here, we assumed H' is symmetric against r_1, r_2. This is a fair assumption as H' is the mutual interaction, and H_1, H_2 are independent.

The J is the energy from the interaction of the particles.

The K is the energy due to EXCHANGE of the particles.

As we can see, the energy J = \epsilon_1 + \epsilon_2 + J'.

The H_T in K can be reduced to H' when the eigen wave function \psi_1 , \psi_2 are in difference orbits. When \psi_1 = \psi_2 , K = \epsilon_1 + \epsilon_2 + K'.

Thus the individual energy can be subtracted, yield,

\begin{pmatrix} J' & K' \\ K' & J' \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} = \Delta E \begin{pmatrix} a \\ b \end{pmatrix}

The solution is

\Delta E = J' + K', (a,b) = (1,1)

This is the symmetric state. For fermion, the spin-part is anti-symmetric. This is the spin-singlet state.

\Delta E = J' - K', (a,b) = (1,-1)

This is the anti-symmetric state. For fermion, the spin-part is symmetric, or spin-triplet state.

It is worth to notice that, in the symmetric state, the energy is higher because the two particles can be as close as possible, that create large energy. In the opposite, in the anti-symmetric state, the two particle never contact each other, thus, the interaction energy reduced.

We can summarized.

The energy J' is the shift of energy due to the mutual interaction.

But the system of identical particle subjects to the exchange symmetry. The exchanged wave function is also a state, that the total wave function has to include the exchanged wave function. This exchanged wave function creates the exchange term or exchange energy K'.

In the system of Fermion, the exchange energy is related to Pauli exclusion principle.

Hartree-Fock for 2-body ground state

Leave a comment

Long time ago, we talked about the mean-field calculation, an touched Hartree-Fock method. In that time, we explained excited-state approach. Now, we explain another approach by variation of the wave functions. This approach is inevitable in atomic physics, because the potential is fixed.

The Hamiltonian is

H = H_1 + H_2 + V_{12}

Since the spin component is anti-parallel, the space part of the total wave function is

\Psi = \phi_1(r_1) \phi_2(r_2)

The Schrodinger equation is

H\Psi = E\Psi

Integrate both side with \phi_1(r_1)

\int dr_1 \phi_1(r_1) H \Psi = \int dr_1 \phi_1(r_1) E \Psi

using the normalization and define \left<x| 1\right> = \phi_1(r_1)

\left<1|H|1\right> \phi_2(r_2) = E \phi_2(r_2)


\left<2|H|2\right> \phi_1(r_1) = E \phi_1 (r_1)

expand H = H_1 + H_2 + V_{12}

(H_2 + \left<1|V_{12}|1\right>) \phi_2(r_2) = (E - \left<1|H_1|1\right>) \phi_2(r_2)

(H_1 + \left<2|V_{12}|2\right>) \phi_1(r_1) = (E - \left<2|H_2|2\right>) \phi_1(r_1)

Now, we have 2 equations,  an initial guess of \Psi = \phi_1(r_1) \phi_2(r_2) ,

The difficulty is that, the \left<1|V_{12}|1\right> contains 2 variables.


Older Entries