Calculate phase shift of square well

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

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 (II)

5 Comments

In the previous post, we tried to derived to Coulomb wave function, and the regular Coulomb wave function F_l(x) is derived, except for the normalization constant, and the mysterious Coulomb phase shift.

From this arxiv article, the Coulomb “Hankel” function is

\displaystyle H_L^{pm}(x) = D_L^{\pm} x^{L+1} e^{\pm i x} U(L+1 \pm i \eta, 2L+2, - \pm 2 i x)

\displaystyle D_L^{\pm} = (-\pm 2i)^{2L+1} \frac{\Gamma(L+1\pm i \eta)}{ C_L \Gamma(2L+2)}

\displaystyle C_L = 2^L \frac{\sqrt{\Gamma(L+1+ i \eta) \Gamma(L+1- i \eta)}}{ e^{\eta \pi/2}\Gamma(2L+2)}

where U(a,b,z) is the confluent hypergeometric function of the second kind.

The regular Coulomb wave function is

\displaystyle F_L(x) = C_L x^{L+1} e^{\pm i x} ~_1F_1(L+1 \pm i \eta, 2L+2, - \pm 2 i x) = \frac{1}{2i} \left( H_L^+ - H_L^- \right)

I am fail to prove the last equality. The irregular Coulomb wave function is

\displaystyle G_L(x) = \frac{1}{2} \left( H_L^+ + H_L^- \right)

The form of the coefficience C_L is different from the form give in this post, but they are the same.

\displaystyle \frac{\sqrt{\Gamma(L+1+ i \eta) \Gamma(L+1- i \eta)}}{ \Gamma(2L+2)} = \frac{\left| \Gamma(L+1+i \eta) \right|}{(2L+1)!}


In Mathematica 11.2, the U(a,b,z) is fail to evaluate when \eta > 2 for l = 0 , and \eta > 0.1 for l = 5 .

Update 20200420, settign WorkingPrecision to 45 can solve the problem. Here are some plots for the Coulomb wave functions with \eta = 10 .

Annotation 2020-04-20 223816.png

Annotation 2020-04-20 223849.png

Annotation 2020-04-20 224003.png

When \eta getting larger, the wave function pushed further.

Annotation 2020-04-20 224017.png

Coulomb wave function

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.

Meaning of phase shift

Leave a comment

For hard sphere with radius R, the phase shift is

\displaystyle \tan \delta_l = - \frac{\hat{j}_l(kR)}{ \hat{n}_l(kR)}

for l = 0 , \displaystyle \delta_0 = - k R

In general, the phase shift has these properties

  • The phase shift is use to calculate the cross section.
  • At very large energy, the phase shift approach 0, \lim_{k\rightarrow \infty} \delta_l = 0
  • When \delta_l = (n+1) \frac{\pi}{2} , it is where resonance take place.
  • At resonance, the slope of the \delta_l(E) is related to the resonance width.
  • The absolute phase shift is defined at \lim_{k\rightarrow \infty} \delta_l = 0 , thus, the absolute phase shift is always positive.
  • The absolute phase shift at k = 0 is equal to n \pi , where n is the number of bound state. ( this is also known as Levinson’s theorem )
  • When \delta_l = n \pi , S_l = 1 and the scattering amplitude for that partial wave is vanished.

Suppose the outer wave take the form

\displaystyle u_l^{II}(k,r) = a_l ( \hat{j}_l(kr) + b_l \hat{n}_l(kr) )

The logarithmic derivative at r = R is

\displaystyle L^{II} = R \frac{1}{u_l^{II}(kR)} \left. \frac{du_l^{II}(kr)}{dr} \right|_{r=R} = kR \frac{\hat{j}_l'(kR) + b_l \hat{n}_l'(kR)}{\hat{j}_l(kR) + b_l \hat{n}_l(kR)}

The inner function will be calculated by numerical method, and the logarithmic derivative at r = R ,

