Laguerre | Nuclear Physics 101

Notes on the original Nilsson paper

June 21, 2020

GoLuckyRyan Basic Laguerre, Nilsson orbital, paper reading, Wigner-D Leave a comment

The paper title is “Binding states of individual nucleons in strongly deformed nuclei” by Sven Gösta Nilsson on 1955.

In the introduction, the total nuclear wave function $\Psi$

$\Psi = \phi_\Omega D_{MK}^I$

where $D_{MK}^I$ is Wigner D-matrix for rotation, $\phi_\Omega$ is the intrinsic motion of all nucleons. The vibration component is skipped. Because of imcompressible nuclear matter, a vibration needs a lot of energy.

The total wave function should preserve parity, so, a complete wave function should be

$\displaystyle \Psi = \sqrt{\frac{2I+1}{16\pi^2}} \left( \phi_\Omega D_{MK}^I + (-1)^{I+K} \phi_{-\Omega} D_{M-K}^I \right)$

Below is a picture of the quantum number.

Annotation 2020-06-21 134519

In large deformation, $J$ is not a good quantum number, but $\Omega$ is always a good one.
At ground state, $K = \Omega$ , so, the rotation angular momentum is perpendicular to the body frame axis.

The single nucleon potential is the usual.

$\displaystyle \frac{H}{\hbar \omega_0} = -\frac{1}{2} \nabla^2 + \frac{1}{2}r^2 - \beta r^2 Y_{20} + C l\cdot s + D l^2$

The basis in Nilsson’s paper is the eigen state of harmonic oscillator in L-S representation. The Hamiltonian is

$\displaystyle \frac{H_0}{\hbar \omega_0} = -\frac{1}{2} \nabla^2 + \frac{1}{2}r^2$

with

$\displaystyle H_0 \left| N l m_l m_s \right> = \hbar \omega_0 \left( N + \frac{3}{2} \right)\left| N l m_l m_s \right>$

The original paper use $\Sigma = m_s , \Lambda = m_l$ . And the function solution is

$\displaystyle \left<r | N l m_l m_s \right> = A r^l e^{-\frac{r^2}{2}} F\left(- \frac{N-l}{2}, l+ \frac{3}{2}, r^2 \right) Y_{lm_l} \chi_{s m_s}$

where $F(a,b,x)$ is the confluent hypergeometric function, and it can be expressed as Laguerre polynomial. So, the solution can be rewrite as

$\displaystyle \left<r | N l m_l m_s \right> = A r^l e^{-\frac{r^2}{2}} L_{\frac{N-l}{2}}^{l+\frac{1}{2}}(r^2) Y_{lm_l} \chi_{s m_s}$

In contrast, in my calculation, I use the $J-\Omega$ presentation, that, the connection is

$\displaystyle \left| N l j \Omega \right> = \sum C_{l m_l s m_s}^{j \Omega } \left| N l m_l m_s \right>$

The different is that, Nilsson needs additional transformation to calculate the $C_{jl}$ coefficient, and the calculation of the $\left< l\cdot s \right>$ and $\left< l^2 \right>$ is a bit complicated, due to L-S is not a good quantum number when $l\cdot s$ interaction was included. Thus, in Nilsson paper, he spent sometimes to talk about the calculation of $\left< l\cdot s \right>$ and $\left< l^2 \right>$ .

Next, Nilsson gives the calculation parameters of $\kappa, \mu$ . And since he is using L-S representation, the Nilsson orbital is expressed in that basis. Here is a comparison between my calculation and Nilsson calculation.

Annotation 2020-06-21 141037

Next, he explained that for $\beta > 0$ , energy increase with increase $\Omega$ due to “surface coupling”. It can be imagine like this:

Annotation 2020-06-21 141313

In above picture, when $\Omega$ is smallest, $l$ is perpendicular to the body axis, so the nucleon has large overlap with the whole nucleus, thus it is most bounded.

For small deformation, the deformation field $r^2 Y_{20}$ is treated as a perturbation, so that $j$ is a good quantum number.

For large deformation, the spin-orbital interaction $l\cdot s, l^2$ is treated as perturbation, and the good quantum number is $N, n_z, m_l , m_s$ , since $\Omega = m_s + m_l$ , so $\Omega$ is also a good quantum number.

In our previous notation, Nilsson orbital is notated as $K[N n_z \lambda]$ or $[N n_z \lambda] K$ , which is equivalent in Nilsson as $[N n_z m_l] \Omega$

The wave function of many nucleons.

After established the single nucleon wave function. The receipt of the construction of many nucleons wave function is

select a set of Nilsson orbitals and form the Slater determent.
Minimum the total energy by adjusting the deformation parameter.

$\displaystyle H_t = \sum_i T_i + \frac{1}{2} \sum_{i \neq j} V_{ij}$

It is interesting that the total wave function is not a mixture of various Slater determents from different combination of Nilsson orbitals, but rather a single Slater determent.

Ground state spin

Since each Nilsson orbital is degenerated to $\pm \Omega$ , which are rotate oppositely. For even-even nucleus, the ground state spin must be zero. For even-odd nucleus, the ground state spin is equal to the $\Omega$ of the last single nucleon.

For odd-odd nucleus, the ground state spin can be $\Omega = |\Omega_p \pm \Omega_n|$ , the p-n interaction decide which $\Omega$ is the ground state.

Decoupling parameter

For odd-A nuclei, the rotational energy is modified by a decoupling factor

$\displaystyle E_r = \frac{\hbar^2}{2 I_M} \left( I(I+1) + a (-1)^{I+\frac{1}{2}} \left(I+\frac{1}{2}\right) \right)$

And

$\displaystyle a = \sum_{j} (-1)^{j - \frac{1}{2}} \left( j+ \frac{1}{2} \right) |C_{jl}^2|$

Magnetic moments, EM transition probability, and ft-value for beta-decay are skipped.

Besides of the skipped material, it turns out the original Nilsson paper did not surprise me. And I still don’t really understand the “decoupling” thing.

Diagonal elements for Nilsson orbital in spherical-spin function

March 30, 2020

GoLuckyRyan Basic Laguerre, Nilsson orbital Leave a comment

In this post, we explain how to calculate the Nilsson orbital using perturbation method by compute the matrix elements using spherical-spin function. In that post, I said I will give the calculation for the perturbation element. Here we go for the diagonal elements

The diagonal matrix element,

$\displaystyle \left<Nljk|r^2 Y_{20}|Nljk\right>$

