## Radial Density of 3D Harmonic Oscillator

From the previous post, we have the radial formula for the 3D harmonic oscillator.

$\displaystyle R_{nl}(r) \\ =\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)$

For $\alpha = 1 fm$, the $R^2$ are plotted below.

Assume all orbits are fully filled.

$\displaystyle \rho_n = \sum_{l}R^2_{nl} (2l+1)$

The red is $n=0$, the pink is $n = 1$, the blue is $n = 2$, the cyan is $n = 3$, and the green is $n = 4$.

Next, we are summing all the shells to get the nucleon density.

In here, the red is the 4He nucleon density function, the pink is the nucleus with all 1p-shell are filled, that is 16O. The blue is the sd-shell are filled, that is 40Ca. The cyan is fp-shell filled.

We can see some systematic trend there. The $n = odd$ nuclei has shape, which is somewhat similar to the Wood-Saxon shape. But for the $n=even$ nuclei, the shape is like a V.

Note that the $\alpha$ is fixed in this calculation. In fact, from Samuel Wong (2004) (Introductory Nuclear Physics),

$\displaystyle \hbar \omega \approx 41 A^{-1/3} \textrm{MeV}$

$\displaystyle \alpha = \sqrt{\frac{\hbar}{m \omega}} = \sqrt{\frac{\hbar^2 c^2}{mc^2 \hbar \omega}} \approx A^{1/6} \textrm{fm}$

## 3D Harmonic oscillator

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.

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$

## Kinematics of Transfer Reaction

Another surprise that no post on the kinematics of transfer reaction!

Suppose the reaction is $A(B, 1)2$, where $A$ is the incident particle, $B$ is the target nucleus that is at rest, $1, 2$ are the outgoing particle.

The Q-value for the reaction is the difference between initial mass and final mass

$Q = m_A+ m_B - m_1-m_2$

When the incident energy is too low, lower then the Q-value (we will see the relationship later) , the reaction is impossible, because the nuclei-1 and 2 cannot be created.

The equation for the reaction is

$P_A + P_B = P_1 + P_2$

where $P_i$ are 4-momenta. The number of freedom is 2. The total unknown is 8 from $P_1, P_2$, the masses of particle-1, and -2 give 2 constrains, the equations give 4 constrains.

$P_A = ( E_A , 0, 0, k_A) = ( m_A + T_A, 0, 0, \sqrt{2 m_A T_A + T_A^2})\\ P_B = (m_B, 0, 0, 0)$

where $T_A$ is the incident energy. We first transform the system into center of momentum (CM) frame. Since the total momentum is zero in the CM frame, we construct a CM 4-momentum

$P_{cm} = (P_A+P_B)/2 = ( m_A+ m_B + T_A, 0, 0, k_A)$

The Lorentz beta is $\beta = k_A/(m_A+m_B+T_A)$.

The invariance mass of the CM 4-momentum is unchanged, which is equal to the total energy of the particles in the CM frame. We can check

$P_A' = ( \gamma E_A - \gamma \beta k_A, 0, 0, -\gamma \beta E_A + \gamma k_A) \\ P_B' = ( \gamma m_B , 0, 0, -\gamma \beta m_B )$

The total energy in CM frame is

$E_{tot} = \gamma( m_A+m_B + T_A) - \gamma \beta k_A = \sqrt{(m_A+m_B+T_A)^2 - k_A^2}$

After the reaction, the momenta of the particle-1 and -2 are the same an opposite direction in the CM frame.Balancing the energy and momentum, the momentum is

$\displaystyle k = \frac{1}{2 E_{tot}} \sqrt{(E_{tot}^2 - (m_1+m_2)^2)(E_{tot}^2 - (m_1-m_2)^2)}$

In the equation, if either of the bracket is negative, the momentum cannot be formed. The smaller term is

