spherical Bessel | Nuclear Physics 101

Variational method for Woods-Saxon potential

August 31, 2021

GoLuckyRyan Basic Hankel, Runge-Kutta, Schrödinger equation, spherical Bessel, Variational method Leave a comment

In the past, we solved the radial wavefunction for the Woods-Saxon potential by Runge-Kutta method. But there is another method based on linear algebra, that all wave function can be expanded as a linear combination of a set of basis. In here, we use a free wave as basis, which is the spherical Bessel function.

to be annoying due to my poor memory , lets state the Schrödinger equation again,

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

Separating the radial part

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

replace $R(r) = u(r) /r$

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

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

simplify

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

We can define a “free wave” Hamiltonian as

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

In this post, we explored the orthogonality of the the Riccati-Bessel function, and it is the eigen state for $H_0$ with eigenvalue of $-k_n^2$ .

$\displaystyle b_k(r) = \sqrt{\frac{2}{pi}} \hat{j_l}\left(k r\right), ~~~ H_0 b_k(r) = -k^2 b_k(r)$

$\displaystyle \left< b_k | b_k' \right> = \delta(k-k')$

Thus, the Riccati-Bessel function can be served as the set of basis functions that span the whole space.

The total Hamiltonian is then

$\displaystyle \left( -\frac{\hbar^2}{2m} \frac{d^2}{dr^2} + \frac{\hbar^2}{2m} \frac{l(l+1)}{r^2} + V(r) \right) u(r) \\~~~~~~~~~~~~~~~~= \left( -\frac{\hbar^2}{2m} H_0 + V(r) \right) u(r) = E u(r)$

Since the solution can be expressed as a linear combination of the basis functions, then

$\displaystyle u(r) = \sum_n a_n b_n(r)$

$\displaystyle \left( -\frac{\hbar^2}{2m} H_0 + V(r) \right) \sum_n a_n b_n(r) = E\sum_n a_n b_n(r)$

Multiple $b_m(r)$ and integrate,

$\displaystyle \sum_{n} a_n \left< b_m \left| -\frac{\hbar^2}{2m}H_0 + V(r) \right|b_n \right> = E a_m$

The above is a matrix equation

$\displaystyle \left< -\frac{\hbar^2}{2m} H_0 + V(r) \right> \vec{a} = E \vec{a}, ~~ \vec{a} = (a_1, a_2,... a_N)$

Everything seems OK. But when I test the theory using Mathematica, it fail. For simplicity, we use the s-orbital, and the basis is

$\displaystyle b_n(r) = \sqrt{\frac{2}{\pi}} \sin\left( \frac{2\pi}{a} n r\right)$

The Woods-Saxon parameters are $V_0 = -50 \mbox{MeV}, V_{so} = 22 \mbox{MeV}, R_0 = R_{so} = 4 \mbox{fm}, a_0 = a_{so} = 0.67 \mbox{fm}$ .

Using Runge-Kutta method, the ground state energy is -36.09 MeV, and the wave function is

We set $a = 40~\mbox{fm}$ , and $n = 1, 2, ...., 40$ . The matrix elements is

$\displaystyle \left<n \left| -\frac{\hbar^2}{2m_n}H_0 + V(r) \right| m \right> = \delta_{nm} \frac{\hbar^2 k^2}{2m_n} + \left<V(r) \right>, ~~ k = \frac{2\pi}{a} n$

where $m_n$ is neutron mass.

After calculate the matrix, find the eigenvalue and eigenstate. The lowest eigen energy is -465.18 MeV, which indicate something wrong. And the wave function, constructed by the eigenstate look likes

Ground state constructed by basis functions

The shape of the wavefunction looks similar to the correct wavefunction, but narrower and larger amplitude.

The next eigen wave function also look like 1s-orbital. The spectrum of the eigenvalues, which are the eigen energies, has 7 bound states.

I don’t know the problems….

One possible reason could be that the basis function is always $b_n(a) = 0$ , but the real solution is not necessary to be $\psi(a) = 0$ . In fact, the 0s-wave function never be zero except $r= 0$ .

I use the same basis and extract the coefficient,

$\displaystyle \psi(r) \approx \sum_{n}^{40} C_n \frac{\sin(k_n r)}{r}, ~ C_n = \int_0^\infty \psi(r) \sin(k_n r) r dr, ~ k_n = \frac{2\pi}{a}n$

And it can reproduce the original wave function at $r < 10 \mbox{fm}$

Since the basis function is a sine wave and we only sum limited term, the long range oscillation can not be suppressed.