where the spherical-spin wave-function is

$\displaystyle |Nljk\rangle = A r^l e^{-\frac{r^2}{2}} L_{\frac{N-l}{2}}^{l + \frac{1}{2}}(r^2) \sum_{m m_s} Y_{lm}(\theta, \phi) \chi_{\frac{1}{2} m_s} C_{lm\frac{1}{2} m_s}^{jk}$

$\displaystyle A = \sqrt{\frac{(\frac{N-l}{2})!(\frac{N+l}{2})! 2^{N+l+2}}{\sqrt{\pi} (N+l+1)!}}$

First, the radial part is

$\displaystyle R_{Nl}(r) = r^l e^{-\frac{r^2}{2}} L_{\frac{N-l}{2}}^{l + \frac{1}{2}}(r^2)$

The radial integral is easy, we can use the integration formula. The radial integration is,

$\displaystyle \int_0^{\infty} r^2 R_{Nl}^2(r) r^2 dr = \int_0^{\infty} r^{2l+4} e^{-r^2} \left(L_{\frac{N-l}{2}}^{l + \frac{1}{2}}(r^2)\right)^2 dr$

When using the integration formula, one has to be careful when changing variable $r^2 \rightarrow r$ . Since $dr^2 = 2 r dr$ , we have to pull $2r$ out to properly do the change of variable.

$\displaystyle \int_0^{\infty} R_{Nl}^2(r) r^2 dr = \frac{1}{2}\int_0^{\infty} r^{2(l+1/2 + 1)} e^{-r^2} \left(L_{\frac{N-l}{2}}^{l + \frac{1}{2}}(r^2)\right)^2 (2 r dr)$

Set $r^2 \rightarrow x$ , $l + 1/2 \rightarrow \alpha$ , and $n = \frac{N-l}{2}$

$\displaystyle = \frac{1}{2}\int_0^{\infty} x^{\alpha+1} e^{-x} \left(L_{N}^{\alpha}(x)\right)^2 (dx) = \frac{1}{2}\frac{(\alpha+n)!}{n!}(2n+\alpha+1)$

Thus,

$\displaystyle \int_0^{\infty} R_{Nl}^2(r) r^2 dr = \frac{1}{2}\frac{( \frac{N+l+1}{2})!}{(\frac{N-l}{2})!}(N+\frac{3}{2})$

The angular-spin part is

$\displaystyle S_{ljk}(\theta,\phi, m_s) = \sum_{m m_s} Y_{lm}(\theta, \phi) \chi_{\frac{1}{2} m_s} C_{lm\frac{1}{2} m_s}^{jk}$

This contains spatial and spin part.

$\displaystyle S_{ljk}^*(\theta,\phi, m_s) Y_{20} S_{ljk}(\theta,\phi, m_s) \\=\sum_{m' m'_s} \sum_{m m_s} Y_{lm'}^*(\theta,\phi) Y_{20} Y_{lm}(\theta, \phi) \chi_{\frac{1}{2} m'_s} \cdot \chi_{\frac{1}{2} m_s} C_{lm'\frac{1}{2} m'_s}^{jk} C_{lm\frac{1}{2} m_s}^{jk}$

The dot-product of the spin part restricted the $m'_s$

