Two-body Coulomb scattering

July 19, 2021

GoLuckyRyan Basic , , Leave a comment

This is a classical problem that is well know. This is more like my personal note. well, this blog is always my personal note. anyway,

The setting of the two-body problem is simple, there are two particles at position \vec{r}_1, \vec{r}_2 with velocity \vec{v}_1 = v \hat{x}, \vec{v}_2 = 0 , mass m_1, m_2 , and charge q_1, q_2 . There is a Coulomb force acting between the two particle and no external force.

The equations of motion for each particle are

\displaystyle m_1 \frac{d^2 \vec{r}_1}{dt^2} = \frac{q_1 q_2 e^2}{|\vec{r}_1 - \vec{r}_2|^{3/2}} \left( \vec{r}_1 - \vec{r}_2 \right)

\displaystyle m_2 \frac{d^2 \vec{r}_2}{dt^2} = \frac{q_1 q_2 e^2}{|\vec{r}_1 - \vec{r}_2|^{3/2}} \left( \vec{r}_2 - \vec{r}_1 \right)

Define the center of mass

\displaystyle \vec{R} = \frac{m_1 \vec{r}_1+ m_2 \vec{r}_2 }{m_1+m_2}, ~~~\displaystyle \frac{d^2\vec{R}}{dt^2} = 0

The particle-1 is moving at velocity v \hat{x} at impact parameter b . Shift to the center of mass frame,

\displaystyle \vec{\rho}_1 = \vec{r}_1 - \vec{R} = \frac{m_2}{m_1+m_2} (\vec{r}_1 - \vec{r}_2)

\displaystyle \vec{\rho}_2 = \vec{r}_2 - \vec{R} = \frac{m_1}{m_1+m_2} (\vec{r}_2 - \vec{r}_1)

The initial velocity at the center of mass frame is

\displaystyle \vec{u}_1 = \frac{m_2}{m_1+m_2} v \hat{x},~~ \vec{u}_2 = - \frac{m_1}{m_1+m_2} v \hat{x}

The total angular momentum with respect to the center of mass \vec{R} is

\displaystyle L = m_1 \vec{y}_1 \times \vec{u}_1 + m_2 \vec{y}_2 \times \vec{u}_2 = \frac{m_1 m_2^2}{(m_1+m_2)^2} vb + \frac{m_2 m_1^2}{(m_1+m_2)^2} vb = \mu v b

where \mu = m_1 m_2 / (m_1+m_2 ) is the reduced mass.

The equation of motion for the relative position is

\displaystyle m_1 m_2 \frac{d^2 (\vec{r}_1 - \vec{r}_2) }{dt^2} = (m_1+m_2) \frac{q_1 q_2 e^2}{ |\vec{r}_1 - \vec{r}_2|^{3/2}} \left( \vec{r}_1 - \vec{r}_2 \right)

\displaystyle \mu \frac{d^2 \vec{r} }{dt^2} = \frac{Z}{r^2} \hat{r}

using polar coordinate,

\displaystyle \vec{r} = r \hat{r} \\ \ddot{\vec{r}} =\dot{r} \hat{r} + r \dot{\theta} \hat{\theta} \\ \ddot{\vec{r}} = (\ddot{r} - r \dot{\theta}^2)\hat{r} + (2\dot{r}\dot{\theta} + r \ddot{\theta}) \hat{\theta}

We have

\displaystyle \ddot{r} - r \dot{\theta}^2 = \frac{Z}{\mu r^2}, ~~~ 2\dot{r}\dot{\theta} + r \ddot{\theta} = 0

From the angular equation,

2\dot{r}\dot{\theta} + r \ddot{\theta} = \frac{1}{r} \frac{r^2 \dot{\theta}}{dt}= 0

And we know that \mu r^2 \dot{\theta } = L , the angular momentum. Define, L = \mu l . Replace \dot{\theta} = l /r^2 , we have

\displaystyle \ddot{r} = \frac{l^2}{r^3} + \frac{Z}{\mu r^2}

The energy equation is

\displaystyle \frac{1}{2}\mu (\dot{r}^2 + (r\dot{\theta})^2 ) + \frac{Z}{r} = T_{cm} = \frac{1}{2}\mu v^2

replace the \dot{\theta} = l/r^2 ,

\displaystyle \dot{r}^2 + \frac{l^2}{r^2} + 2 \frac{Z}{\mu r} =2 \frac{T_{cm}}{\mu}

Differentiate the above equation will get back the equation of motion.

At the closest distance between the two particle r = b' , \dot{r} = 0 , l = vb

\displaystyle \frac{v^2 b^2}{b'^2} + 2 \frac{Z}{\mu b'} = v^2

Solve for b', we get back the Coulomb correction for Glauber model