Another possible reason is that the basis is not normalizable. This root at the orthogonality of Bessel function that

$\displaystyle \int_0^\infty J_\nu^2(kx) x dx = \infty$

An-another possible reason is that, the basis is for unbound state. But in this post, we used a not-normalizable, unbounded basis to construct the bound state wave function. It can form the solution correctly at short-range, but it is periodic. This is the nature of series expansion, always periodic.

In this sense, we are also expanding the solution into a series, that the series only produce periodic function, therefore, a correct way to do is using the spherical Hankel transform, which is transform the Schrodinger equation into momentum space.

Orthogonality of Spherical Bessel function J

August 30, 2021

GoLuckyRyan Basic Bessel function, Fourier Series, Fourier transform, Hankel, spherical Bessel 4 Comments

The orthogonality of Bessel function for finite range is

$\displaystyle \int_{0}^{a} J_{\nu}\left( \alpha_{\nu n} \frac{x}{a} \right) J_{\nu}\left( \alpha_{\nu m} \frac{x}{a} \right) x dx = \frac{a^2}{2} J_{\nu+1}^2(\alpha_{\nu n}) \delta_{nm}$

where $\alpha_{\nu n}$ is the n-th root of the Bessel function $J_{\nu}$ . The above Bessel function is scaled to length $a$ with n-th node. For example, $a = 1, \nu = 0$

Simply the integration by setting $k_n = \alpha_{\nu n} / a$

$\displaystyle \int_{0}^{a} J_{\nu} (k_n x)J_{\nu}(k_m x) x dx = \frac{a^2}{2} J_{\nu+1}^2(\alpha_{\nu n}) \delta_{nm}$

The relation between Riccati-Bessel function, spherical Bessel function, and Bessel function are:

$\displaystyle \hat{j_l}(x) = x j_l(x), ~~ j_l(x) = \sqrt{\frac{\pi}{2x}} J_{l+1/2}(x)$

The orthogonal condition of the spherical Bessel function should be

$\displaystyle \int_{0}^{a} j_l( \kappa_n x) j_l(\kappa_m x) x^2 dx , ~~ \kappa_n = \frac{\alpha_{l+1/2, n}}{a}$

From the Bessel Function orthogonality,

$\displaystyle \int_{0}^{a} J_{l+1/2} (\kappa_n x)J_{l+1/2}(\kappa_m x) x dx = \frac{2}{\pi} \kappa_n \int_0^a x^2 j_l(\kappa_n x) j_l(\kappa_m x) dx$

$\displaystyle \int_0^a j_l (\kappa_n x) j_l( \kappa_m x) x^2 dx = \frac{a^3}{2} j_{l+1}^2(\alpha_{l+1/2, n} ) \delta_{nm}$

In term of Riccati-Bessel function, $\hat{j_l}(kx) = kx j_{l}(kx)$ ,

$\displaystyle \int_0^a \hat{j_l} (\kappa_n x) \hat{j_l}( \kappa_m x) dx = \kappa_n^2 \frac{a^3}{2} j_{l+1}^2(\alpha_{l+1/2, n} ) \delta_{nm} = \frac{a}{2} \hat{j_{l+1}}(\alpha_{l+1/2, n}) \delta_{nm}$

Using the orthogonality, the spherical Bessel function can be basis that span the whole space with $f(a) = 0$ . We can also normalized the spherical Bessel function. For Riccati-Bessel function, it span the whole space with $f(0)=0, f(a)=0$ . Finite-range function expands as Bessel function is called Fourier-Bessel series.

The infinite range orthogonality of Bessel function J is (according to wiki and many reference, why?)

$\displaystyle \int_0^\infty J_\alpha( u x) J_\alpha( vx) x dx = \frac{1}{u}\delta(u-v), ~~ \alpha > -\frac{1}{2}$

in term of spherical Bessel function

$\displaystyle \int_0^\infty j_{\alpha}(u x) j_\alpha(vx) x^2 dx = \frac{\pi}{2u^2} \delta(u-v), ~~ \alpha > -1$

and using $\hat{j_l}(kx) = kx j_{l}(kx)$ , in term of Riccati-Bessel function is

$\displaystyle \int_0^\infty \hat{j_{\alpha}}(u x) \hat{j_\alpha}(vx) dx = \frac{\pi}{2} \delta(u-v), ~~ \alpha > -1$

For example, the spherical Bessel J of order zero,

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

Let’s integrate it

$\displaystyle \int_0^\infty \frac{\sin(x)}{x} \frac{\sin(x)}{x} x^2 dx = \int_0^\infty \sin^2(x) dx$

