Variational method for Woods-Saxon potential

August 31, 2021

GoLuckyRyan Basic , , , , 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

The ground state radial function.

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

Asymptotic normalization coefficient

May 2, 2021

GoLuckyRyan Basic , , , , Leave a comment

The Asymptotic normalization coefficient or the ANC is thought to be an alternative to the spectroscopic factor. As the name suggested, this is the coefficient for asymptote of the radial wave function.

For a single nucleon adding/removal reaction, the reaction probability is

\displaystyle T^2 = \left|\left< \phi_b \Psi_B| V | \phi_a \Psi_A \right>\right|^2

Assume the interaction V only acts on the nucleon that being add or remove, thus

\displaystyle \left< \Psi_B | \Psi_A \right> = \phi(r)

which is the bound state wave function. Due to the nuclear interaction, the bound state wave function most probably not normalized. In other point of view, suppose the nucleus B = A + 1, the wave function of nucleus B could be

\displaystyle \left| \Psi_B \right> = \beta_{00} [\phi_{0} \times \Psi_A(0) ]_{J_B} + ... + \beta_{ij} [\phi_{i} \times \Psi_A(j) ]_{J_B}

where \phi_i is an orbital, \Psi_A(j) is the j-th excited state of nucleus A, and the square bracket is the angular coupling, antisymmetric, and normalization operator. As the wave function of nucleus B must be normalized, Thus,

\displaystyle \sum_{ij} \beta_{ij}^2 = 1 .

And the bound state can be approximated

\displaystyle \phi(r) \approx \beta_{00} \phi_{0} .

The spectroscopic factor of orbital \phi_0 in this A + 1 = B reaction is

\displaystyle \beta_{00}^2 = \int r^2 \phi^2(r) dr

At far away distance, the nucleus potential is very weak or effectively zero. The Schrodinger equation becomes

\displaystyle \left( - \frac{\hbar^2}{2m} \nabla^2 + \frac{Q_1Q_2e^2}{r} \right) \Phi = E \Phi

where e^2 = 1.44~\textrm{MeV.fm}. The Coulomb potential is still here because it is a long range force. Separate the radial and angular 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} + \frac{Q_1Q_2 e^2}{r} \right) R(r) = E R(r)

This radial equation was solved before in this post. There are two solutions for u(r) = R(r)/r , one is bound (in mathematics) and one is unbound. I state the mathematically-bounded solution in here

\displaystyle F_l(r, \eta) = \frac{2^l e^{-\pi \eta /2} \left| \Gamma(l+1+i\eta) \right|}{(2l+1)!} (kr)^{l+1} e^{ikr} {_1F_1(l+1+i\eta, 2l+2, -2ikr)},

\displaystyle k^2 = \frac{2mE}{\hbar^2}

\displaystyle \eta = \frac{k}{2E} Q_1 Q_2 e^2 = \frac{1}{\hbar} \sqrt{\frac{m}{2E}} Q_1 Q_2 e^2

There is another Coulomb wave function G_l(r) , although there is no simple form of the function, the two Coulomb wave function, one behave like a sine wave (F_l(r) ) another one behave like a cosine wave (G_l(r)), they can be combined into a complex function, the Coulomb Hankel function, (it is different from the Hankel function of the first and second kind )

\displaystyle H_l^{\pm}(x, \eta) = D_l^{\pm} x^{l+1} e^{\pm i x} U(l+1\pm \eta, 2l+2, \mp 2ix)

\displaystyle D_l^{\pm} = (\mp 2i) ^{2l+1} \frac{\Gamma(l+1\pm i \eta)} {C_l }

\displaystyle C_l = \frac{2^l}{2^{\eta \pi /2}} \sqrt{\Gamma(l+1+i\eta) \Gamma(l+1-i\eta)}

where U(a,b,z) is the confluent hypergeometric function of the second kind. In Mathematica, the function is given is built-in

HypergeometrixU[l+1+i eta, 2l+2, -2 i r]

And

\displaystyle F_l(r, \eta) = \frac{1}{2i} ( H_l^+(kr, \eta) - H_l^-(kr, \eta) )

\displaystyle G_l(r, \eta) = \frac{1}{2} ( H_l^+(kr, \eta) + H_l^-(kr, \eta) )

We can see, the analogy

\displaystyle F_l(x, \eta) \sim \sin(x) = \frac{1}{2i} ( e^{ix} - e^{-ix} ) \\ G_l(x, \eta) \sim \cos(x) = \frac{1}{2} ( e^{ix} + e^{-ix} )

The long range (unbound, scattering) behaviour of the Coulomb wave function is