\displaystyle L_{I} = R \frac{(u_l^I)'}{u_l^I}

matching the boundary condition,

\displaystyle L^{I} = L^{II} = kR \frac{\hat{j}_l'(kR) + b_l \hat{n}_l'(kR)}{\hat{j}_l(kR) + b_l \hat{n}_l(kR)}

solve for b_l

\displaystyle b_l = - \frac{kR \hat{j}_l'(kR) - L^I \hat{j}_l(kR)}{kR \hat{n}_l'(kR) - L^I \hat{n}_l(kR)} = \tan (\delta_l)

If Coulomb potential included, simply replace \hat{j}_l \rightarrow F_l and \hat{n}_l \rightarrow G_l .


In more general, the phase shift could indicate the attractive/repulsive of the potential.

  • positive phase shift = attractive potential
  • negative phase shift = repulsive potential

When the potential is attractive ( repulsive), the wave-number k will be larger (smaller) than free plane wave inside the potential, thus, the wave will “advance” more (less), resulting positive (negative) phase shift. For example, I use a Woods-Saxon potential with depth only -1 MeV. A proton s-wave radial wave function is

Annotation 2020-04-16 164023.png

The blue curve is the radial wave function, the orange curve is \hat{j}_0 . The red curve is the potential / 10, and the green line is the energy / 100.

The phase shift is 7 deg.

If I flip the potential to be positive

Annotation 2020-04-16 164257.png

The phase shift is -7 deg.


However, there is an ambiguity of the phase shift, because -7 deg = +  353 deg. Thus, the phase shift need to be measure at various incident energies and infer from the data. For example, using high incident energy so that the wave length inside the potential will not be different from outside. That will give a general ideal on the attractive/repulsive nature of the potential. However, high energy also means the s-wave component is small.

But the S-matrix does not have ambiguity.

S-matrix and scattering amplitude

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

Phase shift of elastics scattering

Leave a comment

I found that the derivation of most “google result” is not clear enough. So here is my derivation. Before process, people may need to review the pervious post.

Most people start from,

u_l = A( \hat{h}_l^- - S_l \hat{h}_l^+ )

that annoying me because it is somehow not “natural”. Why there is a “minus” sign? Why the \hat{h}_l^- is the first term? For my self, a more natural way is,

u_l = a \hat{h}_l^+ + b \hat{h}_l^-

where a, b are complex numbers, but that is still not so natural, because in numerical calculation, for simplicity, there is no complex number, we only have,

u_l = \alpha \hat{j}_l + \beta \hat{n}_l

The first term is alway there as it is the free solution and bounded at r = 0. the second term is caused by the potential.


The goal is to find a solution take the form

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

where the first term is free wave and the second term is scattered wave. The solution for elastics scattering is

\displaystyle \psi = \sum C_l P_l (\cos\theta) \frac{u_l}{kr} \rightarrow \sum C_l P_l(\cos\theta) (\alpha \hat{j}_l + \beta \hat{n}_l)

we used the substitution,

\displaystyle R_l(r) = \frac{u_l(\rho)}{\rho},  \rho = kr.

The radial function can be solved using Rungu-Kutta method on the equation,

\displaystyle \frac{d^2}{d\rho^2} u_l = \frac{2 m_\mu}{\hbar^2} (V-E) u_l + \frac{l(l+1)}{\rho^2}

and the solution of u_l at far away is,

u_l \rightarrow  \alpha \hat{j}_l + \beta \hat{n}_l.

the arrow means r \rightarrow \infty. So, the problem is how to rewrite the solution. In the way, we will see how the phase shift or the S-matrix was found.


The free solution is the spherical wave,

\displaystyle e^{i \vec{k} \cdot \vec{r}} = \sum_l (2l+1) i^l P_l(\cos\theta) j_l(kr)

The spherical Bessel function j_l(kr) cna be express as Heankel function

h_l^{\pm} = n_l \pm i j_l \rightarrow e^{\pm i (kr - l \frac{\pi}{2})}

The + sign is outgoing wave.


\displaystyle u_l \rightarrow (\alpha \hat{j}_l + \beta \hat{n}_l)

\displaystyle = \frac{\alpha}{2i} (\hat{h}_l^{+} - \hat{h}_l^{-}) + \frac{\beta}{2}(\hat{h}_l^{+} + \hat{h}_l^{-})

\displaystyle = \frac{\alpha - i \beta}{2i} (\hat{h}_l^{+} - \hat{h}_l^{-}) + \beta \hat{h}_l^{+}

\displaystyle = (\alpha - i \beta ) \left( \frac{\hat{h}_l^+ - \hat{h}_l^-}{2i} + \frac{\beta}{\alpha - i \beta} \hat{h}_l^+\right)

\displaystyle = (\alpha - i \beta ) \left( \hat{j}_l + \frac{\beta}{\alpha - i \beta} \hat{h}_l^+\right)

Since the u_l should be normalized, we can se \alpha = \cos \delta and \beta = \sin\delta.

\displaystyle \frac{\beta}{\alpha - i \beta } = \sin(\delta) e^{i\delta}

We put u_l back

\displaystyle \psi \rightarrow \sum_l C_l P_l (cos\theta)(\alpha - i \beta ) \left( \hat{j}_l + \sin(\delta) e^{i\delta} \hat{h}_l^+\right)

By setting

\displaystyle C_l = A i^l \frac{2l+1}{\alpha - i \beta} ,

we have the first term is the free wave function. In the second term, \hat{h}_l^+ \rightarrow e^{i(kr - l \frac{\pi}{2}}) / kr . Notice that

e^{i l \frac{\pi}{2}} = i^{-l}

That cancel the i^l term in C_l . And we have

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


some people will write the u_l as \hat{h}_l^{\pm} and the S-matrix,

\displaystyle u_l = \frac{\alpha + i \beta} {2i} \hat{h}_l^+ - \frac{\alpha - i \beta}{2i} \hat{h}_l^-

\displaystyle = -\frac{\alpha - i \beta}{2i} \left( \hat{h}_l^- - \frac{\alpha + i \beta}{\alpha - i \beta} \hat{h}_l^+ \right)

\displaystyle = A' (\hat{h}_l^- - S_l \hat{h}_l^+)

where

\displaystyle S_l =\frac{\alpha + i \beta}{\alpha - i \beta} = e^{2i\delta} .

Remember that this is the S-matrix for elastics scattering.

 

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.

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