which is infinite, and the right-hand side is a delta function $\delta(0) = \infty$ .

The problem with the above formulae is, that the orthogonal condition is OK. but it should be also normal, i.e. orthonormal, to be a normalizable basis.

The Hankel Transform is

$\displaystyle F_\nu(k) = \int_0^\infty f(r) J_\nu(kr) r dr , ~~ \nu > -\frac{1}{2}, k \geq 0$

and this can change to spherical-Bessel function as,

$\displaystyle \psi (k) = \sqrt{\frac{2}{\pi}}\int_0^\infty \psi(r) j_l(kr) r^2 dr$

The foundation for the transformation is the orthogonality of the Bessel function. Substitute $f(r)$ as Bessel function.

$\displaystyle \int_0^\infty J_0(h r) J_0(kr) r dr =\frac{1}{h} \delta (k-h)$

So, what do I miss?

For example, $J_\nu(x)$ is oscillating around zero, the square of it is always positive, and with a weighting $x$ , it approach to $\sin^2(x)$ with a phase shift,

$\displaystyle J_\nu^2(x) x \xrightarrow{x\rightarrow \infty} \frac{2}{\pi} \sin^2(x + \theta_\nu)$

And the integrate of $\sin^2(x)$ is $\pi/2$ . Thus, the integral of Bessel function with weight $x$ is an infinite sum for between each zeros, and each terms is ~1. How this infinite sum becomes 1?

Recall the plane wave expansion,

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

As we show it before, it is a special case of the Hankel transform, or belong the general Fourier Transform.

Hankel Transform

May 6, 2021

GoLuckyRyan Basic Bessel function, momentum, spherical Bessel Leave a comment

I recently found out that the transformation from radial wave function to momentum radial wave function is the Hankel Transform in special case.

The Hankel transform is defined as

$\displaystyle F_\nu(k) = H_\nu(f(r)) = \int_{0}^{\infty} f(r) J_\nu(kr) r dr$

where $J_\nu(x)$ is the Bessel function J.

Since the Bessel function and the spherical Bessel function is related as

$\displaystyle j_l(x) = \sqrt{\frac{\pi}{2 x}} J_{l+1/2}(x)$

Set $\nu = l + 1/2$

$\displaystyle H_{l+1/2}(f(r)) = \sqrt{\frac{2k}{\pi}} \int_{0}^{\infty} f(r) \sqrt{r} \sqrt{\frac{\pi}{2 k r}} J_{l+1/2}(kr) r dr \\ = \sqrt{\frac{2k}{\pi}} \int_{0}^{\infty} f(r) \sqrt{r} j_l(kr) r dr$

Thus the transformation from spatial function to momentum function is

$\displaystyle \psi(k) = \sqrt{\frac{2}{\pi}}\int_0^\infty \psi(r) j_l(kr ) r^2 dr = \frac{1}{\sqrt{k}} H_{l+1/2}( \psi(r) \sqrt{r} )$

The Hankel transform is a built-in function in Mathemtica

HankelTransform[ f(r) sqrt[r], r, k, l+1/2 ] / sqrt[k]

And also the inverse Hankel transform

InverseHankelTransfrom[ F[k] sqrt[k], k, r, l+1/2 ] / sqrt[r]

But the built-in Hankel transform does not support digitized data.

Extracting phase shift

March 5, 2021

GoLuckyRyan Basic phase shift, Riccati-Bessel, spherical Bessel Leave a comment

After solving the Schrodinger equation numerically, we have radial function in numerical form. The radial function is a combination of the spherical Bessel function $j_n , y_n$ :

$\displaystyle R_l(r) = a_l j_l(kr) + b_l y_l(kr)$

or in Ricatti-Bessel function

$\displaystyle u_l(r) = a_l \hat{j_l}(kr) + b_l \hat{y_l}(kr)$

A textbook method is matching the logarithmic derivative at the boundary of the external region.

$\displaystyle \left[ \frac{1}{u_l(r)}\frac{du_l(r)}{dr} \right]_{r\rightarrow R^+} = \left[ \frac{1}{u_l(r)} \frac{du_l(r)}{dr}\right]_{r\rightarrow R^+}$

Another method is using the asymptotic property of the Ricattie-Bessel function:

$\displaystyle \hat{j_l}(r \rightarrow \infty) = \sin \left( r - l\frac{\pi}{2} \right) \\ \hat{y_l}(\rightarrow \infty) = (-1)^{l+1} \cos \left( r + l\frac{\pi}{2} \right)$