\displaystyle F_l(r\rightarrow \infty, \eta) = \sin \left( k r - l \frac{\pi}{2} - \eta \log(2kr) + \sigma_l \right)

\displaystyle G_l(r\rightarrow \infty, \eta) = \cos \left( k r - l \frac{\pi}{2} - \eta \log(2kr) + \sigma_l \right)

where \sigma_l = \arg(\Gamma(l+1+i\eta)) is the Coulomb phase shift.

The ANC C is the coefficient between the bound state wave function and the Coulomb wave function at r \rightarrow \infty . Since it is the bound state long range behaviour, k = i \kappa , we have to use the Coulomb Hankel function.

For neutron, \eta = 0

\displaystyle H_0^\pm(x, 0) = e^{\pm i x}

We can see that, when x = i \kappa r , the $latex H_0^\pm (i \kappa r, 0 ) = e^{-\kappa r} is a bounded and real solution.

This is the same for higher l that the decay factor e^{-\kappa r } appeared in the Hankel function.

In the following, we use the neutron and simplify the Coulomb wave function for few l . When no charge, \eta = 0 . And since it is a bound state,

In fact, the solution for \eta = 0 is the Spherical Bessel function J. As the k = i \kappa is pure imaginary, and we know that the Spherical Bessel function J and Y are unbound for pure imaginary position.

The solution is the Spherical Hankel function h_l^\pm. This is the Coulomb Hankel function divided by the radius. In Mathematica, the built-in function is

SphericalHankelH1[l, i kappa r]

In fact,

\displaystyle h_l^+(x) = \frac{-1}{x} H_l^+(x, 0) \\ h_l^-(x) = \frac{(-1)^l}{x} H_l^-(x,0)

In this paper N.K. Timofeyuk, PRC 88, 044315 (2013), in equation 3, the Whittaker function should be the Whittaker W-function. In Mathematica

WhittakerW[ - i eta, l + 1/2, 2 kappa r ]

The difference between the Whittaker W-function and the Hankel function with complex argument is the normalization factor

\displaystyle W_{-i \eta, l+1/2}(2 \kappa r ) = (2 \kappa r)^{l+1} e^{-\kappa r} U(l+1+i \eta, 2l+2, 2 \kappa r)

Since the ANC is the proportional factor between the bound state wave function and the Whittaker W-function (or the Coulomb wave function, or the spherical Hankel function ), so the normalization factor is important.

As we know that the neutron should behave as Spherical Hankel function, thus, we checked that

\displaystyle W_{0, l+1/2}(2 \kappa r) = (-1) i^l ( \kappa r ) h_{l}^+(i \kappa r) = \sqrt{\frac{2 \kappa r}{\pi}} K_{l+1/2}(\kappa r)

where K_{\alpha} is the modified Bessel function of the 2nd kind. Here are a list of the Whittaker W-function for \eta = 0

\displaystyle W_{0, 1/2}( 2 x) = e^{-x}

\displaystyle W_{0, 3/2}( 2 x) = \frac{e^{-x}}{x} (1+x)

\displaystyle W_{0, 5/2}( 2 x) = \frac{e^{-x}}{x^2} (3 + 3x+x^2)

\displaystyle W_{0, 7/2}( 2 x) = \frac{e^{-x}}{x^3} (15 + 15x+ 6x^2 + x^3)

\displaystyle W_{0, 9/2}( 2 x) = \frac{e^{-x}}{x^4} (105 + 105x+ 45x^2 + 10x^3 + x^4)

\displaystyle W_{0, 11/2}( 2 x) = \frac{e^{-x}}{x^5} (945 + 945x+ 420x^2 + 105x^3 + 15x^4+x^5)

For me, the problem of the ANC is that, the Coulomb Hankel function or the Whittaker W-function cannot be normalized. And without a proper normalization, how can we compare the magnitude of two wave functions??

For example, I calculated the 1s1/2 bound state wave function from a Woods-Saxon potential for neutron. The parameters are V_0 = -50.4, R_0 = 3.1246, a_0 = 0.67, V_{SO} = 19.6, R_{SO} = 3.0238, a_{SO} = 0.66, m_\mu = 884.297, the energy is -2.227 MeV and \kappa = 0.3181 . And here is the comparison

\displaystyle \phi(r) = C \frac{W_{0, 1/2}(2 \kappa r )}{\kappa r}

The Woods-Saxon bound state is supposed to be a pure wave function that has SF = 1 or ANC = 1. But since the Whittaker cannot be normalized, and depends on the “definition”, where the \kappa appears in the denominator or not, the ANC can be different. And for this example the ANC = 0.7, for a pure state.

Finite Spherical Square Well

November 15, 2017

GoLuckyRyan Basic , , , , 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.