\displaystyle b' = \eta' + \sqrt{\eta'^2 + b^2}, ~~ \eta' = \frac{Z}{2 T_{cm}}

where \eta' is the half distance of closest approach.

Back to the equation of motion. it may be not possible to solve it analytically in term of time. But it can be solve it in term of angle by replacing r = 1/u .

\displaystyle \frac{dr}{dt} = - \frac{1}{u^2} \frac{du}{d\theta} \frac{d\theta}{dt} = - r^2 \dot{\theta} u'(\theta) = - l u'

\displaystyle \frac{d^2 r}{dt} = - l u'' \dot{\theta} = - l^2 u^2 u''

The equation of motion becomes,

\displaystyle - l u^2 u'' = l u^3+ \frac{Z}{\mu} u^2 \rightarrow u'' + u = - \frac{Z}{\mu l^2}

The solution for this type of differential equation is

\displaystyle \frac{1}{r(\theta)} = A \cos\theta + B \sin\theta - \frac{Z}{\mu l^2}

or

\displaystyle \frac{1}{r} = \frac{Z}{\mu l^2 } \left( e \cos(\theta - \phi) - 1 \right) = \kappa \left( e \cos(\theta - \phi) - 1 \right)

in order to solve the equation of motion completely, the Initial condition has to be specified. Let our locus be like this:

The initial condition is \theta \rightarrow \pi, r \rightarrow \infty, r \sin(\theta) = b, dr/dt = v

The last condition can be rewrite as

\displaystyle \frac{dr}{dt} = - l \frac{du}{d\theta} = v \Rightarrow \frac{du}{d\theta} = -\frac{1}{b}

So, we solve the coupled equations:

\displaystyle ~~~~~0 = \kappa ( e \cos(\theta - \phi) -1 ) \\ - \frac{1}{b} = - \kappa e \sin(\theta-\phi)

\displaystyle \frac{1}{b'} = \frac{Z}{\mu l^2 } \left( e \cos(- \phi) - 1 \right)

The solution is

\displaystyle e = \frac{\sqrt{1 + b^2 \kappa^2}}{ b \kappa } = \sqrt{\frac{b^2}{\eta'^2}+1}

\displaystyle \phi = \frac{\pi}{2} + \tan^{-1}(\kappa b) = \tan^{-1}(-\kappa b, 1)

proton-proton Coulomb scattering in the relative position. The energies of the upper locus variate from 1, 2, 3, …, 9 MeV and the impact parameters are fixed at 1 fm. The energies of the lower locus are fixed at 5 MeV, and the impact parameters variates from 0.1, 0.2, … 0.9 fm.

From the above calculation. We can see there is a “shadow” behind the target. I generated many locus with different impact parameters with fixed energy of 5 MeV for proton-proton Coulomb scattering.

How to calculate the shadow???

With the equation of locus

\displaystyle \frac{1}{r} = \kappa \left( e \cos(\theta- \phi) -1 \right), \\~ \kappa = \frac{Z}{\mu l^2} = \frac{\eta'}{b^2}, ~~ e = \sqrt{\frac{b^2}{\eta'^2}+1}, ~ \phi = \frac{\pi}{2} + \tan^{-1}(\frac{\eta'}{b}) .

We can derived the Rutherford scattering.

The deflection angle is

\displaystyle \Theta = 2 \phi - \pi =2 \tan^{-1}\left( \frac{Z}{\mu v^2 b} \right)

express in term of b

\displaystyle \frac{Z}{\mu v^2 b} = \tan\left(\frac{\Theta}{2} \right)

\displaystyle b = \frac{Z}{\mu v^2} \cot\left( \frac{\Theta}{2} \right) = \eta' \cot\left( \frac{\Theta}{2} \right)

The last step used 2T_{cm} = \mu v^2 . The differential cross section is the function of a uniform ring with impact parameter.

\displaystyle 2\pi b db = - \frac{d\sigma}{d\Omega} d\Omega

with the cylindrical symmetry, d\Omega = 2\pi \sin \Theta d\Theta . The minus sign in the above equation reflect the fact that the change of impact parameter result smaller deflection angle.

\displaystyle \frac{d\sigma}{d\Omega} = - \frac{b}{\sin \Theta} \frac{db}{d\Theta} = \frac{b}{\sin\Theta} \frac{\eta'}{2 \sin^2(\Theta/2)}

Using \sin(x) = 2 \sin(x/2) \cos(x/2)

\displaystyle \frac{d\sigma}{d\Omega} = \frac{\eta'^2}{4 \sin^4 (\Theta/2)}

Using the trajectory equation, the nearest distance happens at \theta = \phi .

\displaystyle \frac{1}{b' } = \kappa ( e - 1) = \frac{\eta'}{b^2} \left( \sqrt{\frac{b^2}{\eta'^2}+1} -1 \right)

rearrange, we get

\displaystyle b' = \eta' + \sqrt{\eta'^2 +b^2}

The only remining thing are the trajectory in CM frame and Lab frame. Recall that, in the CM frame,

\displaystyle \vec{\rho}_1 = \frac{m_2}{m_1+m_2} (\vec{r}_1 - \vec{r}_2) , ~ \vec{\rho}_2 = \frac{m_1}{m_1+m_2} (\vec{r}_2 - \vec{r}_1).

Thus, simply multiple a factor, we can get the CM frame trajectory.

\displaystyle \rho_1 = \frac{m_2}{m_1+m_2} r

Also, the impact parameter with respect to the origin in the CM frame will be smaller,

\displaystyle b \rightarrow \frac{m_2}{m_1+m_2} b, ~~ \kappa = \frac{Z}{\mu l^2} \rightarrow \frac{Z}{\mu l^2} \left(\frac{m_1+m_2}{m_2} \right)^2 = \kappa \left(\frac{m_1+m_2}{m_2} \right)^2

The over all effect is \kappa \rightarrow \frac{m_1+m_2}{m_2} \kappa , the trajectory is

\displaystyle \frac{1}{\rho} = \left(\frac{m_1+m_2}{m_2} \right) \kappa \left( e \cos(\theta -\phi) -1 \right)

so the differential cross section is

\displaystyle \frac{d\sigma_{cm}}{d\Omega} = \left(\frac{m_2}{m_1+m_2}\right)^2 \frac{d\sigma}{d\Omega}

We can see, if the target mass is much bigger than the projectile mass, the correction is small. And we can see that the target will not recoil from the center, so \vec{r} = \vec{r}_1 - \vec{r}_2 \approx \vec{r}_1 . And for the proton-proton Coulomb scattering, in the CM frame, the cross section becomes quarter smaller. This is because the target will recoil and move, and that reduced the cross section, i.e. the field from the target is not rigid and becomes soft.

In order to see/verify the locus or the mass correction in CM frame. I solve the coupled equations using Mathematica,

\displaystyle m_1 \frac{dx_1^2}{dt^2} = \frac{Z}{ \left((x_1-x_2)^2 + (y_1-y_2)^2\right)^{3/2}} (x_1 - x_2)

\displaystyle m_1 \frac{dy_1^2}{dt^2} = \frac{Z}{ \left((x_1-x_2)^2 + (y_1-y_2)^2\right)^{3/2}} (y_1 - y_2)

\displaystyle m_2 \frac{dx_2^2}{dt^2} = \frac{Z}{ \left((x_1-x_2)^2 + (y_1-y_2)^2\right)^{3/2}} (x_2 - x_1)

\displaystyle m_2 \frac{dy_2^2}{dt^2} = \frac{Z}{ \left((x_1-x_2)^2 + (y_1-y_2)^2\right)^{3/2}} (y_2 - y_1)

These set of equations is same in all frame. Different frame has different initial conditions.

The CM frame initial condition are:

\displaystyle (x_1, y_1)(0) = \left(-\infty, \frac{m_2}{m_1+m_2} b \right), (x_1', y_1')(0) = \left(\frac{m_2}{m_1+m_2} v, 0 \right)

\displaystyle (x_2, y_2)(0) = \left(+\infty, -\frac{m_1}{m_1+m_2} b \right), (x_2', y_2')(0) = \left(-\frac{m_1}{m_1+m_2} v, 0 \right)

In Mathematica, it is, we cannot put infinity in the numerical solve, so, I set the initial position very far away.

In the below calculation, x_1(0) = -150, x_1(0) = 150, v = 0.0462127 . Since we use the nuclear unit, the unit of the speed is speed-of-light. The value corresponding to 1, 3, and 5 MeV proton energy. The impact parameter is 2 fm.

The solid locus are from the numerical solution. The dashed line is from the equation:

\displaystyle \frac{1}{\rho} = \frac{m_1+m_2}{m_2} \kappa ( e \cos(\theta-\phi) - 1)

I compare the target-fixed frame (upper), and CM frame (lower) trajectory in below:

For the Lab frame, the initial conditions are

\displaystyle (x_1, y_1)(0) = \left(-\infty, b \right), (x_1', y_1')(0) = \left( v, 0 \right)

\displaystyle (x_2, y_2)(0) = \left(0, 0 \right), (x_2', y_2')(0) = \left(0, 0 \right)

In the below simulation, the energy of proton is 1 MeV, the initial position of the projectile locates at -2000 fm.

Since the Coulomb force has infinite range, even at -2000 fm away, it still push the target proton away.

I guess it is almost most of the things about Two-body Coulomb scattering. We studied the locus in 3 frames: The target-fixed frame, the CM frame, and the Lab frame. There are few things to add into the calculation, e.g. lattice effect of the target that the target is bounded by a quadratic well, and screening of charge by electron. I will do the Rutherford differential cross section in Lab frame in other post. The idea is already in this post.

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.

Scattering with Coulomb potential

April 21, 2020

GoLuckyRyan Basic, Scattering , Leave a comment

I only stay the result in here. The result is very similar to the case without Coulomb potential.

The general solution for Coulomb potential is

\displaystyle \Phi_c( r \rightarrow \infty )= e^{i \vec{k} \cdot \vec{r}} e^{-i \eta \ln(k(r-z))} + f_c(\theta) \frac{e^{ikr}}{r} e^{-i \eta \ln(k(r-z))}

\displaystyle f_c(\theta) = - \frac{\eta}{ 2k \sin^2(\theta/2)} e^{-i \eta \ln(\sin^2(\theta/2))} e^{2i \sigma_0 } ,

\displaystyle \sigma_l = \arg( \Gamma(1+l+i\eta) )

In partial wave,

\displaystyle \Phi_c = A \frac{1}{2kr} \sum_{l=0}^\infty (2l+1) i^{l+1} P_l(\cos\theta) e^{i\sigma_l} \left( H_l^{-} - H_l^+\right)

Scattering with Coulomb and short-range potential, the potential is

\displaystyle V(r) = V_c(r) + V_N(r)

The solution takes the form

\Psi(r) = \Phi_c + \Psi_N

The asymptotic form of \Psi_N is

\displaystyle \Psi_N \rightarrow A f_N(\theta) \frac{e^{ikr}}{r} e^{-i\eta \ln(2kr)}

\displaystyle f_N(\theta) = \frac{1}{2ik} \sum_{l=0}^\infty (2l+1) P_l(\cos\theta) e^{2i\sigma_l} (S_l -1 )

S_l = e^{2i \delta_l}

\displaystyle \tan \delta_l = - \frac{kR F_l'(kR) - F_l(kR) L^I}{kR G_l'(kR) -G_l(kR) L^I}

The cross section is

\displaystyle \frac{d\sigma_c}{d\Omega} = |f_C(\theta) + f_N(\theta)|^2

Coulomb wave function (II)

April 17, 2020

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

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.

Coulomb Energy for two protons system

March 5, 2020

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

 

 

Coulomb energy for generic charge distribution

March 4, 2020

GoLuckyRyan Basic , Leave a comment

Classically, for an arbitrary charge distribution \rho(r,\theta, \phi) , the potential energy at position \vec{r}_0 is

\displaystyle V(r_0) = \int\int\int{\frac{\rho(r,\theta,\phi)}{|\vec{r}-\vec{r}_0|} r^2 \sin\theta d\phi d\theta dr}

Suppose the charge distribution is spherical symmetric, i.e. \rho(r,\theta,\phi) = \rho(r ). And the distance |\vec{r}-\vec{r}_0|^2 = r^2 + r_0^2 - 2 rr_0 \cos\theta .

\displaystyle V(r_0) = 2\pi\int{\rho(r)}r^2\int{\frac{1}{\sqrt{r^2 + r_0^2 - 2 rr_0 \cos\theta}}\sin\theta d\theta} dr

The angular part,

\displaystyle \int_{0}^{\pi}{\frac{1}{\sqrt{r^2 + r_0^2 - 2 rr_0 \cos\theta}}\sin\theta d\theta } \\ = \int_{1}^{-1}\frac{1}{\sqrt{r^2+r_0^2 - 2 r r_0 x}} dx \\ = \frac{-1}{2 rr_0} \int_{1}^{-1} \frac{d(r^2+r_0^2 - 2rr_0 x)}{\sqrt{r^2+r_0^2 - 2 r r_0 x}} \\ = \frac{1}{2 rr_0} \int_{(r-r_0)^2}^{(r+r_0)^2} y^{-\frac{1}{2}} dy = \frac{r + r_0 - |r - r_0|}{rr_0}

The last step is a bit tricky. So, the angular part is

\displaystyle \int_{0}^{\pi}{\frac{1}{\sqrt{r^2 + r_0^2 - 2 rr_0 \cos\theta}}\sin\theta d\theta } = \frac{r + r_0 - |r - r_0|}{rr_0} = \frac{2}{r_>}

Recall that

\displaystyle \frac{1}{|\vec{r}-\vec{r}_0|} = \sum_{L}^{\infty}\sum_{M=-L}^{L} \frac{4\pi}{2L+1} \frac{r_<^L}{r_>^{L+1}} Y_{LM}^*(\Omega)Y_{LM}(\Omega_0)

For L = 0 , that reduced to 4\pi/r_> .

Thus, for spherical arbitrary charge distribution, the Coulomb potential is,

\displaystyle V(r_0) = 4\pi\int{\rho(r)} \frac{r^2}{r_>} dr

Annotation 2020-04-07 081009.png

Lets redo the Coulomb energy for the uniform charge sphere here.

The charge distribution is

\displaystyle \rho_0(x) = Z \frac{3}{4\pi} \frac{1}{R^3} , x < R

The potential at position y is

\displaystyle V(y) = 4\pi \int_0^y \rho_0 (x) \frac{x^2}{y}dx + 4\pi \int_y^\infty \rho_0(x) x dx

For y < R

\displaystyle V(y) = Z \frac{3}{R^3} \left( \int_0^y \frac{x^2}{y}dx +  \int_y^R x dx \right) = Z \frac{3}{R^3} \left( \frac{y^2}{3} + \frac{R^2}{2} - \frac{y^2}{2} \right) \\ = \frac{Z}{R} \left( \frac{3}{2} - \frac{y^2}{2 R^2} \right)

For y > R

\displaystyle V(y) = 4\pi \int_0^R \rho_0 (x) \frac{x^2}{y}dx = \frac{Z}{y}

The result is the same.

The Classical Coulomb energy for another charge distribution \rho_0(\vec{r}_0) is

\displaystyle W = \int V(r_0) \rho_0(\vec{r}_0) d\vec{r}_0

Quantum Mechanically, the Coulomb energy between two wavefunction,

\displaystyle W = \left<\psi_1(r_1) \psi_2(r_2)|V_{12}| \psi_1(r_1) \psi_2(r_2)  \right> \\ = \int e \psi_1^2(r_1) \left(\int \psi_2^2(r_2) \frac{e}{|\vec{r}_1 - \vec{r}_2|} d\vec{r}_2 \right) d\vec{r}_1  \\ = \int \rho_1(r_1) V(r_1) d\vec{r}_1

We can see the expression is consistent with the Classical case.

For fermion N-body, the total wave function (in the 3rd part of this post) is,

\displaystyle \Psi = \frac{1}{\sqrt{N!}}\sum_{n}^{N!} (-1)^{p_n} P_n \phi_a(1)\phi_b(2)...\phi_z(N)  = (N!)^{-1/2}\sum_{n}^{N!} (-1)^{p_n} P_n \phi_H

where P_n is the permutation operator with parity p_n.

The Column energy is

\displaystyle W = \left<\Psi| V_{12} + V_{23} +.... V_{N-1,N} |\Psi \right> = \frac{N(N-1)}{2} \left<\Psi|V_{12}|\Psi \right>

\displaystyle \left<\Psi|V_{12}|\Psi \right> = (N!)^{-1} \sum_{i}^{N!} \sum_{j}^{N!} (-1)^{p_i+p_j+2} \left< P_i(\phi_H^*) |V_{12}|  P_j(\phi_H) \right>

Since the Coulomb operator V_{12} only act on particle 1 and 2, thus, the wave function for particle 3, 4,… N has to be the same. Thus, any none-zero product is exchange between particle 1 and 2. therefore, for any given permutation P_i the other permutation must be P_j = P_i - P_{12} P_i

\displaystyle W = (2(N-2)!)^{-1} \sum_{i}^{N!} \left< P_i(\phi_H^*) |V_{12}|  P_i(\phi_H) \right> - \left< P_i(\phi_H^*) |V_{12}|P_{12}P_i(\phi_H) \right>

In the N! permutation, for a fixed orbital of particle 1 and 2, there are (N-2)! way to permute the N-2 particle. Thus

\displaystyle W = \frac{1}{2} \sum_{a}^{N}\sum_{b\neq a}^{N-1} \left< ab|| ab \right> - \left<ab||ba\right>

or

\displaystyle W = \sum_{a}^{N}\sum_{b>a}^{N} \left< ab|| ab \right> - \left<ab||ba\right>

Here we use

\displaystyle \left<\phi_a^*(1)\phi_b^*(2)|V_{12}|\phi_h(1)\phi_k(2) \right> = \left<ab||hk\right>

for N=2

\displaystyle W = \sum_{a=1}^{2}\sum_{b=2}^{2} \left< ab|| ab \right> - \left<ab||ba\right> = \left< 12|| 12 \right> - \left<12||21\right>

which is the result we expect.

Coulomb Displacement Energy

February 12, 2020

GoLuckyRyan Basic , Leave a comment

I went back to my home town for visa renewal and read few books on history ( american revolution, early islam, modern chinese). After I got the visa a back to US, a lot of work has to be done and waiting for me….

The Coulomb displacement energy is the energy difference between isobaric analog states with same isospin. For example, the ground state energy different between 15O and 15N. The only difference between 15O and 15N is the proton in 0p1/2 shell replaced by a 0p1/2 neutron. Both ground state is in isospin T=1/2. The Coulomb displacement energy is defined as,

\displaystyle \Delta = BE(A, Z-1) - BE(A,Z) \\ ~~~~ = M(A,Z) - M(A,Z-1) + (m_n - m_p)

The binding energies are

BE(^{15}\textrm{O}) = 111.95535 ~\textrm{MeV}
BE(^{15}\textrm{N}) = 115.4919 ~\textrm{MeV}

The difference of the binding energy is 3.54 MeV.

Let’s calculate the Coulomb energy based on proton charge distributed uniformly in a sphere with radiu R=1.25 A^{1/3}. For A = 15, R = 3.083 \textrm{fm}. The Coulomb energy for uniform sphere is

\displaystyle E_c = \frac{3}{5}\frac{e^2}{R} Z(Z-1), e^2 = 1.44 ~\textrm{MeV} \cdot \textrm{fm}

Thus the Coulomb energies of 15O and 15N are

E_c(^{15}\textrm{O}) = 15.69 \textrm{MeV}
E_c(^{15}\textrm{N}) =   9.81 \textrm{MeV}

The difference is 4.11 MeV, about 0.5 MeV difference than the binding energy difference!

Should the electron mass add into the Coulomb displacement energy ??? or the electron mass be taken into account in further calculation? say, the Q-value for charge exchange reaction.

Simple model for 4He and NN-interaction

July 25, 2018

GoLuckyRyan Basic , , , , , , , Leave a comment

Starting from deuteron, the binding energy, or the p-n interaction is 2.2 MeV.

From triton, 3H, the total binding energy is 8.5 MeV, in which, there are only 3 interactions, two p-n and one n-n. Assume the p-n interaction does not change, the n-n interaction is 4.1 MeV.

The total binding energy of 3He is 7.7 MeV. The p-p interaction is 3.3 MeV.

Notices that we neglected the 3-body force in 3H and 3He. And it is strange that the n-n and p-p interaction is stronger then p-n interaction.

In 4He, the total binding energy becomes 28.3 MeV. I try to decompose the energy in term of 2-body, 3-body, and 4-body interaction.

If we only assume 2-body interaction, the interaction strength from n-p, n-n, and p-p are insufficient. One way to look is the 1-particle separation energy.

The neutron separation energy is 20.6 MeV = 2(p-n) + (n-n).
The proton separation energy is 19.8 MeV = 2(p-n) + (p-p).
The total energy is 28.3 MeV = 4(p-n) + (n-n) + (p-p).

There is no solution for above 3 equations. Thus, only consider 2-body interaction is not enough.

The neutron separation energy is 20.6 MeV = 2(p-n) + (n-n) + 2(n-n-p) + (n-p-p)
The proton separation energy is 19.8 MeV = 2(p-n) +(p-p) + 2(n-p-p) + (n-n-p).
The total energy is 28.3 MeV = 4(p-n) + (n-n) + (p-p) + 2(n-n-p) + 2(n-p-p) + (n-n-p-p).

Assuming the 2-body terms are the same in 2H, 3H, and 3He, the (n-p-p) and (n-n-p) is 4.03 MeV, which is strange again, as the Coulomb repulsion should make the (n-p-p) interaction smaller then the (n-n-p) interaction. The (n-n-p-p) interaction is -4.03 MeV.

Lets also add 3-body force in 3H and 3He.

The neutron separation energy of 3H is 6.3 MeV = (p-n) + (n-n) + (n-n-p)
The toal energy of 3H is 8.5 MeV = 2(p-n) + (n-n) + (n-n-p)
The toal energy of 3He is 7.7 MeV = 2(p-n) + (p-p) + (n-p-p)
The neutron separation energy of 4He is 20.6 MeV = 2(p-n) + (n-n) + 2(n-n-p) + (n-p-p)
The proton separation energy is 4He 19.8 MeV = 2(p-n) +(p-p) + 2(n-p-p) + (n-n-p).
The total energy of 4He is 28.3 MeV = 4(p-n) + (n-n) + (p-p) + 2(n-n-p) + 2(n-p-p) + (n-n-p-p).

We have 6 equations, with 6 unknown [ (p-n) , (n-n) , (p-p) , (n-n-p), (n-p-p), and (n-n-p-p)]. Notice that the equation from proton separation energy of 3He is automatically satisfied. The solution is

(p-n) = 2.2 MeV
(p-p) = – 4.7 – (n-n) MeV
(n-n-p)  = 4.1 – (n-n) MeV
(n-p-p) = 8 + (n-n) MeV
(n-n-p-p) = 0 MeV

It is interesting that there is redundant equation. But still, the (p-n) interaction is 2.2 MeV, and 4-body (n-n-p-p) becomes 0 MeV. Also the (n-p-p) is more bound than (n-n-p) by 3.9 + 2(n-n) MeV. If (n-p-p) should be more unbound, than (n-n) must be negative and smaller than -1.95 MeV.

Since the interaction strength has to be on the s-orbit (mainly), by considering 2H, 3H, 3He, and 4He exhausted all possible equations (I think). We need other way to anchor either (n-n), (p-p), (n-p-p), and (n-n-p) interactions.

Use the Coulomb interaction, the Coulomb interaction should add -1.44 MeV on the NN pair (assuming the separation is a 1 fm). Lets assume the (n-n) – (p-p) = 1.44 MeV

(n-n) = -1.63 MeV
(p-p) = – 3.07 MeV
(n-n-p)  = 5.73 MeV
(n-p-p) = 6.37  MeV

The (p-p) is more unbound than (n-n) as expected, but the (n-p-p) is more bound than (n-n-p) by 0.64 MeV. This is surprising! We can also see that, the 3-body interaction play an important role in nuclear interaction.

According to this analysis, the main contribution of the binding energies of 3H and 3He are the 3-body force.

In 3H:  (n-n) + 2(n-p) + (n-n-p) = -1.6 + 4.4 + 5.7 = 8.5 MeV
In 3He: (p-p) + 2(n-p) + (n-p-p) = -3.1 + 4.4 + 6.4 = 7.7 MeV

Worked on the algebra, when ever the difference  (n-n) – (p-p)  > 0.8 MeV, the (n-p-p) will be more bound that (n-n-p). Thus, the average protons separation should be more than 1.8 fm. I plot the interactions energies with the change of Coulomb energy below.

nn-pp.PNG

The (n-n) and (p-p) are isoscalar pair, where tensor force is zero. While the (n-p) quasi-deuteron is isovector pair. Thus, the difference between (n-n) and (n-p) reflect the tensor force in s-orbit, which is 3.8 MeV. In s-orbit, there is no spin-orbital interaction, therefore, we can regard the tensor force is 3.8 MeV for (n-p) isovector pair.

Following this method, may be, we can explore the NN interaction in more complex system, say the p-shell nuclei. need an automatic method. I wonder the above analysis agreed present interaction theory or not. If not, why? 

3D Harmonic oscillator

January 19, 2018

GoLuckyRyan Basic , , Leave a comment

Set x = r/\alpha The Schrodinger equation is

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

in Cartesian coordinate, it is,

\displaystyle  -\frac{\hbar^2}{2m}\left( \frac{d^2}{dx^2}+\frac{d^2}{dy^2}+\frac{d^2}{dz^2} \right) \Psi + \frac{1}{2} m \omega^2 (x^2+y^2+z^2) \Psi = E \Psi

We can set the wave function to be \Psi(r, \Omega)  = X(x) Y(y) Z(z)

\displaystyle  \left( -\frac{\hbar^2}{2m}\frac{d^2X}{dx^2} + \frac{1}{2} m \omega^2 x^2 X \right) YZ + \\ \left( -\frac{\hbar^2}{2m}\frac{d^2Y}{dy^2} + \frac{1}{2} m \omega^2 y^2 Y \right) XZ + \\ \left( -\frac{\hbar^2}{2m}\frac{d^2Z}{dz^2} + \frac{1}{2} m \omega^2 z^2 Z \right) XY = E XYZ

we can see, there are three repeated terms, we can set

\displaystyle -\frac{\hbar^2}{2m}\frac{d^2X}{dx^2} + \frac{1}{2} m \omega^2 x^2 X = E_x X

We decoupled the X, Y, Z. Each equation is a quadratic equation with energy

\displaystyle E_x, E_y, E_z = \left(\frac{1}{2} + n \right) \hbar \omega

and

\displaystyle E_x + E_y + E_z = \left(\frac{3}{2} +  n_x + n_y + n_z  \right) \hbar \omega = \left(\frac{3}{2} +  n  \right) \hbar \omega = E

The number of states for each energy level is

\displaystyle C^{n_x+n_y+n_z+2}_2 = C^{n+2}_2 = \frac{(n+2)!}{n!2!}

The first few numbers of states are 1, 3, 6, 10, 15, 21, 28, … The accumulated numbers of states are 1, 4, 10, 20, 35, 56, 84, … Due to the spin-state, the accumulated numbers of particles are 2, 8, 20, 40, 70, 112, 168, … The few magic numbers are reproduced.

The wave function is the product of the Hermite functions H_n(x) and exponential function

\Phi(x,y,z) = N H_{n_x} (x) H_{n_y}(y) H_{n_z}(z) \exp(-r^2/2)

If we simply replace (x,y,z) \rightarrow r( \cos(\phi) \sin(\theta), \sin(\phi) \sin(\theta) , cos(\theta) ) , we can see the ground state consists of s-orbit, the 1st excited state consists of p-orbit, and the 2nd excited state consists of d-orbit.

To see more clearly, we can project the function onto spherical harmonic, which is fixed angular momentum, i.e.

\langle \Phi(x,y,z) | Y_{lm}(\theta,\phi) \rangle = C_{lm} R_{nl}(r).

where C_{lm} is coefficients and R_{nl}(r) is radial function.

To have a better understanding, the radial function has to be solved. The procedure is very similar to solving Coulomb potential.

\displaystyle \left(-\frac{\hbar^2}{2m}\left(\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{d}{dr}\right)\right) + \frac{\hbar^2}{2m} \frac{L^2}{r^2} + \frac{1}{2} m \omega^2 r^2 \right)\Psi = E \Psi

separate the radial part and angular part.

\displaystyle \left( \frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{d}{dr}\right) - \frac{l(l+1)}{r^2} - \frac{m^2 \omega^2}{\hbar^2} r^2\right) R = - (2n+3) \frac{m \omega}{\hbar} R

Set \alpha^2 = \frac{\hbar}{m \omega}

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

Set x = r/\alpha

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

Set u(x) = x R(x)

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

as usual, the short range behaviour is r^{l+1}, long range behaviour is \exp(-x^2/2) , as stated in the Cartesian coordinate. Thus, we set

\displaystyle u(x) = f(x) \exp\left(-\frac{x^2}{2}\right) x^{l+1}

\displaystyle x\frac{d^2f}{dx^2} + (2(l+1) -2x^2) \frac{df}{dx} + 2x(n-l) f(x) = 0

with change of variable y = x^2 , the equation becomes

\displaystyle y\frac{d^2f}{dy^2} + \left( l + \frac{1}{2} + 1 - y \right)\frac{df}{dy} + \frac{n-l}{2} f(y) = 0

This is our friend, the Laguerre polynomial! In the Laguerre polynomial, (n-l)/2 \geq 0 must be non-negative integer. Now we set k = (n-l)/2 , than the energy is

\displaystyle E = \hbar \omega \left( 2k + l + \frac{3}{2} \right)

In order to have n, l, k are integer, when

n = 0 \rightarrow l = 0 \\ n = 1 \rightarrow l = 1 \\ n= 2 \rightarrow l = 0, 2 \\ n = 3 \rightarrow l = 1, 3 \\ n = 4 \rightarrow l = 0,2,4

The overall solution without a normalization factor is

\displaystyle \Psi_{nlm}(r, \theta, \phi) = r^l \exp\left(-\frac{r^2}{2\alpha^2}\right) L_{k}^{l+\frac{1}{2}}\left( \frac{r^2}{\alpha^2} \right) Y_{lm}(\theta, \phi)

The normalization constant can be calculate easily using the integral formula of Laguerre polynomial.

\displaystyle \int R^2(r) r^2 dr = 1

\displaystyle N^2 \int r^{2l} \exp\left(-\frac{r^2}{2\alpha^2}\right) \left( L_{k}^{l+\frac{1}{2}}\left( \frac{r^2}{\alpha^2} \right)\right)^2 r^2 dr  = 1

change of variable y = r^2/\alpha^2 \rightarrow r=\alpha \sqrt{y}

\displaystyle dr = \frac{\alpha}{2\sqrt{y}}dy

\displaystyle N^2 \frac{\alpha^{2l+3}}{2} \int y^{l+\frac{1}{2}} \exp(-y) \left( L_k^{l+\frac{1}{2}} \right)^2 dy = 1

The integration is

\displaystyle N^2 \frac{\alpha^{2l+3}}{2} \frac{(k+l+\frac{1}{2})!}{k!} = 1

\displaystyle N^2 = \frac{2}{\alpha^{2l+3}} \frac{k!}{(k+l+\frac{1}{2})!}

we can use

\displaystyle \left(n+\frac{1}{2}\right) = \frac{(2n+1)!}{2^{2n+1}n!}\sqrt{\pi}

Thus,

\displaystyle N^2 = \frac{1}{\sqrt{\pi}\alpha^{2l+3}} \frac{k! (k+l)! 2^{2k+2l+3}}{(2k+2l+1)!}

replace k = (n-l)/2 . The total wave function is

\displaystyle \Psi_{nlm}(r, \theta, \phi) \\ =\sqrt{ \frac{1}{\sqrt{\pi}\alpha^{2l+3}} \frac{(\frac{n-l}{2})! (\frac{n+l}{2})! 2^{n+l+2}}{(n+l+1)!}} r^l \exp\left(-\frac{r^2}{2\alpha^2}\right) L_{k}^{l+\frac{1}{2}}\left( \frac{r^2}{\alpha^2} \right) Y_{lm}(\theta, \phi)

Here are some drawing of the square of the wave functions. From the below is n = 0, 1, 2, 3, from left to right, are s-orbit, p-orbit, d-orbit, f-orbit.

3-D Harmonic Oscillator.png

With the LS coupling, the spatial function does not affected, unless the coupling has spatial dependence. With the LS coupling, the good quantum numbers are n, l, j, m_j

 

 

Older Entries