The asymptotic form of the radial solution is

$\displaystyle u_l(r \rightarrow \infty) = a_l \sin \left( kr - l\frac{\pi}{2} \right) + b_l (-1)^{l+1} \cos \left(kr + l\frac{\pi}{2} \right)$

We can find the coefficient $a_l , b_l$ using the orthonormal property of the Ricattie-Bessel function.

Since the free wave solution is $u_l(r) = \hat{j_l}(kr)$ , so, the phase shift $\phi_l$ is

$\displaystyle u_l(r \rightarrow \infty) \propto \sin \left( kr - l\frac{\pi}{2} + \phi_l \right)$

Using the analog of

$\displaystyle a \sin(x) + b\cos(x) = \sqrt{a^2+b^2} \sin(x + \phi) , \tan(\phi) = \frac{b}{a}$

We can also calculate the phase shift using the coefficient of $a_l, b_l$ , assuming they are real.

We can define $h_l =\sqrt{a_l^2+b_l^2}$ , and set $a_l = h_l \cos(\phi_l) , b_l = h_l \sin(\phi_l)$ , so that $\tan(\phi_l) = \frac{b}{a}$ , thus the asymptotic form of the radial solution is

$\displaystyle u_l(r \rightarrow \infty) = h_l \sin \left( kr - l\frac{\pi}{2}- \phi_l \right)$

In this post, the unitarity of the S-matrix is defined as

$\displaystyle S_l = \frac{a+ib}{a-ib} = e^{2i\delta_l}$

We can show that $\delta_l = \phi_l$ ,

$\displaystyle S_l = \frac{a+ib}{a-ib} = \frac{a^2-b^2}{a^2+b^2} - i \frac{2ab}{a^2+b^2} = \cos(2\delta_l) + i \sin(2\delta_l) = e^{2i\delta_l}$

so, that

$\displaystyle \cos(2\delta_l) = \frac{a^2-b^2}{a^2+b^2}, \sin(2\delta_l) = \frac{2ab}{a^2+b^2}$

Thus

$\displaystyle \tan(2\delta_l) = \frac{2 \tan(\delta_l)}{1-\tan^2(\delta_l)} = \frac{2ab}{a^2-b^2} = \frac{2\frac{b}{a}}{1-\left(\frac{b}{a}\right)^2} = \frac{2 \tan(\phi_l)}{1-\tan^2(\phi_l)}$

S-matrix and scattering amplitude

April 16, 2020

GoLuckyRyan Basic, Scattering phase shift, Riccati-Bessel, S-matrix, scattering, spherical Bessel Leave a comment

In this post and this post talk about the scattering and phase shift. But the relation is not shown clearly.

The central idea is solving the Schrödinger equation

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

rewrite,

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

The homogeneous solution ( i.e. $V(r) = 0$ ) is a plane wave $e^{i\vec{k} \cdot \vec{r} }$ . Thus, the general solution is

$\displaystyle \Psi = A \left(e^{i\vec{k} \cdot \vec{r}} + f(\theta) \frac{e^{ikr}}{r} \right)$

Here is the outline:

Solve the Schrödinger equation in partial wave
1. the inner part is solved numerically
2. the outer part is asymptotically approach to free plane wave with a phase shift ( or the S-Matrix)
matching the boundary condition
compare the outer part with the general form of the solution:
$\displaystyle \Psi = A \left(e^{i\vec{k} \cdot \vec{r}} + f(\theta) \frac{e^{ikr}}{r} \right)$
that give out the scattering amplitude $f(\theta)$
The elastic cross section is $|f(\theta)|^2$

Another way to solve the Schrödinger equation is solve it in a partial wave,

$\displaystyle \Psi_{lm} = \frac{u_l(kr)}{kr} Y_{lm}(\theta, \phi)$

where the radial solution $u_l$ satisfy,

$\displaystyle \left( \frac{d^2}{d \rho^2} +1 - \frac{l(l+1)}{\rho^2} \right) u_l = U(\rho) u_l$

with $\rho = kr$ , $U(\rho ) = V(r)/E$ . The total solution is

$\displaystyle \Psi = \sum_l C_l Y_{lm}(\theta, \phi) \frac{u_l(kr)}{kr}$

As very far away, (the arrow stands for very long distance)

$\displaystyle u_l(kr) \rightarrow a \hat{j}_l + b \hat{n}_l$

where $\hat{j}_l(r)$ is Riccati-Bessel function. The spherical Bessel-function can be expressed into spherical Hankel function.

