Orthogonality of Spherical Bessel function J

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

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.

Common functions expressed as Hypergeometric function

Leave a comment

General Hypergeometric function can be expressed in power series

\displaystyle {}_pF_q(a_1, a_2,... a_p ; b_1, b_2, ... b_q; z) = \sum_n \frac{(a_1)_n(a_2)_n ... (a_p)_n}{(b_1)_n (b_2)_n ... (b_q)_n} \frac{z^n}{n!}

where (a)_n is Pochhammar symbol,

\displaystyle (a)_n = \frac{\Gamma(a+n)}{\Gamma(a)} = a(a+1)...(a+n-1)

The General hypergeometric function satisfies the following differential equation,

\displaystyle \frac{d}{dz}[(\theta_{b_1}-1)(\theta_{b_2}-1)... (\theta_{b_q}-1)]y = [\theta_{a_1}\theta_{a_2}...\theta_{a_p}] y

where

\displaystyle \theta_{a} = z\frac{d}{dz} + a


For p = q = 0 , the differential equation becomes

\displaystyle \frac{d}{dz} y = y  \Rightarrow  y = {}_2F_1(;;z) = \exp(z)


For p = 0, q = 1,

\displaystyle \frac{d}{dz}\left( z\frac{d}{dz} + c -1 \right) y = y \Rightarrow \displaystyle z\frac{d^2y}{dz^2} + c \frac{dy}{dz} -y = 0


For p = 1, q = 0

\displaystyle \frac{d}{dz} y= \left(z\frac{d}{dz}+a \right)y \Rightarrow \displaystyle (z-1)\frac{d}{dz} y + ay = 0


For p = 1 = q

\displaystyle \frac{d}{dz}\left( z\frac{d}{dz} + c -1 \right) y = \left(z\frac{d}{dz}+a \right) y \Rightarrow \displaystyle z\frac{d^2y}{dz^2} + (c-z) \frac{dy}{dz} - ay = 0


The Gauss Hypergeometric function is p = 2, q = 1,

\displaystyle {}_2F_1(a,b;c;z) =\sum_n \frac{(a)_n(b)_n}{(c)_n} \frac{z^n}{n!}

which satisfies,

\displaystyle x(1-x) \frac{d^2y}{dx^2} + (c - (a+b+1)x)\frac{dy}{dx} - aby = 0


There are some interesting expression for Pochhammar symbol

\displaystyle (-n)_{k} = (-n)(-n+1)...(-n+k-1) \\ = (-1)^k (n)(n-1)...(n-k+1) \\ = (-1)^k \frac{n!}{(n-k)!}

when k = n

(-n)_n = (-1)^n n!

when k = n + r, r>0

(-n)_{n+r} = 0


Here are list of common function into hypergeometric function

{}_0F_0(; ; z) = e^z

{}_1F_0(-a; -z) = (1+z)^a

\displaystyle {}_0F_1\left(;\frac{1}{2}; -\frac{z^2}{4} \right) = \cos(z)

\displaystyle {}_0F_1\left(;\frac{3}{2}; -\frac{z^2}{4} \right) = \frac{1}{z} \sin(z)

\displaystyle {}_0F_1\left(;a+1; -\frac{z^2}{4} \right) = \frac{2^a}{z^a} \Gamma(a+1) J_a(z)

where J_a(z) is Bessel function of first kind, which satisfies

\displaystyle z^2 \frac{d^2y}{dz^2} + z \frac{dy}{dz} + (z^2 - a^2)y = 0

\displaystyle {}_0F_1\left(; \frac{1}{2}; \frac{z^2}{4} \right) = \cosh(x)

\displaystyle {}_0F_1\left(;\frac{3}{2}; \frac{z^2}{4} \right) = \frac{1}{z} \sinh(z)

\displaystyle {}_0F_1\left(;a+1; \frac{z^2}{4} \right) = \frac{2^a}{z^a} \Gamma(a+1) I_a(z)

where I_a(z) is modified Bessel function of first kind, which satisfies

\displaystyle z^2 \frac{d^2y}{dz^2} + z \frac{dy}{dz} - (z^2 + a^2)y = 0

\displaystyle {}_1F_1\left(\frac{1}{2}; \frac{3}{2}; -z^2 \right) = \frac{\sqrt{\pi}}{2z} Erf(z)

where Erf(z) is error function

Erf(z) = \int_0^z \exp(-t^2) dt

\displaystyle {}_2F_1\left(-a,a; \frac{1}{2}; \sin^2(z) \right) = \cos(2az)

\displaystyle {}_2F_1\left(\frac{1}{2}+a, \frac{1}{2}-a; \frac{3}{2}; \sin^2(z) \right) = \frac{\sin(2az)}{2a \sin(z)}