$\displaystyle = \sum_{m' m m_s} Y_{lm'}^*(\theta,\phi) Y_{20} Y_{lm}(\theta, \phi) C_{lm'\frac{1}{2} m_s}^{jk} C_{lm\frac{1}{2} m_s}^{jk}$

And since $m+m_s = k$ and $m'+m_s =k$ , therefore $m = m'$ ,

$\displaystyle = \sum_{m m_s} Y_{lm}^*(\theta,\phi) Y_{20} Y_{lm}(\theta, \phi) \left(C_{lm\frac{1}{2} m_s}^{jk}\right)^2$

And the integration of spherical harmonic gives,

$\displaystyle \int Y_{lm}^*(\theta,\phi) Y_{20} Y_{lm}(\theta, \phi) = \sqrt{\frac{5}{4\pi}} C_{20l0}^{l0} C_{20lm}^{lm}$

$\displaystyle C_{20l0}^{l0} = -\frac{l+1}{\sqrt{4l^3+8l^2+l-3}} , l > 0$

$\displaystyle C_{20l0}^{l0} = -\frac{l(l+1)-3m^2}{\sqrt{l(4l^3+8l^2+l-3)}} , l > 0$

Sum up everything,

$\displaystyle \left<Nljk|r^2 Y_{20}|Nljk\right> \\ = \frac{(\frac{N-l}{2})!(\frac{N+l}{2})! 2^{N+l+2}}{\sqrt{\pi} (N+l+1)!} \frac{1}{2}\frac{( \frac{N+l+1}{2})!}{(\frac{N-l}{2})!}(N+\frac{3}{2}) \sum_{m m_s} \sqrt{\frac{5}{4\pi}} C_{20l0}^{l0} C_{20lm}^{lm} \left(C_{lm\frac{1}{2} m_s}^{jk}\right)^2 \\ = \sqrt{\frac{5}{4}} \frac{(\frac{N+l}{2})! ( \frac{N+l+1}{2})! 2^{N+l+1}}{\pi (N+l+1)!}(N+\frac{3}{2}) \sum_{m m_s} C_{20l0}^{l0} C_{20lm}^{lm} \left(C_{lm\frac{1}{2} m_s}^{jk}\right)^2$

Note: I haven’t numerically check the formula. ( may be later )

For the off-diagonal element. The angular-spin part should be similar. The difficulty is the radial part, we have to evaluate the most general orthogonal relation of the Laguerre polynomial with weighting $x^{\frac{1}{2}(\alpha+\alpha'+2)} e^{-x}$ .

We only knew the Laguerre polynomil is orthogonal with respect to $r^\alpha e^{-r}$ , i.e.

$\displaystyle \int_0^{\infty} L_n^\alpha L_m^\alpha r^\alpha e^{-r} dr = \frac{(\alpha+n)!}{n!} \delta_{nm}$

But not this.

$\displaystyle \int_0^{\infty} L_n^\alpha L_m^\beta r^\frac{\alpha+\beta}{2} e^{-r} dr$

or this

$\displaystyle \int_0^{\infty} L_n^\alpha L_m^\beta r^{(\frac{\alpha+\beta}{2}+k)} e^{-r} dr$

But we are quite sure the last one, with $\alpha \neq \beta$ will not give zero, otherwise, The Nilsson orbital will be very simple and boring.

Axial Harmonic Oscillator – Nilsson Orbit (III)

January 26, 2018

GoLuckyRyan Basic Laguerre, Nilsson orbital 1 Comment

One of the problem (or difficulty) is the conversion between cylindrical quantum number $|Nn_zm\rangle$ and spherical quantum number $|Nlm\rangle$ . In the limit of $\delta = 0$ , all state with same $N$ has same energy and the conversion is not necessary.

With out LS coupling, the energy level can be found by the approximation

$\displaystyle E(N, n_z) = \hbar\omega \left( \frac{3}{2} + \left(1 + \frac{\delta}{3} \right) N - \delta n_z \right)$

We can see from the formula, for $\delta>0$ , the largest $n_z$ has lowest energy. And the conversion to spherical harmonics can be done using projection method.

With LS coupling, the conversion becomes complicated. One issue is the coupling with the spin. Since the total angular momentum $J = L + S$ is not a conservative quantity due to deformation. In the body-fixed frame, the additional good quantum number is the z-projection of $J$ , which denote as $K = |m_l \pm 1/2 |$ , the quantum state becomes

$|Nn_z m_l K\rangle$

Note that each state will have only 2 degeneracy with negative $K$ . On the spherical basis, the coupling with $J$ is a standard textbook content. the eigen state is

$|N l j mj \rangle = |N l j K\rangle$

The conversion get complicated, because there is no clear rule on how $| Nn_z m_l K\rangle \rightarrow |Nl j m_j \rangle$ . When we want to calculate the energy level using cylindrical basis with the LS coupling, since the $J$ is not clear, even through we can write down the formula

$\displaystyle E( N, n_z, l, j) = \hbar \omega_z \left(\frac{1}{2} + n_z \right) + \hbar \omega_\rho \left(1 + N - n_z \right) \\ + a \frac{1}{2} ( j(j+1) - l(l+1) - s(s+1)) + b l^2$

It is difficult to know the $j , l$ .

The go around this conversion, we can diagonal the Hamiltonian using spherical basis. The Hamiltonian is can be written as

$\displaystyle H = - \frac{\hbar^2}{2m} \nabla^2 + \frac{m\omega}{2}r^2 \left( 1 - \frac{2}{3}\sqrt{\frac{16\pi}{5}} Y_{20}(\theta, \phi) \delta \right) + a L\cdot S + b L^2$

The first 2 terms is the spherical harmonic oscillator. The energy is $\hbar\omega( 3/2 + N)$ . The $Y_20$ restricts the coupling between $l, l \pm 2$ (see here.)

For example, when using the basis on $N = 1$ , which are

$\displaystyle |N l j m_j \rangle = \left( |1 1 \frac{3}{2} \frac{3}{2}\rangle , |1 1 \frac{1}{2} \frac{1}{2}\rangle, |1 1 \frac{3}{2} \frac{1}{2}\rangle \right)$

at $\delta = 0.5$ , the matrix is

$\begin{pmatrix} 2.633 & 0 & 0 \\ 0 & 2.233 & 0 \\ 0 & 0 & 2.233\end{pmatrix}$

The eigen values and eigen state are

$\displaystyle 2.633 , |1 1 \frac{3}{2} \frac{3}{2}\rangle$

$\displaystyle 2.233 , |1 1 \frac{1}{2} \frac{1}{2}\rangle , |1 1 \frac{3}{2} \frac{1}{2}\rangle$

The lower level has mixed $j$ !

Of course, a more complete diagonalization should be admixture from other major shell as well. When we disable the LS coupling, using 7 major shells, the Nilsson diagram look like this

This is consistence with $omega_0 = const.$ , i.e. no volume conservation. With LS coupling switched on.

Although it is not so clear, we can see the state avoid “crossing”. The s-state are unaffected by the LS coupling.

One deflect in the plot is that, the slope seem to be incorrect. And also, I tried to add volume conservation, but the result does not consistence….. (sad).

Axial Harmonic Oscillator – Nilsson Orbit

January 23, 2018

GoLuckyRyan Basic Harmonic Oscillator, Laguerre, Nilsson orbital Leave a comment

The Hamiltonian is

$\displaystyle H = -\frac{\hbar^2}{2m}\nabla^2 +\frac{m}{2}(\omega_\rho^2(x^2+y^2)+\omega_z z^2 )$

Use cylindrical coordinate, the Schrodinger equation is

$\displaystyle \left(-\frac{\hbar^2}{2m}\left(\frac{d^2}{dz^2} + \frac{1}{\rho}\frac{d}{d\rho}\left(\rho\frac{d}{d\rho}\right) + \frac{1}{\rho^2}\frac{d^2}{d\phi^2}\right) +\frac{m}{2}(\omega_\rho^2(\rho^2)+\omega_z z^2 ) \right) \Psi = E \Psi$

Where the energy is

$\displaystyle E = \hbar\omega_\rho(1+ n_x + n_y) + \hbar\omega_z\left(\frac{1}{2}+ n_z\right)$

Note that $n_x+n_y \neq n_\rho$ . One of the reason is there are 2 degree of freedom, it cannot be solely expressed into 1 parameter.

Set $\Psi = Z P \Phi$ ,

$\displaystyle -\frac{\hbar^2}{2m}\left(\frac{d^2Z}{dz^2} P \Phi + \frac{1}{\rho}\frac{d}{d\rho}\left(\rho\frac{dP}{d\rho} \right) Z \Phi + \frac{1}{\rho^2}\frac{d^2\Phi}{d\phi^2} Z P\right) \\ + \frac{m}{2}\omega_\rho^2\rho^2 ZP\phi+\frac{m}{2}\omega_z z^2 ZP\Phi= E ZP\Phi$

$\displaystyle -\frac{\hbar^2}{2m}\left(\frac{d^2Z}{dz^2} \frac{1}{Z} + \frac{1}{\rho}\frac{d}{d\rho}\left(\rho\frac{dP}{d\rho} \right) \frac{1}{P} + \frac{1}{\rho^2}\frac{d^2\Phi}{d\phi^2} \frac{1}{\Phi}\right) \\ + \frac{m}{2}\omega_\rho^2\rho^2 +\frac{m}{2}\omega_z z^2 = E$

The angular part, we can set

$\displaystyle \frac{d^2\Phi}{d\phi^2} \frac{1}{\Phi} = -m_\phi^2$

The solution is $\Phi(\phi) = \exp(i m_\phi \phi)$

The z-part is the 1D Harmonic oscillator, the solution is the Hermite function.

Thus the rest is

$\displaystyle -\frac{\hbar^2}{2m}\left( \frac{1}{\rho}\frac{d}{d\rho}\left(\rho\frac{dP}{d\rho} \right) \frac{1}{P} - \frac{m_\phi^2}{\rho^2}\right) + \frac{m}{2}\omega_\rho^2\rho^2 = \hbar\omega_\rho(1+ n_x + n_y)$

rearrange

$\displaystyle \frac{1}{\rho}\frac{d}{d\rho}\left(\rho\frac{dP}{d\rho} \right) - \frac{m_\phi^2}{\rho^2}P + \frac{m^2\omega_\rho^2}{\hbar^2}\rho^2 P + \frac{m\omega_\rho}{\hbar}(1+ n_x + n_y) P = 0$

Set a dimensionless constant $\alpha^2 = \frac{\hbar}{m \omega_\rho}$ , and $x = \rho/\alpha$

$\displaystyle \frac{1}{x}\frac{d}{dx}\left(x\frac{dP}{dx}\right) + \left( 2(1+n_x+n_y) - \frac{m_\phi^2}{x^2} - x^2 \right) P = 0$

Using the normalisation formula, $\int P^2 x dx = 1$ , set $u = P \sqrt{x}$

$\displaystyle x^2\frac{d^2u}{dx^2}+ \left( 2(1+n_x+n_y) - \frac{\frac{1}{4}-m_\phi^2}{x^2} - x^2 \right) u = 0$

Because of the long-range and short-range behavior, set

$\displaystyle u(x) = f(x) x^{m_\phi+\frac{1}{2}} \exp\left(-\frac{x^2}{2}\right)$

$\displaystyle x\frac{d^2f}{dx^2} + (1+2m_\phi-2x^2)\frac{df}{dx} + 2(n_x+n_y-m_\phi)x f= 0$

Set $y = x^2$

$\displaystyle y \frac{d^2f}{dy^2} + (m_\phi+1-y)\frac{df}{dy} + \frac{n_x+n_y-m_\phi}{2} f= 0$

Define $n_\rho = \frac{n_x+n_y-m_\phi}{2} \rightarrow n_x+n_y = 2 n_\rho + m_\phi = N - n_z$ .

$\displaystyle y \frac{d^2f}{dy^2} + (m_\phi+1-y)\frac{df}{dy} + n_\rho f= 0$

This is our friend again! The Laguerre polynomial! The complete solution is

$\displaystyle \Psi_{n_z n_\rho m_\phi}(z, \rho, \phi) \\ = \sqrt{\frac{1}{\alpha_z \alpha^2}}\sqrt{\frac{ n_\rho !}{2^n_z n_z! (m_\phi + n_\rho)!\sqrt{\pi^3}}} H_{n_z}\left(\frac{z}{\alpha_z}\right) \\ \exp\left(- \frac{1}{2}\left(\frac{z^2}{\alpha_z^2}+\frac{\rho^2}{\alpha^2}\right)\right) \left(\frac{\rho}{\alpha}\right)^{m_\phi} L_{n_\rho}^{m_\phi} \left(\frac{\rho^2}{\alpha^2}\right) \exp(i m_\phi \phi)$

where $\alpha_z^2 = \hbar/m/\omega_z$

The energy is

$\displaystyle E = \hbar\omega_z \left(\frac{1}{2} + n_z \right) + \hbar\omega_\rho \left( 2n_\rho + m_\phi + 1 \right)$

The quantum number $m_\phi$ has same meaning as $m_l$ in spherical case.

The notation for the state is

$|Nn_z m_l \rangle$

with the spin, the only good quantum number is the z-component of the total angular momentum $J = L+S$ long the body axis, i.e. $K = m_\phi \pm 1/2$ , thus, the state is

$|Nn_z m_l K \rangle$

Note that the total angular momentum $J^2$ is not a good quantum number in deformation, as the rotational symmetry is lost. However, the quantum number $K$ is linked with the angular momentum of the Nilsson single particle orbit. It is because when a particle has $m_j = K$ , the angular momentum must at least $j \geq K$ .

The above is a general solution for the harmonic oscillator in cylindrical coordinate. When $\omega_z = \omega_\rho$ , it reduce to spherical case.

According to P. Ring & P. Schuck (2004) (The Nuclear Many-Body Problem, P.68), the Hamiltonian can be expressed as a quadruple deform field by setting

$\displaystyle \omega_\rho^2 = \omega^2 \left(1+\frac{2}{3} \delta\right) \\ \omega_z^2 = \omega^2 \left(1-\frac{4}{3} \delta\right)$

$\displaystyle H = -\frac{\hbar^2}{2m}\nabla^2 + \frac{1}{2} m \omega^2 \left(r^2- \frac{2}{3}\sqrt{\frac{16\pi}{5}} \delta r^2 Y_{20}(\theta,\phi) \right)$

This Hamiltonian has similarity with the deformation

$\displaystyle R(\theta, \phi) = R_0 \left(1 + \alpha_00 + \sum_{\lambda=1}^{\infty}\sum_{\mu=-\lambda}^{\lambda} \alpha_{\lambda \mu} Y_{\lambda \mu}(\theta, \phi) \right)$

take the first quadruple term, and calculate $R^2(\theta, \phi)$

$R^2(\theta, \phi) = R_0^2 (1 + 2 \beta Y_{20}(\theta,\phi) )$

Compare, we have

$\displaystyle \beta = \frac{1}{3}\sqrt{\frac{16\pi}{5}}\delta \approx 1.05689 \delta$

The ratio

$\displaystyle \frac{\omega_z}{\omega_\rho} = \left(\frac{\alpha_\rho}{\alpha_z} \right) = \sqrt{\frac{1-\frac{4}{3}\delta}{1+\frac{2}{3}\delta}} = 1- \delta + \frac{1}{6}\delta^2 - \frac{5}{18}\delta^3...$

We need volume conservation

$\omega_x \omega_y \omega_z = \omega_0^3$

Thus,

$\displaystyle \omega = \omega_0 \left(\frac{1}{(1+2/3\delta)(1-4/3\delta)} \right)^{\frac{1}{6}} \approx \omega_0 \left(1+\frac{2}{9}\delta^2 + \frac{8}{81}\delta^3 \right)$

The energy is

$\displaystyle E = \hbar\omega \left( \sqrt{1-\frac{4}{3}\delta}\left(\frac{1}{2} + n_z \right) + \sqrt{1+\frac{2}{3}\delta} \left( N - n_z + 1 \right) \right) \\ \approx \hbar \omega_0 \left( \frac{3}{2} + \left(1 + \frac{1}{3}\delta\right) N - \delta n_z \right)$

Here is some plots with various $\delta = 0, 0.3, 0.5$

Next time, I will add the LS coupling and L-square term to recreate the Nilsson diagram. Also I will expand the solution from cylindrical coordinate into spherical coordinate. This unitary transform is the key to understand the single particle-ness.

3D Harmonic oscillator

January 19, 2018

GoLuckyRyan Basic Coulomb, Hermite, Laguerre Leave a comment

Set $x = r/\alpha$ The Schrodinger equation is

$\displaystyle \left(-\frac{\hbar^2}{2m} \nabla^2 + \frac{1}{2} m \omega^2 r^2 \right)\Psi = E \Psi$

in Cartesian coordinate, it is,

$\displaystyle -\frac{\hbar^2}{2m}\left( \frac{d^2}{dx^2}+\frac{d^2}{dy^2}+\frac{d^2}{dz^2} \right) \Psi + \frac{1}{2} m \omega^2 (x^2+y^2+z^2) \Psi = E \Psi$

We can set the wave function to be $\Psi(r, \Omega) = X(x) Y(y) Z(z)$

$\displaystyle \left( -\frac{\hbar^2}{2m}\frac{d^2X}{dx^2} + \frac{1}{2} m \omega^2 x^2 X \right) YZ + \\ \left( -\frac{\hbar^2}{2m}\frac{d^2Y}{dy^2} + \frac{1}{2} m \omega^2 y^2 Y \right) XZ + \\ \left( -\frac{\hbar^2}{2m}\frac{d^2Z}{dz^2} + \frac{1}{2} m \omega^2 z^2 Z \right) XY = E XYZ$

we can see, there are three repeated terms, we can set

$\displaystyle -\frac{\hbar^2}{2m}\frac{d^2X}{dx^2} + \frac{1}{2} m \omega^2 x^2 X = E_x X$

We decoupled the X, Y, Z. Each equation is a quadratic equation with energy

$\displaystyle E_x, E_y, E_z = \left(\frac{1}{2} + n \right) \hbar \omega$

and

$\displaystyle E_x + E_y + E_z = \left(\frac{3}{2} + n_x + n_y + n_z \right) \hbar \omega = \left(\frac{3}{2} + n \right) \hbar \omega = E$

The number of states for each energy level is

$\displaystyle C^{n_x+n_y+n_z+2}_2 = C^{n+2}_2 = \frac{(n+2)!}{n!2!}$

The first few numbers of states are 1, 3, 6, 10, 15, 21, 28, … The accumulated numbers of states are 1, 4, 10, 20, 35, 56, 84, … Due to the spin-state, the accumulated numbers of particles are 2, 8, 20, 40, 70, 112, 168, … The few magic numbers are reproduced.

The wave function is the product of the Hermite functions $H_n(x)$ and exponential function

$\Phi(x,y,z) = N H_{n_x} (x) H_{n_y}(y) H_{n_z}(z) \exp(-r^2/2)$

If we simply replace $(x,y,z) \rightarrow r( \cos(\phi) \sin(\theta), \sin(\phi) \sin(\theta) , cos(\theta) )$ , we can see the ground state consists of s-orbit, the 1st excited state consists of p-orbit, and the 2nd excited state consists of d-orbit.

To see more clearly, we can project the function onto spherical harmonic, which is fixed angular momentum, i.e.

$\langle \Phi(x,y,z) | Y_{lm}(\theta,\phi) \rangle = C_{lm} R_{nl}(r)$ .

where $C_{lm}$ is coefficients and $R_{nl}(r)$ is radial function.

To have a better understanding, the radial function has to be solved. The procedure is very similar to solving Coulomb potential.

$\displaystyle \left(-\frac{\hbar^2}{2m}\left(\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{d}{dr}\right)\right) + \frac{\hbar^2}{2m} \frac{L^2}{r^2} + \frac{1}{2} m \omega^2 r^2 \right)\Psi = E \Psi$

separate the radial part and angular part.

$\displaystyle \left( \frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{d}{dr}\right) - \frac{l(l+1)}{r^2} - \frac{m^2 \omega^2}{\hbar^2} r^2\right) R = - (2n+3) \frac{m \omega}{\hbar} R$

Set $\alpha^2 = \frac{\hbar}{m \omega}$

$\displaystyle \left( \frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{d}{dr}\right) + \frac{2n+3}{\alpha^2} - \frac{l(l+1)}{r^2} - \frac{r^2}{\alpha^4}\right) R = 0$