$E_{tot}^2 - (m_1+m_2)^2 = (m_A+m_B+T_A)^2 - k_A^2 - (m_1+m_2)^2 \\ =(m_A+m_B+T_A)^2 - k_A^2 - (m_A+m_B-Q)^2 \\ = 2m_BT_A + 2(m_A + m_B)Q - Q^2$

If we set the line for this quantity be zero,

$\displaystyle 2m_BT_A + 2(m_A + m_B)Q - Q^2 = 0 \\ T_A = \frac{Q^2 - 2(m_A+m_B)}{2m_B} \\ = \frac{1}{2m_B} ( (Q-m_A-m_B)^2 - (m_A+m_B)^2) \\ = \frac{1}{2m_B} ( (m_1+m_2)^2 - (m_A+m_B)^2)$

Thus, if the $T_A$ is smaller then that value, the reaction is impossible.

The 4-momenta of the outgoing particles can be constructed.

$P_1' = ( \sqrt{m_1^2 - k^2} , k \cos(\phi) \sin(\theta) , k \sin(\phi) \sin(\theta) , k \cos(\theta)) \\ P_2' = ( \sqrt{m_2^2 - k^2} , -k \cos(\phi) \sin(\theta) ,- k \sin(\phi) \sin(\theta) , -k \cos(\theta))$

Back to Lab frame,

$P_1 = \begin{pmatrix} \gamma \sqrt{m_1^2 + k^2} + \gamma \beta k \cos(\theta) \\ k \cos(\phi) \sin(\theta) \\ k \sin(\phi) \sin(\theta) \\ \gamma \beta \sqrt{m_1^2 - k^2} + \gamma k \cos(\theta)\end{pmatrix}$

For simplicity, set $\phi = 0$,

$P_1 = \begin{pmatrix} \gamma \sqrt{m_1^2 + k^2} + \gamma \beta k \cos(\theta) \\ k \cos(\phi) \sin(\theta) \\ 0 \\ \gamma \beta \sqrt{m_1^2 - k^2} + \gamma k \cos(\theta) \end{pmatrix}$

We can see, the locus in the momentum space is

$( k cos(\phi) \sin(\theta) , 0, \gamma \beta \sqrt{m_1^2 - k^2} + \gamma k \cos(\theta))$

This is an ellipse with length $k (1, \gamma)$, and centered at $(0, \gamma \beta \sqrt{m_1^2 + k^2})$. If the relativistic effect is small, it is a circle.

When the ellipse long axis $\gamma k$ is longer then the center position $\gamma \beta \sqrt{m_1^2 + k^2}$, the scattering angle can be to whole circle, i.e.

$\displaystyle \gamma \beta \sqrt{m_1^2 + k^2} < \gamma k \rightarrow \gamma \beta m_1 < k$

in here, we can treat $\gamma \beta m_1 = k_\beta(1)$, this is the momentum “gain” from shifting from CM frame to Lab frame. By compare $k_\beta(1)$ or $k_beta(2)$ to $k$ we can know the span of scattering angle.

Here is an example,

## Finite Spherical Square Well with spin-orbital potential

Last time, we studied the finite spherical square well and calculated the energy levels for difference angular momentum. In that case, the good quantum numbers (quantities that is commute with the Hamiltonian, i.e. conserved, invariance with time) are the principle quantum number $n$, which dictated the number of node and the angular momentum $l$, and the spin state $s$. A eigen state can be denoted as $|nl m_l s m_s\rangle$. The degeneracy for each state is $2(2l+1)$.

With the spin-orbital potential

$\displaystyle V_{so} = a L \cdot S$

The angular momentum $m_l$ and spin $m_s$ are not good quantum numbers. The sum, total angular momentum $j = l + s$ and it z-component $m_j$ are. The new eigen state is $|n j m_j l s\rangle$.

The Hamiltonian inside the well is

$\displaystyle H = H_0 + V_{so} = -\frac{\hbar^2}{2m} \frac{1}{r^2} \frac{d}{dr}r^2\frac{d}{dr} - |V_0| + \frac{\hbar^2}{2m} \frac{1}{r^2}L^2 + a L\cdot S$

