Hartree method for Helium ground state

Leave a comment

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.

The radial part,

\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

The radial part

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)

Helium.PNG


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

 

Advertisements

Integrations of Laguerre polynomial

Leave a comment

Since the Laguerre polynomial is deeply connected to the hydrogen-like electron orbital. The radial solution is

\displaystyle R_{nl}(r) = A \frac{1}{r} \exp(- \frac{x}{2}) x^{l+1}L_{n-l-1}^{2l+1}(x)

\displaystyle x = \frac{2Z}{na_0} r

From the normalization condition, we can get the coefficient A

\displaystyle \int_0^{\infty} R_{nl}^2(r) r^2 dr

Beside of calculating the normalization factor, the general expectation value also need to evaluation of the integration of the Laguerre polynomial.

In here, we will show the calculation for 3 integrals:

\displaystyle \int_0^{\infty} x^{\alpha+m} e^{-x} L_n^\alpha(x) dx,  n,\alpha\geq0

\displaystyle \int_0^{\infty} x^{\alpha+k} e^{-x} (L_n^\alpha(x))^2 dx

\displaystyle \int_0^{\infty} x^{\alpha +k} e^{-x} L_m^{\beta}(x) L_n^\alpha(x) dx


Using the Rodrigues formula

\displaystyle L_n^\alpha(x) = \frac{1}{n!} x^{-\alpha} e^x \frac{d^n}{dx^n}(e^{-x}x^{n+\alpha})

\displaystyle \int_0^{\infty} x^{\alpha+m} e^{-x} L_n^\alpha(x) dx = \frac{1}{n!} \int_0^{\infty} x^{m} \frac{d^n}{dx^n}(e^{-x}x^{n+\alpha})dx

For m=0,

\displaystyle \int_0^{\infty} x^{\alpha+m} e^{-x} L_n^\alpha(x) dx \\ = \frac{1}{n!} \int_0^{\infty} \frac{d^n}{dx^n}(e^{-x}x^{n+\alpha})dx \\ = \frac{1}{n!} (F_{n-1}(\infty)-F_{n-1}(0) )

\displaystyle F_{n-1}(x) = \int_0^{\infty} \frac{d^n}{dx^n}(e^{-x}x^{n+\alpha})dx \\ = \sum_{i=0}^{n-1} (-1)^i \begin{pmatrix} n-1 \\ i \end{pmatrix} \frac{(n+\alpha)!}{(\alpha+1+i)!} x^{\alpha+1+i} e^{-x}

It is obvious that F(\infty) = 0 due to the e^{-x}. Since \alpha>0, then x^{\alpha+1+i} = 0 at x = 0. Thus,

\displaystyle \int_0^{\infty} x^{\alpha+m} e^{-x} L_n^\alpha(x) dx = 0 ,  m=0

Using integration by path,

\displaystyle \int_0^{\infty} x^{m} d \left(\frac{d^{n-1}}{dx^{n-1}}(e^{-x}x^{n+\alpha}) \right)dx \\ = [ x^{m} \frac{d^{n-1}}{dx^{n-1}}(e^{-x}x^{n+\alpha}) ]_0^{\infty} - (m) \int_0^{\infty} x^{m-1} \frac{d^{n-1}}{dx^{n-1}}(e^{-x}x^{n+\alpha})dx

For the same reason, the first term is zero for \alpha + m \geq 0. Repeat the integration by path k times,

\displaystyle \int_0^{\infty} x^{m} d \left(\frac{d^{n-1}}{dx^{n-1}}(e^{-x}x^{n+\alpha}) \right)dx \\ = (-1)^k \frac{(m)!}{(m-k)!} \int_0^{\infty} x^{m-k} \frac{d^{n-k}}{dx^{n-k}}(e^{-x}x^{n+\alpha})dx

When m \geq n || 0 > m, k can go to n, then,

\displaystyle \int_0^{\infty} x^{m} d \left(\frac{d^{n-1}}{dx^{n-1}}(e^{-x}x^{n+\alpha}) \right)dx \\ = (-1)^{n} \frac{(m)!}{(m-n)!} \int_0^{\infty} x^{m-n}(e^{-x}x^{n+\alpha})dx \\ = (-1)^{n} \frac{(m)!}{(m-n)!} (\alpha+m)!

When n > m > 0, k can only go to m. then,

\displaystyle \int_0^{\infty} x^{m} d \left(\frac{d^{n-1}}{dx^{n-1}}(e^{-x}x^{n+\alpha}) \right)dx \\ = (-1)^{m} (m)! \int_0^{\infty} \frac{d^{n-k}}{dx^{n-k}}(e^{-x}x^{n+\alpha})dx \\ = (-1)^{m} (m)!  \frac{d^{n-m-1}}{dx^{n-m-1}}(e^{-x}x^{n+\alpha}) = 0

For the similar result as F_{n-1}(x), the result is zero.

To summarize, the following formula suitable for all m > n || 0>m