Set $x = r/\alpha$

$\displaystyle \left( \frac{1}{x^2}\frac{d}{dx}\left(x^2\frac{d}{dx}\right) + (2n+3) - \frac{l(l+1)}{x^2} - x^2\right) R = 0$

Set $u(x) = x R(x)$

$\displaystyle \frac{d^2 u }{d x^2} + \left( (2n+3) - x^2 - \frac{l(l+1)}{x^2} \right) u = 0$

as usual, the short range behaviour is $r^{l+1}$ , long range behaviour is $\exp(-x^2/2)$ , as stated in the Cartesian coordinate. Thus, we set

$\displaystyle u(x) = f(x) \exp\left(-\frac{x^2}{2}\right) x^{l+1}$

$\displaystyle x\frac{d^2f}{dx^2} + (2(l+1) -2x^2) \frac{df}{dx} + 2x(n-l) f(x) = 0$

with change of variable $y = x^2$ , the equation becomes

$\displaystyle y\frac{d^2f}{dy^2} + \left( l + \frac{1}{2} + 1 - y \right)\frac{df}{dy} + \frac{n-l}{2} f(y) = 0$

This is our friend, the Laguerre polynomial! In the Laguerre polynomial, $(n-l)/2 \geq 0$ must be non-negative integer. Now we set $k = (n-l)/2$ , than the energy is