$\displaystyle \hat{h}_l^{\pm}(r) = \hat{n}_l(r) \pm i \hat{j}_l(r) \rightarrow e^{\pm i (r - l \frac{\pi}{2}) }$

Follow below manipulation,

$\displaystyle u_l(kr) \rightarrow a \hat{j}_l + b \hat{n}_l = \frac{a}{2i}\left(\hat{h}_l^+ - \hat{h}_l^- \right) + \frac{b}{2}\left(\hat{h}_l^+ + \hat{h}_l^- \right) \\ ~~~~~~~~~~~~~~~~~~~~~~~~~= \frac{a+ib}{2i} \hat{h}_l^+ - \frac{a-ib}{2i} \hat{h}_l^- \\~~~~~~~~~~~~~~~~~~~~~~~~~= \frac{a-ib}{2i} \left( S_l \hat{h}^+ - \hat{h}^- \right), ~~S_l = \frac{a + i b}{ a - i b }$

Drop the $a-ib$ ,

$\displaystyle u_l(kr) \rightarrow \frac{1}{2i} \left( S_l \hat{h}^+ - \hat{h}^- \right) = \frac{1}{2i} \left( S_l e^{i(kr - l \pi/2)} - e^{-i(kr - l\pi/2)} \right)$

The total wave function is the sum of all partial waves,

$\displaystyle \Psi = \sum_l C_l Y_{lm}(\theta, \phi) \frac{u_l(kr)}{kr}$

At far distance, $m = 0$ and $Y_{lm} \rightarrow P_l(cos\theta)$ .

$\displaystyle \Psi \rightarrow \sum_l C_l P_l(\cos\theta)\frac{1}{2i} \frac{ S_l e^{i(kr - l \pi/2)} - e^{-i(kr - l\pi/2)} }{kr}$

Separate out $e^{\pm ikr}$

$\displaystyle \Psi \rightarrow \frac{e^{-ikr}}{r} \sum_l P_l(\cos\theta) \frac{ -C_l e^{i l \pi/2}}{2ik} + \frac{e^{ikr}}{r} \sum_l P_l(\cos\theta) \frac{ C_l S_l e^{- i l \pi/2}}{2ik}$

The plane wave solution can be expanded into partial wave,

$\displaystyle e^{i\vec{k} \cdot \vec{r}} = \sum_l (2l+1) i^l P_l(\cos\theta) \frac{\hat{j}_l(kr)}{kr}$

Thus,

$\displaystyle \hat{j}_l = \frac{\hat{h}_l^+ - \hat{h}_l^-}{2i} \rightarrow \frac{e^{i(r-l\pi/2)} - e^{-i(r-l\pi/2)}}{2i}$

substitute into the plane wave

$\displaystyle e^{i\vec{k} \cdot \vec{r}} \rightarrow \sum_l (2l+1) i^l P_l(\cos\theta) \frac{e^{i(kr-l\pi/2)} - e^{-i(kr-l\pi/2)}}{2ikr}$

substitute into $\Psi =A( e^{i\vec{k} \cdot \vec{r}} + f(\theta) \frac{e^{ikr}}{r} )$ and collect $e^{\pm ikr}$ .

$\displaystyle \Psi \rightarrow A\left(\sum_l (2l+1) i^l P_l(\cos\theta) \frac{e^{i(kr-l\pi/2)} - e^{-i(kr-l\pi/2)}}{2ikr} + f(\theta) \frac{e^{ikr}}{r} \right)$

$\displaystyle \Psi \rightarrow \frac{e^{- ikr}}{r} \sum_l (2l+1) i^l P_l(\cos\theta) \frac{-A e^{il\pi/2}}{2ik} + \\ \frac{e^{ikr}}{r} \left( A f(\theta) + \sum_l (2l+1) i^l P_l(\cos\theta) \frac{Ae^{-il\pi/2}}{2ik}\right)$

Compare the coefficient of $\frac{e^{\pm ikr}}{r}$

$\displaystyle \frac{ -C_l e^{i l \pi/2}}{2ik} = (2l+1) i^l \frac{-A e^{il\pi/2}}{2ik}$

Thus,

$\displaystyle \Rightarrow C_l = A (2l+1) i^l$

And,

$\displaystyle \sum_l P_l(\cos\theta) \frac{ C_l S_l e^{- i l \pi/2}}{2ik} = A f(\theta) + \sum_l (2l+1) i^l P_l(\cos\theta) \frac{Ae^{-il\pi/2}}{2ik}$