\displaystyle \int_0^{\infty} x^{\alpha+m} e^{-x} L_n^\alpha(x) dx = (-1)^n (\alpha+m)! \begin{pmatrix} m \\ n \end{pmatrix}


To evaluate the integral,

\displaystyle \int_0^{\infty} x^{\alpha+k} e^{-x} (L_n^\alpha(x))^2 dx

We need the know that the Laguerre polynomial can be expressed as,

\displaystyle L_n^\alpha(x) = \sum_{i=0}^{n} (-1)^i \begin{pmatrix} n+\alpha \\ n-i \end{pmatrix} \frac{x^i}{i!}

Thus, we can evaluate

\displaystyle \int_0^{\infty} x^{\alpha+k} e^{-x} x^i L_n^\alpha(x) dx

Then combine everything,

\displaystyle \int_0^{\infty} x^{\alpha+k+i} e^{-x} L_n^\alpha(x) dx = (-1)^n (\alpha+i+k)! \begin{pmatrix} i+k \\ n \end{pmatrix}

put inside the series expression,

\displaystyle \int_0^{\infty} x^{\alpha+k} e^{-x} (L_n^\alpha(x))^2 dx = \sum_{i=0}^{n} (-1)^{i+n} \begin{pmatrix} n+\alpha \\ n-i \end{pmatrix} \begin{pmatrix} i+k \\ n \end{pmatrix} \frac{(\alpha+k+i)!}{i!}

In general, i + k \geq n || 0 > i+k .

For k = 0 , the sum only has 1 term for i = n

\displaystyle \int_0^{\infty} x^{\alpha} e^{-x} (L_n^\alpha(x))^2 dx = \frac{(\alpha+n)!}{n!}

Using the formula, for k=1, we have

\displaystyle \int_0^{\infty} x^{\alpha+1} e^{-x} (L_n^\alpha(x))^2 dx = \frac{(\alpha+n)!}{n!} (2n+\alpha +1)


To evaluate the expectation value of radial function

\displaystyle \langle R_{nl} | \frac{1}{r^m} | R_{nl} \rangle

We have to calculate the integral,

\displaystyle \int_0^\infty e^{-x} x^{2l+2-m} \left( L_{n-l-1}^{2l+1} (x) \right)^2 dx

For m = 0, we get the normalization constant,

\displaystyle A^2 \frac{na_0}{2Z} \int_0^{\infty} e^{-x} x^{2l+2} \left( L_{n-l-1}^{2l+1}(x) \right)^2 dx = 1

\displaystyle \langle \frac{1}{r} \rangle = \frac{1}{n^2} \frac{Z}{a_0}

\displaystyle \langle \frac{1}{r^2} \rangle = \frac{2}{n^3(2l+1)} \frac{Z^2}{a_0^2}

In general,

\displaystyle \langle \frac{1}{r^m} \rangle = \left(\frac{Z}{a_0}\right)^m \frac{2^{m-1}(n-l-1)!}{n^{m+1}(n+l)!} G_{nlm}

\displaystyle G_{nlm} = \sum_{i=0}^{n-l-1} (-1)^{n-l-1+i} \begin{pmatrix} n+l \\ n-l-1+i \end{pmatrix} \begin{pmatrix} 1-m+i \\ n-l-1 \end{pmatrix} \frac{(2l+2-m+i)!}{i!}


For the integral, for m \neq n

\displaystyle \int_0^{\infty} x^{\alpha +k} e^{-x} L_m^{\beta}(x) L_n^\alpha(x) dx

When k = 0, without lost of generality, when m < n , the integral is zero.

When k \neq 0, m<n

\displaystyle \int_0^{\infty} x^{\alpha +k} e^{-x} L_m^{\beta}(x) L_n^\alpha(x) dx \\ = \sum_{i=0}^{m} \frac{(-1)^i}{i!}\begin{pmatrix} m+\beta \\ m-i \end{pmatrix} \int_0^{\infty} x^{\alpha+k+i} e^{-x} L_n^\alpha(x) dx \\ = \sum_{i=0}^{m} \frac{(-1)^{i+n}}{i!}\begin{pmatrix} m+\beta \\ m-i \end{pmatrix} \begin{pmatrix} k+i \\ n \end{pmatrix} (\alpha+k+i)!

The sum only non-zero when m\geq i > n-k or -k > i \geq 0

This shows that the Laguerre polynomial is orthogonal with weighting e^{-x} x^{\alpha}.

 

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

n=2.PNG

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

ws.PNG

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}}

Molecular orbital

Leave a comment

Using 2 hydrogen atom, we started from 2 atomic orbital (AO),

\phi_1(\vec{r}) = \phi(\vec{r} - \vec{R_1})  

\phi_2(\vec{r}) = \phi(\vec{r} - \vec{R_2})  

The atomic orbital \phi is the ground state of hydrogen atom. From these 2 AO, we can construct linear combination of atomic orbitals (LCAO),

\psi_1(\vec{r}) = A_1( \phi_1(\vec{r}) + \phi_2(\vec{r}))