$\displaystyle E = \hbar \omega \left( 2k + l + \frac{3}{2} \right)$

In order to have $n, l, k$ are integer, when

$n = 0 \rightarrow l = 0 \\ n = 1 \rightarrow l = 1 \\ n= 2 \rightarrow l = 0, 2 \\ n = 3 \rightarrow l = 1, 3 \\ n = 4 \rightarrow l = 0,2,4$

The overall solution without a normalization factor is

$\displaystyle \Psi_{nlm}(r, \theta, \phi) = r^l \exp\left(-\frac{r^2}{2\alpha^2}\right) L_{k}^{l+\frac{1}{2}}\left( \frac{r^2}{\alpha^2} \right) Y_{lm}(\theta, \phi)$

The normalization constant can be calculate easily using the integral formula of Laguerre polynomial.

$\displaystyle \int R^2(r) r^2 dr = 1$

$\displaystyle N^2 \int r^{2l} \exp\left(-\frac{r^2}{2\alpha^2}\right) \left( L_{k}^{l+\frac{1}{2}}\left( \frac{r^2}{\alpha^2} \right)\right)^2 r^2 dr = 1$

change of variable $y = r^2/\alpha^2 \rightarrow r=\alpha \sqrt{y}$

$\displaystyle dr = \frac{\alpha}{2\sqrt{y}}dy$