The energy for $H_0$ dose not change, as the $L^2$ is still a good quantum number. The spin-orbital coupling is evaluated like

$\displaystyle L\cdot S = \frac{1}{2}(J^2 - L^2 - S^2)$

then,

$\displaystyle \langle n j m_j l s| L\cdot S |n j m_j l s\rangle = \frac{1}{2} ( j(j+1) - l(l+1) - s(s-1) )$

Since the spin has only 2 spin states, thus,

$\displaystyle \langle L\cdot S \rangle \\ = \frac{1}{2} l , j = l+\frac{1}{2} \\ = -\frac{1}{2}(l+1), j = l - \frac{1}{2}$

Note that, the mean energy of the L-states are unchanged. The degeneracy for $j_+ = l+1/2$ is $2(l+1)$, for $j_- = l-1/2$ is $2l$. Also, the s-states are not affected, because $l = 0$.

Thus, the wave vector inside the well $k$ becomes

$\displaystyle k = \sqrt{\frac{2m}{\hbar^2 }(|V_0|- \beta \langle L\cdot S\rangle) - |E|)}$

the wave vector outside the well $k'$

$\displaystyle k' = i \kappa = i \sqrt{\frac{2m}{\hbar^2 }(\beta \langle L\cdot S\rangle) + |E|)}$

And the rest is matching the boundary condition.

We can also add one more term, an additional $\gamma L^2$ term so that the centroid of the L-state is shifted. The energy for this potential is $\gamma l(l+1)$.

Here is the energy levels for $\beta = -1$ and $\gamma = -0.5$.

We can see, we reproduced the shell ordering 0s1/2, 0p3/2, 0p1/2, 0d5/2, 1s1/2, 0d3/2, 0f7/2, 1p3/2, 1p1/2, 0f5/2….

The magic number 2, 8, 20, and 40 are reproduced. May be, if I adjust the strength of the spin-orbit coupling, my get the magic number 28 and 50. Since the energy level of the s-states are fixed, I need to adjust the shift of the centroid of the others, and the splitting.

Here is the result of $\beta = -1.2, \gamma = -0.2$

It seems that the magic number of 28, which is from the isolation of 0f7/2, is recreated. And the isolation of 0g9/2, created the magic number of 50.

## 1-D Harmonic Oscillator

It is quite surprised that This is not here!

The Hamiltonian is

$\displaystyle H=\frac{P^2}{2m} + \frac{m \omega^2}{2} X^2$

Lets define creation operator $A^\dagger$ and destruction operator $A$, or the ladder operators, such that,

$\displaystyle A^\dagger =\sqrt{\frac{1}{2m}} P + i \sqrt{\frac{m \omega^2}{2}} X$

$\displaystyle A =\sqrt{\frac{1}{2m}} P - i \sqrt{\frac{m \omega^2}{2}} X$

Then we can see

$\displaystyle H = A^\dagger A + \frac{\hbar \omega}{2} = A A^\dagger - \frac{\hbar \omega}{2}$

and the commutative relations

$\displaystyle A^\dagger A - A A^\dagger = -\hbar \omega$

$\displaystyle HA^\dagger - A^\dagger H = \hbar \omega A^\dagger$

$\displaystyle HA - A H = -\hbar \omega A$

From the last two equations, we can see

$\displaystyle H |n\rangle = E_n |n\rangle \\ A H|n\rangle = (HA+\hbar \omega A) |n\rangle = E_n A|n\rangle \\ HA|n\rangle = (E_n-\hbar \omega) A|n\rangle$

This shows that, if $|n\rangle$ is a state with energy $E_n$, then $A|n\rangle$ is also a state with energy $E_n - \hbar \omega$. Similar for $A^\dagger |n\rangle$, the energy will be $E_n + \hbar \omega$.