$\displaystyle f(\theta) = \sum_l P_l(\cos\theta) (2l+1) i^l (S_l -1) \frac{e^{- i l \pi/2}}{2ik}$

And since $e^{-i l \pi/2} = (-i)^l$ ,

$\displaystyle f(\theta) = \frac{1}{k} \sum_l P_l(\cos\theta) (2l+1)\frac{(S_l -1)}{2i}$

And

$\displaystyle S_l = \frac{a+i b}{a - ib} = e^{2i\delta_l } \Rightarrow \frac{S_l -1 }{2 i} = e^{i \delta_l } \sin(\delta_l )$

Finally,

$\displaystyle f(\theta) = \frac{1}{k} \sum_l P_l(\cos\theta) (2l+1) e^{i \delta_l } \sin(\delta_l )$

The elastics cross section ( because the energy is same as incoming and outgoing wave ) is

$\displaystyle \frac{d\sigma}{d\Omega} = |f(\theta)|^2$

The total elastics cross section is

$\displaystyle \sigma = 2\pi \int |f(\theta)|^2 d(\cos\theta)$

Using

$\displaystyle \int P_l(\cos\theta) P_{l'}(\cos\theta) d(\cos\theta) = \frac{2}{2l+1} \delta_{ll'}$

$\displaystyle \sigma = 2\pi \int |f(\theta)|^2 d(\cos\theta) = \frac{4\pi}{k^2} \sum_l (2l+1) \sin^2(\delta_l )$

At $\theta = 0 , P_l(1) = 1$ , $e^{i \delta_l } = \cos \delta_l + i \sin \delta_l$ ,

$\displaystyle f(0) = \frac{1}{k} \sum_l (2l+1) ( \cos \delta_l \sin\delta_l + i \sin^2\delta_l )$

$\displaystyle \sigma = \frac{4\pi}{k^2} \sum_l (2l+1) \sin^2(\delta_l ) = \frac{4\pi}{k} Im(f(0))$

Coulomb Energy for two protons system

March 5, 2020

GoLuckyRyan Basic Coulomb, Coulomb Energy, spherical Bessel, spherical harmonic Leave a comment

From this post and this post, we already have the framework to calculate the energy. The 1-body intersection is zero, as we do not interested on the kinematics energy. And also, suppose we already have the solution for the wave function for a fixed mean-field. Although it is unrealistic, it is a good starting point.

In fact, we have to solve the Schrödinger equation (using Hartree-Fock or numerical solution) with a realistic NN-interaction with and without Coulomb potential, and subtract the difference to get the Coulomb energy.

We have to evaluate the integral

$\displaystyle W = \left<12||12\right> - \left<12||21\right>$

When the 2 protons are in same orbital, the total wavefunction is

$\displaystyle \Psi(1,2) = \frac{1}{\sqrt{2}}\left(\uparrow \downarrow - \downarrow \uparrow \right) \phi(r_1) \phi(r_2)$

The exchange term is zero as the Coulomb operator does not act on the spin part.

$\displaystyle \left<12||12\right> = \sum_{L}^{\infty} \sum_{M=-L}^{L} \frac{4\pi}{2L+1}\int_{0}^{\infty} \int_{0}^{\infty} \phi^2(x) \phi^2(y) \frac{r_<^L}{r_>^{L+1}} x^2 y^2 dx dy \\ \int Y_{lm}^*(1) Y_{LM}^*(1) Y_{lm}(1) d\Omega_1 \int Y_{lm}^*(2) Y_{LM}(2) Y_{lm}(2) d\Omega_2$

For the angular part, from this post or this post, we have

$\displaystyle \int Y_{l_1m_1}^* Y_{LM}^* Y_{l_2m_2} d\Omega = \sqrt{\frac{(2L+1)(2l_1+1)}{4\pi(2l_2+1)}} C_{L0l_10}^{l_20} C_{LMl_1m_1}^{l_2m_2}$

For $l_1 = l_2$ , the Clebsh-Gordon coefficient dictated that $2l \geq L$ , $L = even$ , $M = 0$ .

We use the potential is 3D spherical infinite well with radius $a = 2.2$ fm. the wave function is spherical Bessel function ,

$\phi_n(r ) = \begin{matrix} A j_n(k r) & r < a \\ 0 & r > a \end{matrix}, j_n(ka) = 0$ ,

where $A$ is the normalization factor.

Screenshot 2020-03-05 at 17.48.19.png

The first non-zero root of spherical Bessel Function take the form,

$r_0(n) = 2.96519 + 0.701921\sqrt{n} + 0.993301 n, n < 500$