\displaystyle {}_2F_1(1,1;2;-z) = \frac{1}{z} \log_e(z+1)

\displaystyle {}_2F_1(\frac{1}{2},-1;\frac{a}{2};z) = 1-  \frac{z}{a}

\displaystyle {}_2F_1\left( \frac{1}{2}, 1; \frac{3}{2}; z^2 \right) = \frac{1}{2z} \log_e \left( \frac{1+z}{1-z} \right) = \frac{1}{z} \tanh^{-1}(z)

\displaystyle {}_2F_1 \left( \frac{1}{2}, 1; \frac{3}{2} ; -z^2 \right) = \frac{1}{z} \tan^{-1}(z)

\displaystyle {}_2F_1 \left( \frac{1}{2}, \frac{1}{2}; \frac{3}{2} ; z^2 \right) = \frac{1}{z}\sin^{-1}(z)

\displaystyle {}_2F_1 \left( \frac{1}{2}, \frac{1}{2}; \frac{3}{2} ; -z^2 \right) = \frac{1}{z}\sinh^{-1}(z)

\displaystyle {}_2F_1 \left( \frac{1}{2}, \frac{1}{2}; \frac{3}{2} ; \frac{1-z}{2} \right) = \frac{1}{\sqrt{2(1-z)}}\cos^{-1}(z)

\displaystyle {}_2F_1\left(-n, n+1; 1; \frac{1-z}{2} \right) = P_n(z)

where P_n(z) is Legendre function, which satisfies

\displaystyle (1-z^2)\frac{d^2y}{dz^2} -2z \frac{dy}{dz} + n(n+1) y = 0

\displaystyle {}_2F_1\left(m-n,m+n+1; m+1; \frac{1-z}{2} \right) \\= (-1)^m\frac{(n-m)!m!2^m}{(n+m)!(1-x^2)^{\frac{m}{2}}} P_n^m(z), m\geq0

where P_n^m(z) is associate Legendre function, which satisfies

\displaystyle (1-z^2)\frac{d^2y}{dz^2} -2z \frac{dy}{dz} + \left(n(n+1) -\frac{m^2}{1-z^2} \right)y = 0

\displaystyle {}_2F_1\left(\frac{1}{2}, \frac{1}{2}; 1; z^2 \right) = \frac{2}{\pi} K(z)

where K(z) is complete elliptic integral of 1st kind

\displaystyle K(z) = \int_0^{\frac{\pi}{2}} \frac{1}{\sqrt{1-z^2 \sin^2(t)}} dt

\displaystyle {}_2F_1\left(-\frac{1}{2}, \frac{1}{2}; 1; z^2 \right) = \frac{2}{\pi} E(z)

and E(z) is complete elliptic integral of 2nd kind

\displaystyle E(z) = \int_0^{\frac{\pi}{2}} \sqrt{1-z^2 \sin^2(t)} dt

Reference

“Notes on hypergeometric functions” by John D. Cook (April 10, 2003)
“Generalized Hypergeometric Series” by W. N. Bailey, Cambridge (1935)
“Handbook of Mathematical Functions” by Abramowitz and Stegun (1964)
“The special functions and their approximations” by Yudell L. Luke v. 1 (1969)
“Concrete Mathematics” by Graham, Knuth, and Patashnik (1994)