\psi_2(\vec{r}) = A_2( \phi_1(\vec{r}) - \phi_2(\vec{r}))

Notice that an electron can be in either \psi_1 or \psi_2, that means the electron is not belong to any atom, but belong to the molecule.

The last step is adding the spin, there are 2 spins \alpha, \beta, The molecular orbital (MO) are,

\chi_1(\vec{r}) = \psi_1(\vec{r}) \alpha

\chi_2(\vec{r}) = \psi_1(\vec{r}) \beta

\chi_3(\vec{r}) = \psi_2(\vec{r}) \alpha

\chi_4(\vec{r}) = \psi_2(\vec{r}) \beta

The Slater determinant can be formed by any 2 of them.

\Psi = \sqrt{2}^{-1} ( \chi_1(\vec{r_1})\chi_2(\vec{r_2}) - \chi_2(\vec{r_1})\chi_1(\vec{r_2}) ) = |\chi_1\chi_2 \rangle

We started by 1 atomic orbital, and there are and 2 atoms,  then 2 LCAOs are formed, and then multiplied by 2 spin states, gives 4 molecular orbitals. From these 4 MO and 2 electrons, we can from 6 different Slater determinant.


If we have different number of AO for different atom, we have AO

\vec{n} = (n_1, n_2,..., n_k)

The LCAO are

\displaystyle \psi_i = \sum_{r=1}^k C_{ir} \phi_r,  r = 1,2,..., n_k

The total number of LCAO are all permutation of the AOs. A simple construction is that only plus or minus sign are used. Then the number of LCAO ( not orthonormal ) is

\displaystyle  K = 2^{-1} \prod_{i=1}^{k} (2n_i) = 2^{k-1} \prod_{i=1}^{k} n_i

The MO are

\chi_{2i-1} = \psi_i \alpha \\ \chi_{2i} = \psi_i \beta , i = 1,2,..., K

Thus the number of Slater determinant are

\displaystyle N = {}_{2K}C_p = \frac{(2K)!}{(2K-p)!p!}


For example, we have 3 atoms, there are (1, 2, 2) AO, there are 16 LCAOs. Thus, we have 32 MO. With 2 electrons, there are 496 Slater determinant.

For C-H molecule, there are 6 electrons for carbon and 1 electron for hydrogen, use 6 AOs for carbon and 1 AO for hydrogen, the number of MO is 24 and 346,104 Slater determinant.

When using all Slater determinant, it is called full CI.

When only increase the AO, but limited number of Slater determinant, The energy is called Hartree limit.

A problem is to construct LCAOs for arbitrary AOs. For example, for (1,1,1) AOs, The number of possible combination with pule or minus sign is 4. However. I cannot construct orthonormal 4 LCAOs, because the system are over constrained (16 constrains with 12 unknown of C_{ij} )

 

 

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.

some work done

Leave a comment

recently, i am working on 3 things, both are calculation.

the 1st thing is a concrete foundation of microwave inducted optical nuclear polarization theory, which is the theory that my experiment is current using. i found some papers talk on the theoretical treatment on this subject. however, they are all start from some un-stated assumptions, for example, they are started with the Hamiltonian of hyperfine interaction, without state clear where the term comes from. i guess it is very elementary.  nevertheless, different papers may use different Hamiltonians. moreover, some of the mathematics techniques or steps were skipped and directly jump to result. Thus. that is my motivation on building a concrete theory, based on dipole -magnetic field interaction.

when i finished the hyperfine interaction, i have to apply a magnetic field and it cause a Zeemen splitting of energy level. and the Zeeman effect depends on the strength of the magnetic field. they are called strong field, intermedia field and weak field. in different field strength, the origina Hamiltonian may or may not be treated as a perturbation. Thus, i have to know the magnitude of the hyperfine structure. at the point, i come up with in-consistance value that cannot smoothy translate from weak field to intermedia field and strong field.

So, i restudy the atomic theory on Hydrogen atom. Start from Bohr model, fine structure to hyperfine structure and calculated the energy level numerically, so that i can get a full and complete picture on this matter. then i can now compare the magnetic field strength. it turns out that it is not that clear when is weak or strong, because it is depends on the total vale of the Hamiltonian, not the coefficient. Thus, for higher orbital angular momentum, a very small magnetic field can be treated as strong field, coz it already large enough to break the coupling. i am still working on the hyperfine Zeeman effect.

Meanwhile, my professor asked me to calculate the filling factor of the sample and NMR coil. the NMR signal strength depends on the geometry of the  coil and sample. and this geometric factor is called filling factor. it is a good name, i thing, since the maximum value is that the sample completely filled the coil. In order to calculated this, i parallel calculated the magnetic field generated by a cylindrical uniformly magnetized rod and the magnetic induction due to multiple from a single coil. the far field of the rod is done and basically, it is dipole, but the near field is complicated, the multiple terms appear and they are convergence slowly.

that’s why i don’t have post for 2 weeks.

Older Entries