$\displaystyle N^2 \frac{\alpha^{2l+3}}{2} \int y^{l+\frac{1}{2}} \exp(-y) \left( L_k^{l+\frac{1}{2}} \right)^2 dy = 1$

The integration is

$\displaystyle N^2 \frac{\alpha^{2l+3}}{2} \frac{(k+l+\frac{1}{2})!}{k!} = 1$

$\displaystyle N^2 = \frac{2}{\alpha^{2l+3}} \frac{k!}{(k+l+\frac{1}{2})!}$

we can use

$\displaystyle \left(n+\frac{1}{2}\right) = \frac{(2n+1)!}{2^{2n+1}n!}\sqrt{\pi}$

Thus,

$\displaystyle N^2 = \frac{1}{\sqrt{\pi}\alpha^{2l+3}} \frac{k! (k+l)! 2^{2k+2l+3}}{(2k+2l+1)!}$

replace $k = (n-l)/2$ . The total wave function is

$\displaystyle \Psi_{nlm}(r, \theta, \phi) \\ =\sqrt{ \frac{1}{\sqrt{\pi}\alpha^{2l+3}} \frac{(\frac{n-l}{2})! (\frac{n+l}{2})! 2^{n+l+2}}{(n+l+1)!}} r^l \exp\left(-\frac{r^2}{2\alpha^2}\right) L_{k}^{l+\frac{1}{2}}\left( \frac{r^2}{\alpha^2} \right) Y_{lm}(\theta, \phi)$

Here are some drawing of the square of the wave functions. From the below is $n = 0, 1, 2, 3$ , from left to right, are s-orbit, p-orbit, d-orbit, f-orbit.

3-D Harmonic Oscillator.png

With the LS coupling, the spatial function does not affected, unless the coupling has spatial dependence. With the LS coupling, the good quantum numbers are $n, l, j, m_j$

Wave function in momentum space

October 31, 2017

GoLuckyRyan Basic Fourier transform, hydrogen, Laguerre, spherical Bessel Leave a comment

The wave function often calculated in spatial coordinate. However, in experimental point of view, the momentum distribution can be extracted directly from the experimental data.

The conversion between momentum space and position space is the Fourier transform

$\displaystyle \phi(\vec{k}) = \frac{1}{\sqrt{2\pi}^3} \int Exp\left(-i \vec{k}\cdot \vec{r} \right) \phi(\vec{r}) d\vec{r}$

Using the plane wave expansion

$\displaystyle \exp(i k\cdot r) = \sum_{l=0}^\infty (2l+1) i^l j_l(kr) P_l(\hat{k}\cdot\hat{r})$

$\displaystyle \exp(i k\cdot r) = 4\pi \sum_{l=0}^\infty \sum_{m=-l}^{l} i^l j_l(kr) Y_{lm}(\Omega_k) Y_{lm}^{*}(\Omega_r)$

Thus,

$\displaystyle \phi(\vec{k}) = \frac{4\pi}{\sqrt{2\pi}^3} \sum_{l=0}^\infty (-i)^l \sum_{m=-l}^{l} \int j_l(k r) Y_{lm}(\Omega_k) Y_{lm}^*(\Omega) \phi(\vec{r}) r^2 dr d\Omega$

where $j_l (x)$ is spherical Bessel function. Usually, due to conservation of angular momentum, the angular part can be separated from the spatial part. Let assume the wave function in position space is

$\phi(\vec{r}) = \psi(r) Y_{l_r m_r}(\Omega)$

$\phi(\vec{k}) = \psi(k) Y_{l_k m_k}(\Omega_k)$

Then we have

$\displaystyle \psi(k) = \frac{4\pi}{\sqrt{(2\pi)^3}} \int j_l(k r) \psi(r) r^2 dr$

$\displaystyle \phi(\vec{k}) = \psi(k) (-i)^l Y_{lm}(\Omega_k)$

where $l = l_r = l_k, m=m_r = m_k$ , due to the orthogonality of spherical harmonics.

For s, p, d, f-state, the spherical Bessel function is

$\displaystyle j_0(x) = \frac{\sin(x)}{x}$

$\displaystyle j_1(x) = \frac{\sin(x)}{x^2} - \frac{\cos(x)}{x}$

$\displaystyle j_2(x) = \sin(x)\frac{2-x^2}{x^3} - \cos(x)\frac{3}{x^2}$

$\displaystyle j_3(x) = \sin(x)\frac{15-6x^2}{x^4} - \cos(x)\frac{15-x^2}{x^3}$

For Hydrogen-like wave function, the non-normalized momentum distribution is

$\displaystyle \psi_{10}(k) = \frac{4 Z^{5/2}}{(k^2 + Z^2)^2}$

$\displaystyle \psi_{20}(k) = \frac{16 \sqrt{2}Z^{5/2}(4k^2-Z^2)}{(4k^2 + Z^2)^3}$

$\displaystyle \psi_{21}(k) =\sqrt{\frac{2}{3}}\frac{64 k Z^{7/2}}{(4k^2 + Z^2)^3}$

$\displaystyle \psi_{30}(k) = \frac{36 \sqrt{3} Z^{5/2} (81k^3-30k^2Z^2+Z^4)}{(9k^2 + Z^2)^4}$

$\displaystyle \psi_{31}(k) = \frac{144 \sqrt{6} Z^{7/2} (9k^3-kZ^2)}{(9k^2 + Z^2)^4}$