The energy for the $n$ -th s-state is

$\displaystyle E_n = \frac{\hbar k^2}{2 m_p^2}, k a = r_0(n)$

In Mathematica, the radial calculation can be done as

a=2.2;
R[r_]:=Piecewise[{{SphericalBesselJ[n, r0[n] r], 0 < r < a}, {0, r > 0}}];
A=(NIntegrate[(R[r])^2 r^2, {r, 0, ∞}]])^(-½);
w = NIntegrate[(R[x]R[y])^2 y^L/x^(L+1) x^2 y^2 , {x, 0, ∞}, {y,0, x}] + NIntegrate[(R[x]R[y])^2 x^L/y^(L+1) x^2 y^2 , {x, 0, ∞}, {y,x, ∞}]
W=1.44 A^4 w

The result for first n s-orbital is

Screenshot 2020-03-05 at 17.59.37.png

I also extract the radius for uniformly charged sphere. The trend is reasonable as higher principle number, the wave function is close to the boundary $a = 2.2 \textrm{fm}$ .

Next time, we will work on Woods-Saxon potential.

Finite Spherical Square Well

November 15, 2017

GoLuckyRyan Basic 2H, deuteron, Hankel, neutron Halo, spherical Bessel Leave a comment

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

In mathematica, the spherical Hankel function is a build in function.

SphericalHankelH1[n, r]