In Wolfram research (http://functions.wolfram.com/functions.html), many functions are listed. We can click to a function, then we click “Representations through more general functions”, then “Through hypergeometric functions”, then we can see how the function looks like.

Very short introduction to Partial-wave expansion of scattering wave function

Leave a comment

In a scattering problem, the main objective is solving the Schrödinger equation

H\psi=(K+V)\psi=E\psi

where H is the total Hamiltonian of the scattering system in the center of momentum, K is the kinetic energy and V is the potential energy. We seek for a solution \psi,

\displaystyle \psi_{k}^{+}(r)=e^{i\vec{k}\cdot \vec{r}}+f(\theta)\frac{e^{ikr}}{kr}

The solution can be decomposed

\displaystyle \psi_{k}^{+}(r)=R_{l}(k,r)Y_{lm}(\theta,\phi)=\frac{u_{l}(k,r)}{kr}Y_{lm}(\theta,\phi)

The solution of u_{l}(k,r) can be solve by Runge-Kutta method on the pdf

\displaystyle \left(\frac{d^2}{d\rho^2} + 1 - \frac{l(l+1)}{\rho^2} \right)u_{l}(k,\rho)=U(\rho)u_{l}(k,\rho)

where \rho=kr, k=\sqrt{2\mu E}/\hbar, \mu=(m_1+m_2)/(m_1 m_2) and U=V/E.


For U = 0, the solution of u_l is

\displaystyle u_{l}(k,r)=\hat{j}_l(\rho) \xrightarrow{r\rightarrow \infty} \sin(r') = \frac{e^{ir'}-e^{-ir'}}{2i}

where r' = kr-l\pi/2 and \hat{j}_l is the Riccati-Bessel function. The free wave function is

\displaystyle \phi_k(r)=e^{i\vec{k}\cdot\vec{r}}=\sum\limits_{l=0} P_l(\cos(\theta)) \frac{2l+1}{2ikr}i^l (e^{ir'}-e^{-ir'})

where P_l(x) is the Legendre polynomial.

Note that, if we have Coulomb potential, we need to use the Coulomb wave instead of free wave, because the range of coulomb force is infinity.


For U\neq 0, the solution of u_l(r<R) can be found by Runge-Kutta method, where R is a sufficiency large that the potential V is effectively equal to 0.  The solution of u_l(r>R) is shifted

\displaystyle u_{l}(k,r>R)=\hat{j}_l(\rho)+\beta_l \hat{n}_l(\rho) \xrightarrow{r\rightarrow \infty} \frac{1}{2i}(S_l e^{ir'}-e^{-ir'})

where S_l is the scattering matrix element, it is obtained by solving the boundary condition at r = R. The scattered wave function is

\displaystyle \psi_k(r)=\sum\limits_{l=0} P_l(\cos(\theta)) (2l+1) i^l \frac{u_l(r)}{kr}

put the scattered wave function and the free wave function back to the seeking solution, we have the f(\theta)

 \displaystyle f(\theta) = \sum\limits_{l=0} P_l(\cos(\theta)) \frac{2l+1}{2ik} (S_l - 1)

and the differential cross section

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


In this very brief introduction, we can see

  • How the scattering matrix S_l is obtained
  • How the scattering amplitude f(\theta) relates to the scattering matrix

But what is scattering matrix? Although the page did not explained very well, especially how to use it.

3-D spherical infinite well

Leave a comment

the potential is

\displaystyle V(r,\theta,\phi) = \begin{pmatrix} 0 & |r|<a \\ \infty & |r| \geq a \end{pmatrix}

The Laplacian in spherical coordinate is:

\displaystyle \nabla^2 = \frac{1}{r^2}\left( \frac{d}{dr}r^2\frac{d}{dr} \right)-\frac{L^2}{r^2}= \frac{d^2}{dr^2} + \frac{2}{r}\frac{d}{dr} - \frac{L^2}{r^2}

since the L is the reduced angular momentum operator, if we set the solution be:

\displaystyle \psi(r,\theta,\phi) = R(r) Y_{lm}(\theta,\phi)

Then the angular part was solved and the radial part becomes:

\displaystyle L^2 Y_{lm} = l(l+1) Y_{ml}

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

\displaystyle k^2 = \frac{2 m E}{\hbar^2}

The radial equation is the spherical Bessel function.

The solution was common written as:

\displaystyle R(r) = j_l( k r) = \left( - \frac{r}{k} \right)^l \left(\frac{1}{ r} \frac{d}{dr}\right)^l \frac{sin(k r)}{kr}

The Boundary condition fixed the k and then the energy,

\displaystyle j_l ( k_{nl} a ) = 0

the all possible root are notated as n. thus the quantum numbers for this system are:

  • n , the order of root
  • l , the angular momentum

We can see in here, the different between Coulomb potential and spherical infinite well:

  • there is no restriction on n and l, therefore, there will be 1s, 1p, 1d, 1f orbit.
  • the energy level also depend on angular momentum, since it determined the order of spherical Bessel function.

we can realized the energy level by the graph of Bessel function. we set some constants be 1, the root are :

k_{nl} a = \displaystyle \frac{1}{\hbar} \sqrt{ 2 m a^2} \sqrt{E_{nl}} = \pi \sqrt{E_{nl}}

Thus, we plot

\displaystyle j_l( \pi \sqrt{E_{nl}})

Since k_{nl} is a scaling factor to “force” the function to be zero at the boundary. Interestingly, the spherical Bessel function is not normalizable or orthogonal with r^2 , i.e.

\displaystyle \int_{0}^{\infty} j_l(r) j_{l'}(r) r^2 dr

is diverged. Of course, the spherical Bessel function is a “spherical wave” that propagating in space, same as plane wave, which is also not normalizable or orthogonal. However, in the infinite spherical well, with the scaling factor k for same l , they are orthogonal! And for difference l , the spherical part are already orthogonal. Thus, the wave function for difference energy are orthogonal. i.e. the eigen states are orthogonal! Power of math! In the following plots, the left is l = 0 and the right is l = 1.

pic1.PNG

In fact, even for Woods-Saxon potential, the numerical solution of the eigen wave-function are also orthogonal for same angular momentum.

Scattering phase shift

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