Here is the plot for momentum distribution $(\psi_{nl}(k))^2 k^2$ .

It is interesting that, the number of node decrease with higher angular momentum. But be-aware that it is only in atomic case, not a universal.

The higher the principle quantum number, the smaller of the spread of momentum. This is because, the spread of position wave function getting larger, and the uncertainty in momentum space will be smaller. This is a universal principle.

We also plot the Hydrogen radial function in here $\psi(r)^2 x^2$ , for reference,

Complete derivation from Schrodinger equation to Laguerre equation for Coulomb potential

October 25, 2017

GoLuckyRyan Basic Coulomb, hydrogen, Laguerre Leave a comment

This is for atomic case, that the eigen energy is $E_n = - \frac{Z^2}{2n^2}$ , and the Coulomb force is attractive.

The Hamiltonian is

$\displaystyle H = -\frac{\hbar^2}{2m}\nabla^2 - \frac{Ze^2}{4\pi\epsilon_0r}$

Separate the radial and angular part. The radial equation is

$\displaystyle \left( -\frac{\hbar^2}{2m}\frac{1}{r^2}\left(\frac{d}{dr} r^2 \frac{d}{dr} \right) - \frac{Ze^2}{4\pi\epsilon_0 r} + \frac{\hbar^2}{2m} \frac{l(l+1)}{r^2} \right) R(r) = E R(r)$

rearrange, using Atomic unit

$\displaystyle \left( -\frac{1}{2}\frac{1}{r^2}\left(\frac{d}{dr} r^2 \frac{d}{dr} \right) - E -\frac{Z}{ r} + \frac{1}{2} \frac{l(l+1)}{r^2} \right) R(r) = 0$

Since the normalization condition of $R(r)$ is $\int_0^\infty R(r)R(r) r^2 dr = 1$ , it is natural to define $u(r) = r R(r)$

$\displaystyle R = \frac{u}{r}, \frac{d}{dr}\left(r^2\frac{dR}{dr}\right) = r \frac{d^2u}{dr^2}$

The radial equation becomes

$\displaystyle -\frac{1}{2}\frac{d^2u}{dr^2} - \left( E + \frac{Z}{ r} - \frac{l(l+1)}{2r^2} \right) u(r)= 0$

Substitute the eigen energy $E_n = - \frac{Z^2}{2n^2}$

$\displaystyle \frac{d^2u}{dr^2} + \left( - \frac{Z^2}{n^2} + \frac{2Z}{ r} - \frac{l(l+1)}{r^2} \right) u(r)= 0$

Set

$\displaystyle x = \frac{2Z}{n} r, \frac{d}{dr} = \frac{d}{dx} \frac{2Z}{n}$

The radial equation becomes

$\displaystyle \left(\frac{2Z}{n}\right)^2\frac{d^2u(x)}{dx^2} + \left( - \frac{Z^2}{n^2} + \frac{4Z^2}{n x} - \frac{4Z^2}{n^2}\frac{l(l+1)}{x^2} \right) u(x)= 0$

Take out the $4Z^2/n^2$

$\displaystyle \frac{d^2u(x)}{dx^2} + \left(\frac{n}{ x} - \frac{1}{4} - \frac{l(l+1)}{x^2} \right) u(x)= 0$

Now, we know the short and long-range behaviour of $u(x)$ , assume it to be

$u(x) = \exp\left(-\frac{x}{2} \right) x^{l+1} y(x)$

Then the equation of $y(x)$ is

$\displaystyle x \frac{d^2y}{dx^2} + \left( 2(l+1) - x\right) \frac{dy}{dx} +\left( n - l- 1 \right) y(x) = 0$

This is the Laguerre equation

$\displaystyle x \frac{d^2y}{dx^2} + \left( \alpha+1 - x\right) \frac{dy}{dx} + m y(x) = 0$

where $\alpha = 2l+1$ and $m = n-l -1$ . Therefore, the solution of the radial equation is,

$R(r) = \frac{1}{r} A \exp\left(-\frac{x}{2} \right) x^{l+1} L_{n-l-1}^{2l+1}(x)$

where $A$ is normalization factor.

Notice that the Laguerre polynomial is only defined for $m \geq 0$ , thus, $n > l$ .

Using generating function to calculate integrals of Laguerre polynomial

October 25, 2017

GoLuckyRyan Math generating function, Laguerre Leave a comment

The generating function of Laguerre polynomial is

$\displaystyle \sum_n^\infty L_n^\alpha(x) t^n = \frac{1}{(1-t)^{\alpha+1}} \exp \left(-\frac{tx}{1-t}\right) =G_\alpha(t, x)$

We are going to evaluate the integral

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

Consider this integral

$\displaystyle \int_0^\infty x^{\alpha+k} e^{-x} G_\alpha(t,x) G_{\beta}(s,x) dx$

It is equal to

$\displaystyle \int_0^\infty x^{\alpha+k} e^{-x} \sum_{n,m}^{\infty} L_n^\alpha(x) L_m^\beta(x) t^n s^m dx \\ = \sum_{n,m}^{\infty} t^n s^m \int_0^\infty x^{\alpha+k} e^{-x} L_n^\alpha(x) L_m^\beta(x) dx \\ = \frac{1}{(1-t)^{\alpha+1}} \frac{1}{(1-s)^{\beta+1}} \int_0^\infty x^{\alpha+k} e^{-x} e^{-\frac{t x}{1-t} - \frac{s x}{1-s}} dx$

It is easy to calculate the last integral,

$\displaystyle \int_0^\infty x^{\alpha+k} e^{-x} e^{-\frac{t x}{1-t} - \frac{s x}{1-s}} dx = \int_0^\infty x^{\alpha+k} \exp \left(-\frac{1- t s}{(1-t)(1-s)}x \right) dx$

Replace variable,

$\displaystyle \int_0^\infty x^{\alpha+k} e^{-x} e^{-\frac{t x}{1-t} - \frac{s x}{1-s}} dx \\ = \frac{(1-t)^{(\alpha+1+k)}(1-s)^{(\alpha+1+k)}}{(1-t s)^{(\alpha+1+k)}}\int_0^\infty y^{\alpha+k} e^{-y} dy$

The integral becomes

$\displaystyle \sum_{n,m}^{\infty} t^n s^m \int_0^\infty x^{\alpha+k} e^{-x} L_n^\alpha(x) L_m^\beta(x) dx = \frac{(1-t)^{k}(1-s)^{(\alpha-\beta+k)}}{(1-t s)^{(\alpha+1+k)}} (\alpha+k)!$