If the potential depth is 60 MeV for proton. Radius is 1 unit for a light nuclei. Set $\hbar = 1, c = 1, m = 1$ . (This is equivalent to radius of $\sqrt{(\hbar c)^2 / m} = 6.44 \textrm{fm}$

The result is follow,

FiniteSquareWell

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,

compare

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

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,

Fermi and Gamow-Teller Transition

February 29, 2016

GoLuckyRyan Basic beta decay, Fermi, Fermi-Golden rule, ft value, Gamow-Teller, half life, log ft, neutron, pion, spherical Bessel, weak force Leave a comment

The beta decay is caused by the weak interaction. The weak interaction is very short range, because the mediate particles, the $W^{\pm}$ and $Z^0$ bosons are 80 GeV and 91 GeV respectively. The effective range is like $10^{-3}$ fm. So, the interaction can assumed to be a delta function and only the coupling constant matter. The Fermi coupling constant is

$1.17 \times 10^{-11} (\hbar c)^2~ \mathrm{MeV^2}$

The fundamental process of beta decay is the decay of quark.

$\displaystyle u \xrightarrow{W^+} d + e^+ + \nu_e$

Since a pion is made from up and down quark, the decay of pion into position and electron neutrino is also due to weak interaction.

The Hamilton of the beta decay is

$\displaystyle H_w(\beta^{\pm})=G_V \tau_{\mp} + G_A \sigma \tau_{\mp}$

where $G_V$ is the vector coupling constant, the term is called Fermi transition. The $\tau_{\pm}$ is the isospin ladder operator. The beta+ decay changes the isospin from +1/2 (neutron) to -1/2 (proton). The $G_A$ is the axial coupling constant, the term is called Gamow-Teller transition. $\sigma$ is spin operator. Because of this operator, the Gamow-Teller transition did not preserve parity.

The $G_A$ is different from $G_V$ , which is caused by the effect of strong interaction. The Goldberger-Trieman relation

$\displaystyle g_A = \frac{G_A}{G_V} = \frac{f_\pi g_{\pi N}}{M_N c^2} = -1.3$

where $f_\pi \sim 93~\textrm{MeV}$ is the pion decay constant. $g_{\pi N} \sim 14 \times 4\pi$ is the coupling constant between pion and nucleon. This, we can see the effect of the strong interaction, in which pion is the meson for strong nuclear force.

The transition probability can be estimated by Fermi-Golden rule

$\displaystyle W(p_e)=\frac{2\pi}{\hbar}|\left< \psi_f|H|\psi_0\right> |^2 \rho(E_f)$

the final state wavefunction

$\displaystyle \left|\psi_f\right> = \frac{1}{\sqrt{V}} e^{ik_e r} \frac{1}{\sqrt{V}} e^{ik_{\nu}r} \left|j_f m_f\right>$

$\displaystyle e^{ikr} = \sum \limits_{L}\sqrt{4\pi (2L+1)} i^L j_L(kr) Y_{L0}(\theta)$

using long wavelength approximation, the spherical Bessel function can be approximated by the first term.

$\displaystyle j_L(kr) \sim \frac{(kr)^L}{(2L+1)!!}$

$\displaystyle \left| \psi_f\right>=\frac{1}{V}(1 + i \sqrt{\frac{4\pi}{3}} Y_{10} + ...) \left|j_f m_f\right>$

The first term 1, or L=0 is called allowed decay, so that the orbital angular momentum of the decayed nucleus unchanged. The higher order term, in which the weak interaction have longer range has very small probability and called L-th forbidden decay.

The density of state is

$\displaystyle \rho(E_f) = \frac{V}{2\pi^2 \hbar^7 c^3} F(Z,E_e)p_e^2 (E_0-E_e) ( (E_0-E_e)^2-(m_{\nu} c^2)^2)^2$

where the $F(Z, E_e)$ is the Fermi function.

The total transition probability is the integration with respect to the electron momentum.

$\displaystyle W = \int W(p_e) dp_e = \frac{m_e^5 c^4}{2 \pi^3 \hbar^7} f(Z,E_0) |M|^2$

where $f(Z,E_0)$ is the Fermi integral. The half-life

$\displaystyle T_{1/2} = \frac{\ln{2}}{W}$

To focus on the beta decay from the interference of the density of state, the ft-value is

$\displaystyle ft = f(Z,E_0) T_{1/2} =\frac{2\pi^3\hbar^7}{m_e^5 c^4} \frac{\ln{2}}{|M|^2}$

The ft-value could be difference by several order.

There is a super-allowed decay from $0^{+} \rightarrow 0^{0}$ with same isospin, which the GT does not involve. an example is

$\displaystyle ^{14}\mathrm{O} \rightarrow ^{14}\mathrm{N} + e^+ + \nu_e$

The ft-value is 3037.7s, the smallest of known.

Fermi	Gamow-Teller
$\Delta S=0$	$\Delta S=1$
$J_f=J_i + L$	$J_f=J_i + L+1$
$T_f=T_i + 1$

transition	L	$\log_{10} ft_{1/2}$	$\Delta J$		$\Delta T$	$\Delta \pi$
transition	L	$\log_{10} ft_{1/2}$	Fermi	GT	$\Delta T$	$\Delta \pi$
Super allowed		3.1 ~ 3.6	$0^+ \rightarrow 0^+$	not exist	0	no
allowed	0	2.9 ~ 10	0	(0), 1	0, 1 ; $T_i=0\rightarrow T_f=0$ forbidden	no
1st forbidden	1	5 ~ 19	(0),1	0, 1, 2	0,1	yes
2nd forbidden	2	10 ~18	(1), 2	2, 3		no
3rd forbidden	3	17 ~ 22	(2), 3	3, 4		yes
4th forbidden	4	22 ~ 24	(3), 4	4, 5		no

The () means not possible if either initial or final state is zero. i.e $1^{-} \rightarrow 0^+$ is not possible for 1st forbidden.

Scattering phase shift

January 10, 2011

GoLuckyRyan advance, Scattering, theory Bessel function, Cross section, phase shift, potential, scattering, spherical Bessel Leave a comment

for a central potential, the angular momentum is a conserved quantity. Thus, we can expand the wave function by the angular momentum wave function:

$\sum a_l Y_{l , m=0} R_l(k, r)$

the m=0 is because the spherical symmetry. the R is the radial part of the wave function. and a is a constant. k is the linear momentum and r is the radial distance.

for free particle, potential equal to zero,

$R_l(k,r) \rightarrow J_{Bessel} (l, kr )$

which is reasonable when r is infinite and the nuclear potential is very short distance. when r goes to infinity,

$J_{Bessel} (l,kr) \rightarrow \frac {1}{kr} sin( k r - \frac{1}{2} l \pi )$

for elastic scattering, the probability of the current density is conserved in each angular wave function, thus,

the effect of the nuclear potential can only change the phase inside the sin function:

$\frac{1}{kr} sin( k r - \frac {1}{2} l \pi +\delta_l )$

with further treatment, the total cross section is proportional to $sin^2(\delta_l)$ .

thus, by knowing the scattering phase shift, we can know the properties of the nuclear potential.

for more detail : check this website

Nuclear Physics 101

Variational method for Woods-Saxon potential

Orthogonality of Spherical Bessel function J

Hankel Transform

Extracting phase shift

S-matrix and scattering amplitude

Coulomb Energy for two protons system

Finite Spherical Square Well

Wave function in momentum space

Fermi and Gamow-Teller Transition

Scattering phase shift

Recent Posts

Calendar

Archives

Keywords

# of Visitors

Login