Suppose the ground state is $|0\rangle$, which is is lowest energy state. If we apply the destruction operator, it will be gone, i.e.

$\displaystyle A|0 \rangle = 0 \\ A^\dagger A|0\rangle = 0 \\ \left(H- \frac{\hbar\omega}{2}\right)|0\rangle \\ H|0\rangle = \frac{\hbar\omega}{2}|0\rangle = E_0 |0\rangle$

Thus, the ground state energy is $\hbar \omega /2$!

We can donate

$\displaystyle A^\dagger |n\rangle = U_n |n+1\rangle \\ A |n\rangle = D_n |n-1\rangle$

Using $\displaystyle A^\dagger A - A A^\dagger = -\hbar \omega$, we have

$\displaystyle (A^\dagger A - A A^\dagger ) |n \rangle \\ = A^\dagger D_n |n-1 \rangle - A U_n|n+1 \rangle \\ = (U_{n-1}D_n- U_n D_{n+1} )|n\rangle = - \hbar \omega |n\rangle$

or

$\displaystyle U_n D_{n+1} = U_{n-1}D_n + \hbar \omega$

Consider above equation on $|0\rangle$, such that $A|0\rangle = 0$, thus,

$\displaystyle (A^\dagger A - A A^\dagger ) |0 \rangle = - U_0 D_{1} )|0\rangle = - \hbar \omega |0\rangle$

Thus, $U_0 D_1 = \hbar \omega$, then, in general, we have

$\displaystyle U_n D_{n+1} = (n+1) \hbar \omega$

Then

$\displaystyle A^\dagger A |n\rangle = U_{n-1} D_n |n\rangle \\ A A^\dagger |n\rangle = U_{n} D_{n+1} |n\rangle$

The simplest solution is

$U_n = \sqrt{\hbar \omega}\sqrt{n+1} \\ D_n =\sqrt{\hbar \omega}\sqrt{n}$

And, if we treat the state $|n\rangle$ contains $n$ oscillators, we have a number operator $N = A^\dagger A$, such that

$A^\dagger A |n \rangle = n \hbar \omega |n \rangle$

some people will normalized the ladder operators with $\sqrt{\hbar \omega}$. For me, I don’t want to make the starting point be so “artificial”. We can summarized the 1-D Harmonic oscillator as

Now, we are going to solve the wave function. One way is use the ladder operator,

$\displaystyle A|0\rangle = 0 \\ \left( \sqrt{\frac{1}{2m}} P - i \sqrt{\frac{m \omega^2}{2}} X \right) \phi_0 = 0 \\ \frac{d\phi_0}{dx} = - \frac{m \omega}{\hbar} x \phi_0 = 0$

The solution is

$\displaystyle \phi_0 = N \exp\left(-\frac{m \omega}{2\hbar}x^2\right)$

Another way is the Schrodinger equation is

$\displaystyle -\frac{\hbar^2}{2m}\frac{d^2\phi}{dx^2} + \frac{1}{2}m\omega^2 x^2 \phi = E \phi = \frac{2n+1}{2}\hbar \omega \phi$

$\displaystyle \frac{d^2\phi}{dx^2} - \frac{m^2 \omega^2}{\hbar^2} x^2 \phi = (2n+1)\frac{m \omega}{\hbar} \phi$

Define $\alpha = \sqrt{\frac{\hbar}{m \omega}}$, the dimension of $\alpha$ is length. Set $u = x/ \alpha$, then

$\displaystyle \frac{d^2\phi}{du^2} + ((2n+1) - u^2) \phi = 0$

Since we know that the function must decay outside the well, set $\phi = f(y) \exp(-y^2/2)$, then

$\displaystyle \frac{d^2f}{du^2} -2 u \frac{df}{du} +2 nf = 0$

The solution is Hermite polynomial $H_n(u)$ of order $n$.

Assume