Expand the right hand side,

$\displaystyle (1-t)^k = \sum_i^k (-1)^i \begin{pmatrix} k \\ i \end{pmatrix} t^i$

$\displaystyle (1-s)^{(\alpha-\beta+k)} = \sum_i^{\alpha-\beta+k} (-1)^i \begin{pmatrix} \alpha-\beta+k \\ i \end{pmatrix} s^i$

$\displaystyle (1-ts)^{-(\alpha+1+k)} = \sum_i^\infty (-1)^i \begin{pmatrix} -(\alpha+k+1) \\ i \end{pmatrix} (t s)^i = \sum_i^\infty \begin{pmatrix} \alpha+k+i\\ i \end{pmatrix} (t s)^i$

We used

$\displaystyle \begin{pmatrix} m \\ n \end{pmatrix} = \prod_{i=0}^{n-1} \frac{m-i}{n-i}$

$\displaystyle (-1)^n\begin{pmatrix} -m \\ n \end{pmatrix} = \prod_{i=0}^{n-1} \frac{m+i}{n-i} = \frac{(m+n-1)!}{n!(m-1)!}$

At the last step. The product

$\displaystyle \frac{(1-t)^{k}(1-s)^{(\alpha-\beta+k)}}{(1-t s)^{(\alpha+1+k)}} \\ = \sum_a^k \sum_{b}^{\alpha-\beta+k} \sum_c^\infty (-1)^{(a+b)} \begin{pmatrix} k \\ a \end{pmatrix} \begin{pmatrix} \alpha-\beta+k \\ b \end{pmatrix} \begin{pmatrix}\alpha+k+c \\ c \end{pmatrix} t^{(a+c)} s^{(b+c)}$

By comparing the order of $t^n s^m$ , the integral can be evaluated. With $a+c = n, b+c = m$

$\displaystyle \int_0^\infty x^{\alpha+k} e^{-x} L_n^\alpha(x) L_m^\beta(x) dx \\ = (\alpha+k)! \sum_a^k \sum_{b}^{\alpha-\beta+k} \sum_c^\infty (-1)^{(a+b)} \begin{pmatrix} k \\ a \end{pmatrix} \begin{pmatrix} \alpha-\beta+k \\ b \end{pmatrix} \begin{pmatrix} \alpha+k+c \\ c \end{pmatrix}$

When $\alpha = \beta$

$\displaystyle \int_0^\infty x^{\alpha+k} e^{-x} L_n^\alpha(x) L_m^\alpha(x) dx \\ = (\alpha+k)! \sum_a^k \sum_{b}^{k} \sum_c^\infty (-1)^{(a+b)} \begin{pmatrix} k \\ a \end{pmatrix} \begin{pmatrix} k \\ b \end{pmatrix} \begin{pmatrix} \alpha+k+c \\ c \end{pmatrix}$

For example, $k = 0$ ,

$\displaystyle \frac{(1-s)^{(\alpha-\beta)}}{(1-t s)^{(\alpha+1)}} = \sum_{b}^{\alpha-\beta} \sum_c^\infty (-1)^{b} \begin{pmatrix} \alpha-\beta \\ b \end{pmatrix} \begin{pmatrix} \alpha+c \\ c \end{pmatrix} t^{c} s^{(b+c)}$

When $\alpha \neq \beta, , n= m = 1$ , only $b = 0, c = 1$

$\displaystyle \int_0^\infty x^{\alpha} e^{-x} L_1^\alpha(x) L_1^\beta(x) dx = (\alpha)! \begin{pmatrix} \alpha-\beta \\ 0 \end{pmatrix} \begin{pmatrix} \alpha+1 \\ 1 \end{pmatrix} = (\alpha+1)!$

When $\alpha = \beta, , n = m$ , only $b = 0, c = n$

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

For $k = 1, \alpha = \beta, n= m$ , there are two terms, $a = b = 0, c = n$ and $a = b = 1, c = n -1$

$\displaystyle \int_0^\infty x^{\alpha+1} e^{-x} (L_n^\alpha(x))^2 \\ = (\alpha+1)! \left(\begin{pmatrix} \alpha+n+1 \\ n \end{pmatrix} + \begin{pmatrix} \alpha+n \\ n-1 \end{pmatrix} \right) = \frac{(\alpha+n)!}{n!}(2n+\alpha+1)$

We recovered all formula. Using generating function, we don’t have to worry about the condition of $n, m, \alpha, \beta, k$ . We only need to fulfill the sum with condition $a +c = n, b+c = m$ and $a, b, c \geq 0$ .

Hartree method for Helium ground state

October 17, 2017

GoLuckyRyan atomic, Basic, Math Clebsch-Gordon, Coulomb, Eigen System, Hartree-Fock, helium, hydrogen, Laguerre 7 Comments

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 a 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 matrixes, it is Hermitian, or symmetric.

We also notice that $G_{ij}^{hk}$ is a matrix of matrixes. The inner matrix runs for $hk$ indexes, while the bigger matrix runs for $ij$ index. This will simplify our minds.

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

In this simple 2-states example, index $h, k$ runs from 1 to 2. And remember the $G_{ij}^{hk}$ is a matrix of matrixes, using the identify $G_{ij}^{hk} = G_{hk}^{ij}$ the sum becomes:

$\sum_{h,k} a_h a_k G_{ij}^{hk} = a_1^2 G_{11} + a1 a_2 (G_{12} + G_{21}) + a_2^2 G_{22}$ ,

where the $G_{ij}$ represent the bigger matrix in $G_{ij}^{hk}$ .

The first trial wave function is 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)$ .

Integrations of Laguerre polynomial

October 15, 2017

GoLuckyRyan Math hydrogen, Laguerre 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$ , 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}$ .

Nuclear Physics 101

Notes on the original Nilsson paper

Diagonal elements for Nilsson orbital in spherical-spin function

Axial Harmonic Oscillator – Nilsson Orbit (III)

Axial Harmonic Oscillator – Nilsson Orbit

3D Harmonic oscillator

Wave function in momentum space

Complete derivation from Schrodinger equation to Laguerre equation for Coulomb potential

Using generating function to calculate integrals of Laguerre polynomial

Hartree method for Helium ground state

Integrations of Laguerre polynomial

Recent Posts

Calendar

Archives

Keywords

# of Visitors

Login