Calculate phase shift of square well

September 7, 2021

GoLuckyRyan Basic , , , 4 Comments

In the book, Introduction to Nuclear Reaction, by C. A. Bertulani, Chapter 2, there is a section that shows an example of phase shift from a square well. I follow the receipt and able to reproduce the result.

The actual method I did is little different. I use Mathematica to solve the radial Schrodinger equation in u(r) = R(r)/r , the differential equation I solve is

\displaystyle \frac{d^2}{dr^2} u(r) = \begin{cases} 0 & r = 0 \\ \frac{2m}{\hbar^2} W(r) u(r) + \frac{l(l+1)}{r^2} u(r) & r > 0 \end{cases}

where W(r) is the Woods-Saxon function, set the diffusiveness parameter to 0.001, the Woods-Saxon shape becomes square well.

In solving the differential equation, I set the right-hand side to be zero when r=0. That solved to problem when l > 0, it has to deal with the 1/r^2 . In fact, since u(0) = 0 , the right-hand side is always zero when r = 0.

For the solving range, I set

\displaystyle r_{max} = 5 \frac{2\pi}{k} + \frac{\sqrt{l(l+1)}}{k}

This is because the “first” peak of the Riccati-Bessel function appear after \sqrt{l(l+1)}/k . This will make sure the range has 5 oscillations in the solution.

Next, I table the solution for the last 2 oscillations, find the maximum, normalize the solution, and then fit with

\displaystyle u(r) \xrightarrow{r \rightarrow \infty} A \hat{j}_l(kr) + B \hat{n}_l(kr)

and extract A, and B.

The phase shift is

\displaystyle \tan(\phi) = \frac{A}{-B} .

The negative is due to the \hat{n}_l is negatively defined in Mathematica.

The elastics cross section is

\displaystyle \sigma_{el} = \frac{4\pi}{k^2} \sum_{l} \sqrt{2l+1} \sin^2(\phi_l)

From the calculation, we found that, at very small energy, only s-wave scattering cross section is non-zero, which is agree with theory.

[20220531] update

The Mathematical notebook can be downloaded here.

Extracting phase shift

March 5, 2021

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

Coulomb wave function

April 17, 2020

GoLuckyRyan Basic , , , , , , Leave a comment

The Coulomb wave function is the solution of a REPULSIVE Coulomb potential. The Schrödinger equation is

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

where e^2 = 1.44~\textrm{MeV.fm}, separate the radial and angular part as usual,

\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)

setting

\displaystyle k^2 = \frac{2mE}{\hbar^2},  \eta =\frac{k}{2E} Q_1Q_2 e^2

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

using the usual substitution R = u / r

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

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

Setting x = kr

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

The short range behaviour is approximated as

\displaystyle  \frac{d^2u(x)}{dx^2} = \frac{l(l+1)}{x^2}u(x)  \Rightarrow u(x) \approx x^{l+1}~\textrm{or}~ x^{-l}

and for long range behaviour, it should approach Riccati–Bessel functions \hat{j}_l, \hat{n}_l with phase shift.

Set u(x) = x^{l+1} \exp(i x) y(x) , the equation of u(x) becomes,

\displaystyle  x \frac{d^2y}{dx^2} + (2L+2 + 2ix) \frac{dy}{dx} + (2i(L+1) - 2\eta) y = 0

Change of variable z = - 2i x ,

\displaystyle  z \frac{d^2y}{dz^2} + (2L+2 - z ) \frac{dy}{dx} - (L+1 + i \eta) y = 0

This is our friend Laguerre polynomial again!!! with \alpha = 2L+1,  n = -(L+1+i\eta) . Since the n is not an integer anymore, we go to a more general case, that is the Kummer’s equation,

\displaystyle z \frac{d^2 w}{dz} + (b-z) \frac{dw}{dz} - a w = 0

The solution of Kummer’s equation is the confluent hypergeometric function

w(z) = _1F_1(a, b, z)

Thus, the solution for the radial function is

\displaystyle u_l(x) =A x^{l+1} e^{ix} _1F_1(L+1+i \eta, 2L+2, -2 i x )

where A is a normalization constant by compare the long range behaviour with Riccati-Bessel function. The full solution is,

\displaystyle u_l(x) =  \\ F_l(x) = \frac{2^l e^{-\pi \eta/2} |\Gamma(l+1+i\eta)| }{(2l+1)!} x^{l+1} e^{ix}~_1F_1(L+1+i \eta, 2L+2, -2 i x )

At long range,

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

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

Using Kummer’s transform

\displaystyle _1F_1(a,b,z) = e^z ~_1F_1(b-a, b, -z)

The solution can be written as

\displaystyle u_l(x) =  \\ F_l(x) = \frac{2^l e^{-\pi \eta/2} |\Gamma(l+1+i\eta)| }{(2l+1)!} x^{l+1} e^{-ix}~_1F_1(L+1-i \eta, 2L+2, 2 i x )

Another solution should behave like \cos and unbound at x = 0 .

Set u(x) = x^{-l} \exp(i x) y(x) , the equation of u(x) becomes,

\displaystyle  x \frac{d^2y}{dx^2} + (-2L + 2ix) \frac{dy}{dx} + (-2i L - 2\eta) y = 0

Change of variable z = -2ix ,

\displaystyle  z \frac{d^2y}{dz^2} + (-2L -z ) \frac{dy}{dx} - ( - L + i \eta) y = 0

\displaystyle u_l(x) =A x^{-l} e^{ix}~ _1F_1(-L+i \eta, -2L, -2 i x )

However, _1F_1(a, b, z) does not exist for non-position b.

We can also transform into Whittaker’s equation by change of variable z = 2ix

\displaystyle  -4 \frac{d^2u(x)}{dz^2} + \left( 1 - \frac{-4l(l+1)}{z^2} - \frac{4i\eta}{z} \right) u(x) =0

\displaystyle  \frac{d^2u(x)}{dz^2} + \left( -\frac{1}{4} +\frac{i\eta}{z}- \frac{l(l+1)}{z^2} \right) u(x) =0

using \displaystyle l(l+1) = (l+1/2)^2 - 1/4

\displaystyle  \frac{d^2u(x)}{dz^2} + \left( -\frac{1}{4} +\frac{i\eta}{z} +  \frac{ 1/4 - (l+1/2)^2}{z^2} \right) u(x) =0

This is the Whittaker’s equation with \kappa = i \eta,  \mu = l+1/2

The solutions are

u_l(x) = e^{-ix} (2ix)^{l+1} ~_1F_1(l+1-i \eta, 2l+2, 2ix)

u_l(x) = e^{-ix} (2ix)^{l+1} U(l+1-i \eta, 2l+2, 2ix)

I still cannot get the second solution, the G_l . According to Wolfram, https://mathworld.wolfram.com/CoulombWaveFunction.html

\displaystyle G_l(x) = \frac{2}{\eta C_0^2(\eta)} F_l(x) \left( \log(2x) + \frac{q_l(\eta)}{p_l(\eta)} \right) + \frac{x^{-l}}{(2l+1) C_l(\eta)} \sum_{K=-l}^\infty a_k^l(\eta) x^{K+l} ,

where q_l, p_l, a_k^l are defined inAbramowitz, M. and Stegun, I. A. (Eds.). “Coulomb Wave Functions.” Ch. 14 in Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, 9th printing. New York: Dover, pp. 537-544, 1972.

S-matrix and scattering amplitude

April 16, 2020

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

  1. 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)
  2. matching the boundary condition
  3. 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)
  4. 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))Â