$\displaystyle H_n(x) = \sum_{i = 0}^{\infty} a_i x^i$

sub into the Hermite equation,

$\displaystyle a_{i+2} = \frac{2(i-n)}{(i+1)(i+2)} a_i$

First,  the even and odd series are separated. And either the odd-series or the even-series are converge, as the ratio

$\displaystyle \lim_{i->\infty} \frac{a_{i+2}}{a_i} = 2$

However, the series can be terminated when $i = n$. Thus, the Hermite polynomial has either even terms or odd terms, but not mixed.

When $n = 0$, $a_2 = 0$ and all later even terms are zero. The equation is $f''(x) - 2 x f'(x) = 0$, which only support even function. Thus, $H_0(x) = a_0$ is a constant. Thus, we can see, when $n$ is even (odd), Hermite polynomial is even (odd) with only even (odd) $n$.

Hermite polynomial are orthogonal with $\exp(-x^2)$, thus, the wave function $\phi_n(x) = H_n(x) \exp(-x^2/2)$ are orthogonal.

The normalized wave function is

$\displaystyle \phi_n(x) = \frac{1}{\sqrt{2^n n! \sqrt{\pi}}} H_n(x) \exp\left(-\frac{x^2}{2}\right)$

or

$\displaystyle \phi_n\left(\frac{x}{\alpha}\right) = \frac{1}{\sqrt{\alpha 2^n n! \sqrt{\pi}}} H_n\left(\frac{x}{\alpha}\right) \exp\left(-\frac{x^2}{2\alpha^2}\right)$

## Finite Spherical Square Well

Hello, everyone, in order to calculate deuteron by Hartree-Fock method, I need a basis. The basis of infinite spherical square well is too “rigid”, that it has to “extension” to non-classical region. Beside of the basis of Wood-Saxon potential. The finite spherical square well is a good alternative. The radial equation is basically the same as infinite spherical square well.

The potential is

$\displaystyle V(r) = -|V_0|, r\leq a$

$\displaystyle V(r) = 0, r > a$

Within the well, the wave vector is

$\displaystyle k = \sqrt{\frac{2m}{\hbar^2} (|V_0 |-|E|) }$

, outside the well, the wave vector is

$\displaystyle k' = i \kappa = i \sqrt{\frac{2m}{\hbar^2} |E| }$

The solution is spherical Bessel function. Since the Bessel function of the first and second kind are oscillating like sin or cosine function. To form a decay function when $r > a$, we need the Hankel function with complex argument.

$\displaystyle h_n( i \kappa r) = h_n(x) = - (i x)^n \left( \frac{1}{x} \frac{d}{dx}\right)^n \left(\frac{\exp{(-x)}}{x} \right)$

To make it real, we need a factor $i^n$.

The boundary conditions are continuity and differential continuity.

$\displaystyle j_l(ka) = A i^l h_l(i\kappa a)$

$\displaystyle \left(\frac{d}{dr}j_l(kr)\right)_{r=a} = A i^l \left(\frac{d}{dr}h_l(i\kappa a) \right)_{r=a}$

These two conditions solved for two parameters $A$ and $E$. However, I cannot find an analytical solution to the energy $E$.

If the potential depth is 60 MeV for proton. Radius is 1 fm for a light nuclei. Set $\hbar = 1, m = 1$.

The result is follow,

The 1st column is 1s, 1p, 1d, 1f, and 1g. The 2nd column is 2s, 2p, 2d, 2f, and 2g. The 3rd column is 3s, 3p, and 3d.

By compare with infinite square well,

The energies are lower in finite well, because the wave functions can spread-out to non-classical region, so that the wave length is longer and energy is lower.

In this example, the 3d and 2g orbits are bounded (of course, all orbits in infinite well are bound). This is not because of the depth of the well, but the boundary of the well. In other word, to bring down an orbit, the wave function has to spread out, that is connected with the neutron-halo.

## Wave function in momentum space

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

or

$\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,