Deuteron Properties from AV18

Leave a comment

Haven’t posted anything in 2 months!! I have been focused to develop a few GUI DAQ programs using Qt6. And, I am buying a house! Anyway~


From the last post, we calculated the Deuteron wave function. The d-wave probability is \int w^2(r) dr and it is 5.8 % (experimental is 6.8%). The binding energy is -2.224574 MeV while the experimental value is -2.224575 MeV, only 1 eV difference!

The full wave function is

\displaystyle \Phi_M(r) = \frac{u(r)}{r} Y_{00} (\Omega) \chi_M + \frac{w(r)}{r} \sum_m C^{1M}_{2m1(M-m)} Y_{2m} (\Omega) \chi_{M-m}

where M= -1, 0, 1, The angular integrated square of wave function is

\displaystyle \int \Phi_M(r) \Phi_M^*(r) d\Omega = \frac{1}{r^2} ( u^2(r) + w^2(r) )

The matter radius or the root-mean-square radius is

\displaystyle \left<r^2 \right> = \frac{1}{2} \sqrt{\int_{0}^{\infty} (u^2(r) + w^2(r) ) r^2 dr } = 1.967~ \textrm{fm}

The factor of 2 is due to transforming from the relative motion to the center of mass motion and the masses of the proton and neutron are almost the same.


The dipole moment of deuteron has 2 contributions, one from the proton and neutron intrinsic spin, and the second part is from the orbital motion of the d-wave. The dipole moment is

\displaystyle \vec{\mu} = \mu_N \left(  \left(g_L \frac{\vec{L}}{2} + g_S \vec{S} \right) \cdot \vec{J} \right) \vec{J}

The orbital part needs to be half in CM frame than the relative frame, or g_l = 1/2.

The total spin \vec{S} is

\displaystyle \vec{\mu_S} = \mu_N \left(  \frac{g_p \vec{s_p} + g_n \vec{s_n}}{S(S+1)} \cdot \vec{S} \right) \vec{S} = \mu_N g_S \vec{S}

The S(S+1) is the magnitude of the vector \vec{S} . It is the same as projecting a classical vector \vec{u} to the direction of \vec{w} , we have to divide by |w|^2 as \vec{w} is not a uni-vector.

\displaystyle \vec{s_p} \cdot \vec{S} = \frac{1}{2} \left( S (S+1) + s_p(s_p+1) - s_n(s_n+1) \right) = 1, S = 1

Thus,

\displaystyle \vec{\mu_S} = \mu_N g_S \vec{S} = \mu_N \frac{g_p + g_n}{2} \vec{S}

The projection to \vec{J} is

\displaystyle \frac{\vec{\mu}}{\mu_N} = \frac{1}{J(J+1)} \left( \frac{1}{2} \vec{L} \cdot \vec{J} + \frac{g_p + g_n}{2} \vec{S} \cdot \vec{J} \right) \vec{J}

For L = 0, S = 1, J = 1, we have

\displaystyle \vec{\mu} = \mu_N \frac{g_p + g_n}{2} \vec{J}

Thus, the magnetic moment is

\displaystyle \mu = \left< \vec{\mu} \cdot \vec{J} \right> = \mu_N ( g_p + g_n)

For L = 2, S = 1, J = 1, we have

\displaystyle \vec{\mu} = \mu_N \frac{1}{2}\left( \frac{3}{4} - \frac{g_p + g_n}{2} \right)\vec{J}

Now, the magnetic moment is the sum of L=0 and L =2 state

\mu = \mu_N (g_p + g_n) \int u^(r) dr + \mu_N \left( \frac{3}{4} - \frac{g_p + g_n}{2} \right) \int w^2(r) dr = 0.84699 ~\mu_N

the experimental value is 0.857406~\mu_N.


The Quadrupole moment

\displaystyle Q = e \sqrt{\frac{\pi}{5}} \left< \Phi_1 |r^2 Y_{20}(\Omega)| \Phi_1 \right>

\displaystyle \left< \Phi_1 |r^2 Y_{20}(\Omega)| \Phi_1 \right> \\= \int \left( \frac{ u^2(r)}{4\pi} + \frac{u(r) w(r)}{\sqrt{10 \pi}} Y_{20} + \frac{w^2(r)}{10} (Y^2_{20} + 3|Y_{21}|^2 + 6 |Y_{22}|^2) \right) r^2 Y_{20}dr d\Omega

\displaystyle \left< \Phi_1 |r^2 Y_{20}(\Omega)| \Phi_1 \right> = \int \left(  \frac{u(r) w(r)}{\sqrt{10 \pi}} - \frac{w^2(r)}{20} \sqrt{\frac{5}{\pi}} \right) r^2 dr

\displaystyle Q = \frac{e}{20}  \int \left(  \sqrt{8} u(r) w(r) - w^2(r) \right) r^2 dr = 0.26966 ~ \textrm{e. fm}^2

The experimental value is 0.2859 e. fm2.

Deuteron Wave function using AV18

Leave a comment

It has been a long journey to this point. The basic facts for deuteron are discussed here, in which the tensor force plays an important role to bound the pn- pair and creating the d-wave. The coupled equation of the deuteron in the C.M. frame is derived here. The AV18 potential is used for the radial form of the central and tensor interaction. Now, the final step is solving it numerically.

The coupled equation is

\displaystyle  y''(r) = g(r) y(r),  \\ g(r) = \frac{1}{\alpha} \begin{pmatrix} V_1(r) - E & \sqrt{8} V_T(r) \\ \sqrt{8} V_T(r) & \alpha \frac{6}{r^2} + V_2(r) - E \end{pmatrix}

where E < 0 , \alpha = \frac{\hbar^2}{2m_\mu} = 41.47~\text{MeV} ,  m_\mu  = \frac{m_p m_n}{m_p + m_n}, and \int u^2 + w^2 dr = 1. The boundary conditions are u(0) = u(\infty) = w(0) = w(\infty) = 0 .

The coupled equation can be solved using Numerov’s method. However, the boundary conditions alone are not sufficient to solve the equation. We need 4 more, which is the derivative of u, w , the ratio of the two, and the eigen energy.

To solve it. we can use Numerov’s method with 2 sets of derivatives, and get 2 sets of solutions.

u(\delta r) = u_a , w(\delta r) = w_a  \rightarrow  u_1(r) , w_1(r)

u(\delta r) = u_b , w(\delta r) = w_b  \rightarrow  u_2(r) , w_2(r)

I used \delta r = 0.01 fm, u_a = w_b = 1, u_b = w_a = 0.1. The actual solution would be

u(r) = a u_1(r) + b u_2(r), \\w(r) = a w_1(r) + b w_2(r)

Next, at r \rightarrow \infty, the coupled equation becomes

\displaystyle  g(r \rightarrow \infty) = \frac{1}{\alpha} \begin{pmatrix} - E & 0 \\ 0 & \alpha \frac{6}{r^2} - E \end{pmatrix}

and decoupled. The asymptotic solutions are

\displaystyle f_s(r) = \exp( - \sqrt{|E|/\alpha} r), \\f_d(r) =   \exp( - \sqrt{|E|/\alpha} r) \left(1 + \frac{3}{r \sqrt{|E|/\alpha}} + \frac{3}{r^2 |E|/\alpha} \right)

Thus, we can relate

\displaystyle u(r) = a u_1(r) + b u_2(r) \rightarrow  \eta f_s(r), \\w(r) = a w_1(r) + b w_2(r) \rightarrow \zeta f_d(r)

and also the derivative

\displaystyle  u'(r) = a u'_1(r) + b u'_2(r) \rightarrow \eta f'_s(r), \\w'(r) = a w'_1(r) + b w'_2(r) \rightarrow \zeta f'_d(r)

Let r = x be the maximum radial distance for the numerical solution. We have

\displaystyle  a u_1(x) + b u_2(x) = \eta f_s(x), \\a w_1(x) + b w_2(x) = \zeta f_d(x), \\ a u'_1(x) + b u'_2(x) = \eta f'_s(x), \\ a w'_1(x) + b w'_2(x) = \zeta f'_d(x)

In matrix form

\displaystyle A = \begin{pmatrix} u_1(x) & u_2(x) & f_s(x) & 0 \\ w_1(x) & w_2(x) & 0 & f_d(x) \\ u'_1(x) & u'_2(x) & f'_s(x) & 0 \\  w'_1(x) & w'_2(x) & 0 & f'_d(x) \end{pmatrix},  \vec{v} = (a, b, -\eta, -\zeta)

such that

\displaystyle A \vec{v} = 0

For a non-zero solution for a null space equation, the determinant of the matrix A must be zero. So, by solving the u_i, w_i with different values of E, the eigen energy can be found when \det(A) = 0 .

The eigen energy is found to be -2.224574206 MeV. (experimental value is -2.224575 MeV. )

Once we have the eigen energy. We can solve the null space equation A\vec{v} = 0, in fact, we found that the 1st column vector and the 2nd column vector of the matrix $latex $ are proportional. Because of that, we can eliminate a in the 4 coupled equation of A\vec{v} = 0, and get

\displaystyle b =  \frac{f_s(x) u'_1(x) - f'_s(x) u_1(x)}{ u_2(x) u'_1(x) - u'_2(x) u_1(x)} \eta

Treat \eta as a free parameter. Eliminate \eta, we have

\displaystyle a = -b \frac{f_d(x) w'_2(x) - f'_d(x) w_2(x)}{ f'_d(x) w_1(x) - f_d(x) w'_1(x) }

or, the A matrix is

\displaystyle A(x) = \begin{pmatrix} 1.255 \times 10^8 & 9.693 \times 10^8 & 0.000964826 & 0 \\ -3.159 \times 10^9 & -2.4395 \times 10^{10} & 0 & 0.00144172 \\ 2.9065 \times 10^7 & 2.24486 \times 10^8 & -0.00022346 & 0 \\  -7.8296 \times 10^8& -6.04724\times 10^9 & 0 & -0.000351821 \end{pmatrix}

We can see that the ratio between column 1 and column 2 is 0.129474.

The null space is

\displaystyle \vec{v} = (0.000918457, -0.000118917, 0.99969, 0.0248831)

Thus, we have the deuteron wave function using AV18 potential in numerical form. Finally, normalized with a constraint

\displaystyle \int u^2(r) + w^2(r) dr = 1

We have the final wave function as

\displaystyle \psi_M(r) = \frac{u(r)}{r} Y_{00}(\Omega) \chi_{1M}+ \frac{w(r)}{r} \sum_{m}Y_{2m}(\Omega) \chi_{1,M-m}, M = -1, 0, 1

The radial functions are

The d-wave is 5.76%, and the experimental value is 6.8% from the numerical solution.


Place holder for Dipole moment, Quadrupole moment, and matter radius.

The dipole vector is

\displaystyle \vec{\mu_d} = g_p \vec{\mu_p} + g_n \vec{\mu_n} + \frac{\vec{l}}{2}

The quadrupole moment is

\displaystyle Q_d = \frac{e}{5} \int\left( \frac{2 u w}{\sqrt{8}} + \frac{w^2}{4} \right) r^2 dr

Equations of the deuteron wavefunction

4 Comments

The total Hamiltonian of a deuteron system is

\displaystyle H = \frac{P_p^2}{2m_p}+\frac{P_n^2}{2m_n} + V(\vec{r_p}-\vec{r_n}) + V_{LS}(\vec{r_p}-\vec{r_n}) L\cdot S + V_T(\vec{r_p}-\vec{r_n}) S_{12}

Change to the center of mass frame

\displaystyle H_c = \frac{P^2}{2m} + V(\vec{r}) + V_{LS}(\vec{r}) L\cdot S + V_T(\vec{r}) S_{12}, ~~ m = \frac{m_pm_n}{m_p+m_n}

The LS coupling is

\displaystyle L\cdot S = \frac{1}{2} \left(j(j+1)-l(l+1)-s(s+1) \right) = \begin{cases} 0 & l=0 \\ -3 & l=2 \end{cases}

We used j = 1, s = 1 in the calculation above.

For a two spin-1/2 system, the tensor operator S_{12} is already discussed here. We recall the result here,

\displaystyle S_{12} = \frac{3(\vec{\sigma_1}\cdot \vec{r})(\vec{\sigma_2}\cdot \vec{r}) - r^2 (\vec{\sigma_1}\cdot \vec{\sigma_2})}{r^2}

in matrix form:

\displaystyle S_{12} = \begin{pmatrix} 3 \cos^2(\theta) -1 & \frac{3}{\sqrt{2}} e^{i \phi} \sin(2\theta) & 3 e^{2i\phi} \sin^2(\theta) \\ \frac{3}{\sqrt{2}} e^{-i \phi} \sin(2\theta) & 2(1-3 \cos^2(\theta) ) & -\frac{3}{\sqrt{2}} e^{i \phi} \sin(2\theta) \\ 3 e^{-2i\phi} \sin^2(\theta) & -\frac{3}{\sqrt{2}} e^{-i \phi} \sin(2\theta) & 3 \cos^2(\theta) -1 \end{pmatrix}

with the spin states are:

\displaystyle \chi_{11}=\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}, ~~ \chi_{10}=\begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix}, ~~\chi_{1-1}=\begin{pmatrix} 0 \\ 0 \\ -1 \end{pmatrix}

We can rewrite it in terms of spherical harmonic,

\displaystyle S_{12} = \frac{4\sqrt{\pi}}{\sqrt{5}}\begin{pmatrix} Y_{2,0} & -\sqrt{3}Y_{2,1} & \sqrt{6} Y_{2,2} \\ \sqrt{3} Y_{2,-1} & -2 Y_{2,0} & \sqrt{3} Y_{2,1} \\ \sqrt{6} Y_{2,-2} & -\sqrt{3} Y_{2,-1} & Y_{2,0} \end{pmatrix}


Given that the deuteron total angular momentum J = 1, and the total spin is S=1, and the orbital angular momentum can be L=0,2. The wave function is spin-spherical harmonic,

\displaystyle \psi_{LS}^{JM} = \sum_{m_L, m_S} C_{Lm_L,Sm_S}^{JM} R_L(r) \chi_{Sm_S} Y_{Lm_L}(\theta, \phi)

We assumed the radial function does not depend on the magnetic sub-state (why?).

The S-wave is L=0

\displaystyle \psi_{01}^{1M} = R_0(r) Y_{00}(\theta,\phi) \chi_{1,M}

The D-wave is L = 2

\displaystyle \psi_{21}^{11} = \frac{R_2(r)}{\sqrt{10}} \begin{pmatrix} Y_{20} \\ -\sqrt{3} Y_{21} \\ \sqrt{6} Y_{22} \end{pmatrix}

\displaystyle \psi_{21}^{10} = \frac{R_2(r)}{\sqrt{10}} \begin{pmatrix} \sqrt{3}Y_{2-1} \\ -\sqrt{2} Y_{20} \\ \sqrt{3} Y_{21} \end{pmatrix}

\displaystyle \psi_{21}^{1-1} = \frac{R_2(r)}{\sqrt{10}} \begin{pmatrix} \sqrt{6}Y_{2-2} \\ -\sqrt{3} Y_{2-1} \\ Y_{20} \end{pmatrix}

The total wave function is

\displaystyle \Psi_{M} = a \psi_{01}^{1M} + b \psi_{21}^{1M}, ~~ a^2+b^2 = 1

Set a R_0(r) = U(r),  b R_2(r) = W(r)

\displaystyle \Psi_{1} = \begin{pmatrix} U(r) Y_{00} + \frac{1}{\sqrt{10}} W(r) Y_{20} \\ - \sqrt{\frac{3}{10}} W(r) Y_{21} \\ \sqrt{\frac{6}{10}} W(r) Y_{22}  \end{pmatrix}

\displaystyle \Psi_{0} = \begin{pmatrix} \sqrt{\frac{3}{10}} W(r) Y_{2-1} \\ U(r) Y_{00} - \sqrt{\frac{2}{10}} W(r) Y_{20} \\ \sqrt{\frac{3}{10}} W(r) Y_{21}  \end{pmatrix}

\displaystyle \Psi_{-1} = \begin{pmatrix} \sqrt{\frac{6}{10}} W(r) Y_{2-2} \\ - \sqrt{\frac{3}{10}} W(r) Y_{2-1} \\ U(r) Y_{00} + \sqrt{\frac{1}{10}} W(r) Y_{20}  \end{pmatrix}


Lets apply the Kinematic operator, the LS operator, and the tensor operator on the \Psi_{1},

\displaystyle \frac{P^2}{2m} \Psi_{1} = -\frac{\hbar^2}{2m} \left( \nabla^2 - \frac{L^2}{r^2} \right) \Psi_{1} = -\frac{\hbar^2}{2m} \left( \frac{1}{r^2} \left( \frac{d}{dr} r^2 \frac{d}{dr} \right) - \frac{L^2}{r^2} \right) \Psi_{1}

Before we proceed, we set

\displaystyle U(r) = \frac{u(r)}{r}, ~~~ \frac{1}{r^2} \left( \frac{d}{dr} r^2 \frac{d}{dr} \right) U(r) = \frac{1}{r} \frac{d^2}{dr^2} u(r) = \frac{u''(r)}{r}

\displaystyle \frac{P^2}{2m} \Psi_{1} = -\frac{\hbar^2}{2m} \frac{1}{r}\begin{pmatrix} u''(r) Y_{00} + \frac{1}{\sqrt{10}} \left( w''(r) Y_{20} - 6 \frac{w(r)}{r^2} Y_{20} \right) \\ - \sqrt{\frac{3}{10}} \left( w''(r) Y_{21} - 6 \frac{w(r)}{r^2} Y_{21} \right) \\ \sqrt{\frac{6}{10}} \left( w''(r) Y_{22} - 6 \frac{w(r)}{r^2} Y_{22} \right) \end{pmatrix}

\displaystyle L\cdot S \Psi_{1} =  - 3 \begin{pmatrix} \frac{1}{\sqrt{10}} W(r) Y_{20} \\ - \sqrt{\frac{3}{10}} W(r) Y_{21} \\ \sqrt{\frac{6}{10}} W(r) Y_{22}  \end{pmatrix}

\displaystyle S_{12}^* \Psi_{1} =\frac{4\sqrt{\pi}}{\sqrt{5}} \begin{pmatrix} U(r) Y_{00} Y_{20} + W(r) \frac{1}{\sqrt{10}} \left(Y_{20}^2-3Y_{2-1}Y_{21}+6Y_{22}Y_{2-2} \right) \\ -\sqrt{3} U(r) Y_{00} Y_{21} + W(r) \frac{1}{2\sqrt{5}} \left( \sqrt{6} Y_{21}Y_{20} + 6 Y_{2-1}Y_{22} \right) \\ \sqrt{6} U(r) Y_{00}Y_{22} + W(r) \frac{1}{\sqrt{10}} \left(-3Y_{21}^2+\sqrt{6}Y_{20}Y_{22} \right)  \end{pmatrix}

Using the equation of product of Spherical Harmonics,

\displaystyle Y_{l_1m_1}Y_{l_2m_2} = \sum_{lm} \sqrt{\frac{(2l_1+1)(2l_2+1)}{4\pi (2l+1)}} C_{l_10l_20}^{l0} C_{l_1m_1l_2 m_2}^{lm} Y_{lm}

The tensor term becomes

\displaystyle S_{12}^* \Psi_{1} =\frac{2}{\sqrt{5}} \begin{pmatrix} U(r) Y_{20} + W(r) \left(\sqrt{10}Y_{00}-\frac{1}{\sqrt{2}}Y_{20} \right) \\ -\sqrt{3} U(r) Y_{21} + W(r) \frac{\sqrt{3}}{\sqrt{2}} Y_{21} \\ \sqrt{6} U(r) Y_{22} - W(r) \sqrt{3} Y_{22}  \end{pmatrix}

In fact,

\displaystyle S_{12}^* \psi_{01}^{1M} = \sqrt{8} \psi_{21}^{1M}

\displaystyle S_{12}^* \psi_{21}^{1M} = \sqrt{8} \psi_{01}^{1M} - 2 \psi_{21}^{1M}


Now we can combine everything

\displaystyle H_c \Psi_{1}= \left(\frac{P^2}{2m}+ V(\vec{r}) + V_{LS}(\vec{r}) L\cdot S + V_T(\vec{r}) S_{12}^* \right) \Psi_{1}

Lets focus on the \chi_{11} component,

\displaystyle -\frac{\hbar^2}{2m} \frac{1}{r} \left( u''(r) Y_{00} + \frac{1}{\sqrt{10}} \left( w''(r) Y_{20} - 6 \frac{w(r)}{r^2} Y_{20} \right) \right) \\ + V(r) \left( U(r) Y_{00} + \frac{1}{\sqrt{10}} W(r) Y_{20} \right)  +3 V_{LS}(r) \frac{W(r)}{\sqrt{10}} Y_{20} \\ + V_T(r) \frac{\sqrt{4}}{\sqrt{5}} \left( U(r) Y_{20} + W(r) \left(\sqrt{10}Y_{00}- \frac{1}{\sqrt{2}}Y_{20}\right) \right) \\ ~~~~~~ = E \left( U(r) Y_{00} + \frac{1}{\sqrt{10}} W(r) Y_{20}  \right)

Change the U(r) = u(r)/r,  W(r) = w(r)/r

\displaystyle \frac{\hbar^2}{2m} \left( u''(r) Y_{00} + \frac{1}{\sqrt{10}} \left( w''(r) Y_{20} - 6 \frac{w(r)}{r^2} Y_{20} \right) \right) \\ - V(r) \left( u(r) Y_{00} + \frac{1}{\sqrt{10}} w(r) Y_{20} \right)  -3 V_{LS}(r) \frac{w(r)}{\sqrt{10}} Y_{20}  \\ -  V_T(r) \frac{\sqrt{4}}{\sqrt{5}} \left( u(r) Y_{20} + w(r) \left(\sqrt{10}Y_{00}- \frac{1}{\sqrt{2}}Y_{20}\right) \right) \\ ~~~~~~ = - E \left( u(r) Y_{00} + \frac{1}{\sqrt{10}} w(r) Y_{20}  \right)

Group the u(r), w(r) , set \hbar^2/2m = \alpha ,

\displaystyle   \left( \alpha  \frac{d^2}{dr^2} - V(r) + E  \right) u(r) Y_{00} - \sqrt{8} V_T(r) Y_{00} w(r) - V_T(r) \frac{2}{\sqrt{5}} u(r) Y_{20} \\ +  \left( \alpha \frac{d^2}{dr^2} - \alpha \frac{6}{r^2} - V(r) +3 V_{LS}(r) + E + 2 V_T(r) \right) \frac{w(r)}{\sqrt{10}} Y_{20}   = 0

The Y_{00} and Y_{20} are orthogonal. So, we have

\displaystyle  \left( \alpha  \frac{d^2}{dr^2} - V(r) + E  \right) u(r)= \sqrt{8} V_T(r)  w(r) — (1)

\displaystyle   \left( \alpha \frac{d^2}{dr^2} - \alpha \frac{6}{r^2} - V(r) +3 V_{LS}(r) + 2 V_T(r) +E \right) w(r) = \sqrt{8}V_T(r)  u(r) — (2)

with normalization

\displaystyle  \int_{0}^{\infty} u^2(r) + w^2(r) dr = 1

For the other components, I use Mathematica to calculate,

\Psi_M \chi_{1m_s} Commentequations
M = 1m_s = 1Group Y_{00}, Y_{20}(1), (2)
M = 1m_s = 0Only Y_{20}(2)
M = 1m_s = -1Take out the e^{-2i\phi}/\sin^2(\theta)(2)
M = 0m_s = 1Take out the e^{i\phi}/\sin(2\theta)(2)
M = 0m_s = 0Group Y_{00}, Y_{20}(1), (2)
M = 0m_s = -1Take out the e^{-i\phi}/\sin(2\theta)(2)
M = -1m_s = 1Take out the e^{2i\phi}/\sin^2(\theta)(2)
M = -1m_s = 0Take out the e^{i\phi}/\sin(2\theta)(2)
M = -1m_s = -1Group Y_{00}, Y_{20}(1),(2)

Spin-spin (Tensor) interaction

9 Comments

The spin-spin interaction, also called as tensor interaction, as the strength depends on the relative direction of the spins. The mathematical form of the interaction is

\displaystyle V_T(r) = \frac{1}{r^3} \frac{3(\vec{\sigma_1}\cdot\vec{r_1})(\vec{\sigma_2}\cdot\vec{r_2})-r^2(\vec{\sigma_1}\cdot\vec{\sigma_2})}{r^2}

For spin-1/2, we have.

\displaystyle \vec{\sigma} = \left(\sigma_x, \sigma_y, \sigma_z \right) = \left( \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}, \begin{pmatrix} 0 & i \\ -i & 0 \end{pmatrix}, \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \right)

The \vec{\sigma} is also called Pauli matrices.


The term \vec{\sigma}\cdot \vec{r} , we can expand,

\displaystyle \vec{\sigma}\cdot \vec{r} = \sigma_x x + \sigma_y y + \sigma_z z = \begin{pmatrix} z & x + i y \\ x-iy & -z \end{pmatrix}

If we express in spherical coordinate,

\displaystyle \vec{\sigma}\cdot \vec{r} =r \begin{pmatrix} \cos(\theta) &  \sin(\theta)e^{i\phi} \\ \sin(\theta)e^{-i\phi} & - \cos(\theta) \end{pmatrix}

The eigen values are 1 and -1, with eigen vectors are

\displaystyle v_{-1} = (\sin(\theta), -e^{-i\phi}\cos(\theta) )

\displaystyle v_{1} = (\cos(\theta), e^{-i\phi}\sin(\theta) )

with

\displaystyle v_i \cdot v_j^* = \delta_{ij}

The v_{1} is a “condensed” vector of the \vec{r'} = (\cos(\theta), \sin(\theta) \sin(\phi), \sin(\theta) \cos(\phi)) = (z, y, x) . Since only the the z-component is known, the x,y-components condensed as a phase \phi.


The term (\vec{\sigma_1}\cdot\vec{r_1})(\vec{\sigma_2}\cdot\vec{r_2}) is a tricky one. It involves 2 spins in two spaces. The product has 9 terms, \sigma_1(i) \sigma_2(j) r_i r_j . The \sigma_1 and \sigma_2 act on particle-1 and particle-2 separately. In fact, when we talk about a 2-spin system, we can form a direct-product space, \sigma_1(i) \sigma_2(j) = \sigma_1(i) \otimes \sigma_2(j) , in matrix from, for example

\sigma_1(x) \otimes \sigma_2(y) = \begin{pmatrix} 0 & \sigma_y \\ \sigma_y & 0 \end{pmatrix} = \begin{pmatrix} 0 & 0 & 0 & i \\ 0 & 0 & -i & 0 \\ 0  & i & 0 & 0 \\ -i & 0 & 0 & 0 \end{pmatrix}

The states in the direct-product space are

\displaystyle \left|\uparrow \uparrow \right> = (1, 0, 0, 0)

\displaystyle \left|\uparrow \downarrow \right> = (0, 1, 0, 0)

\displaystyle \left|\downarrow \uparrow \right> = (0, 0, 1, 0)

\displaystyle \left|\downarrow \downarrow \right> = (0, 0, 0, 1)

So, the eigen states for total spin 1 are (1, 0, 0, 0), \frac{1}{\sqrt{2}}(0, 1, 1, 0), (0,0,0,1)  , and for total spin 0 are \frac{1}{\sqrt{2}} (0, 1, -1, 0) . We set

\displaystyle \chi = \left( (1, 0, 0, 0), \frac{1}{\sqrt{2}}(0, 1, 1, 0), (0,0,0,1) , \frac{1}{\sqrt{2}} (0, 1, -1, 0) \right)

In similar fashion, the term \vec{\sigma_1}\cdot \vec{\sigma_2} = \sigma_1(x)\otimes\sigma_2(x) + \sigma_1(y)\otimes\sigma_2(y)+\sigma_1(z)\otimes\sigma_2(z)

\displaystyle \vec{\sigma_1}\cdot \vec{\sigma_2} = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & -1 & 2 & 0 \\ 0 & 2 & -1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}

Combine everything, we have

\displaystyle \frac{3(\vec{\sigma_1}\cdot\vec{r_1})(\vec{\sigma_2}\cdot\vec{r_2})-r^2(\vec{\sigma_1}\cdot\vec{\sigma_2})}{r^2} =\\ \frac{1}{r^2}\begin{pmatrix} 3z^2-r^2 & 3z(x+iy) & 3z(x+iy) & 3(r^2 - 2y^2-z^2+2ixy) \\ ... & r^2 - 3z^2 & r^2-3z^2 & -3z(x+iy) \\ ... & ... & r^2-3z^2 & -3z(x+iy) \\ ... & ... & ... & 3z^2 - r^2 \end{pmatrix}

Next, it is easy to calculate the matrix elements of the spin-spin interaction with the states \chi.

\displaystyle \left<\chi| V_T | \chi \right> = \frac{1}{r^2}\begin{pmatrix} 3z^2-r^2 & 3\sqrt{2} z(x+iy) & 3(x+iy)^2 & 0 \\ ... & 2(r^2 - 3z^2) & -3\sqrt{2}z (x+iy) & 0 \\ ... & ... & 3z^2-r^2 & 0 \\ ... & ... & ... & 0 \end{pmatrix}

We can see that, all matrix elements with spin-0 are 0. It is what we expect, as the spin-0 result no tensor force, as no spin in the system. For the diagonal elements:

\displaystyle \left< \chi_{Sm_S}| V_T | \chi_{Sm_S} \right> = \begin{cases} 2\left(1- 3 \frac{z^2}{r^2} \right), & S=1, m_S= 0 \\ 3 \frac{z^2}{r^2} -1, & S=1, m_S=\pm 1\\ 0, & S=0 \end{cases}


We can also calculate the matrix element of the \vec{\sigma_1}\cdot \vec{\sigma_2},

\left< \chi | \vec{\sigma_1}\cdot \vec{\sigma_2} |\chi \right> = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & -3 \end{pmatrix}

we can also calculate in other way. First, we convert to the spinor that we familiar,

s_i  = \frac{1}{2} \sigma_i

where the eigen value are 1/2, which is the correct eigen value for spin-1/2 system. And set

S = s_1 + s_2

so that the eigen values are 1 and 0.

\vec{\sigma_1}\cdot \vec{\sigma_2} = 4 s_1 \cdot s_2 =  2\left( S^2 - s_1^2 - s_2^2 \right)

Thus, we get

\displaystyle \left< \chi| \vec{\sigma_1}\cdot \vec{\sigma_2} | \chi\right> = \begin{cases} 1, & S=1\\ -3, & S=0 \end{cases} = 4 S - 3

This agrees with the matrix method.


The diagonal result can be written in spherical coordinate,

\displaystyle \left< \chi_{Sm_S}| V_T | \chi_{Sm_S} \right> = \begin{cases} 2\left(1- 3 \cos^2(\theta) \right), & S=1, m_S= 0 \\ 3 \cos^2(\theta) -1, & S=1, m_S=\pm 1\\ 0, & S=0 \end{cases}

We can see the angular dependent of the interaction. Let plot the 3 \cos^2(\theta) -1 ) for m_S = \pm 1 in below

At \theta = 0, \pi the force is strongest. This can picture as two magnets, aligned to the z-axis (m_S=\pm1). When the are aligned tail to head, or \theta = 0, \pi, the attraction is highest. At \theta = \pi/2 , when their head-to-head, tail-to-tail, they are repulsive.

Consider this distortion effect in deuteron. The proton and neutron couples to spin-1 and both are in the s-orbitals. Due to the tensor force, this s-orbital is no longer an eigenstate and it squeezes the orbital along the z-axis due to attraction and expands the x-y plan due to the repulsion. The total wave function would be oblate. Thus, the total wave function has to be mixed with the next orbital with the same parity, which is the d-orbital. And this is the case for m_s = \pm1.

For the case of m_s = 0 , the situation reversed, the tensor force expands the wavefunction along the z-axis and squeezes it on the x-y plane, making the total wavefunction be prolate. And since the magnitude is 2 times stronger for the m_s = 0 state than that of the m_s = \pm1 state. The deformation would be stronger.


[update 20230718]

In the case the total spin S= 1, m_S=0, the situation is like this:

This geometrical illustration explained why the tensor force change sign.

However, the x-y component of the total spin vector is unknown, so, the total spin can be rotated 90 degrees and look like this:

In this case, all tensor force is repulsive. Thus, the 2-spin system is not an axial deformed but also deformed on the x-y plane. However, after taking an average of all possible rotations, the strength on the x-y plane should be zero, and the deformation is axial.

We can see the tensor force, no matter how the spin is orientated, whenever tip to toe, it is attractive. whenever in parallel, repulsive. The m_S or the z-axis is artificially imposed so that we can do the calculation and make the tensor force operational. Or to say, the tensor force is rotationally symmetric.


The reason why the strength for the m_s = 0 case is twice that of the m_s = 1 is due to the internal degree of freedom on the x-y plane. When we sum up the strength, for m_s = 1, \theta = 0, \phi = 0, only 1 configuration is possible and has a strength of 2, but for m_s = 0, there are 4 configurations that have a strength of -1 each and sum up to -4. Similar for the \theta = \pi/2 as illustrated below.

But the problem is, the spin can rotate in any direction on the x-y plane, not just 4. Should it be 2\pi more instead of 2? Any better explanation?


In this post, the radial component of the tensor force AV18 p-n potential is negative, which explained why most nuclei tend to be prolated (longer in the z-axis) for axial deform.

Since different m_s states should degenerate without any external magnetic field. So that the different m_s states should be equally occupied. And the strength of the m_s = 0 state is twice that of the m_s = \pm 1 states, thus, the deformation effect along the z-axis of the tensor force should be canceled, but not on the x-y plane. As we discussed, for the m_s = 0 state, the x-y plane has no force on average due to rotation, but that is not true for the m_s = \pm1 states, the x-y plane is always attractive ( given the radial tensor force is negative), thus, the x-y plane is contracted, and make the whole nuclei prolate.

Why triton is unstable?

Leave a comment

There may be a simple answer….. if someone know, please leave a massage.


3H is unstable, it will undergo beta- decay, emit an electron and electron anti-neutrino and becomes 3He, which is stable. The Q value is very small, only 18.6 keV, assuming that the anti-neutrino has no mass. The half-life is 12.3 years.

Given that the deuteron is the only stable pn pair. 3H can be made by adding a neutron in a deuteron. And 3He can be made by adding a proton in a deuteron. Since the Coulomb force between protons, it is natural to think that 3He would be less bound than 3H. And it is,

BE(^{3}He) = 7.718043 MeV
BE(^{3}H) = 8.481798 MeV

The difference is 0.76376 MeV.

The masses are

M(^{3}He) = m_p + m_n + m_p - BE(^{3}He) = 2808.39 MeV
M(^{3}H) = m_p + m_n + m_n - BE(^{3}H) = 2808.92 MeV

The mass different between proton and neutron is 1.2933 MeV.

m_n - m_p - m_e = 0.782333 MeV

The binding energy usually connects to the stability. Higher binding energy should be more stable, and lesser mass for the same atomic number. It seems that the triton is a special case. Triton and 3He against this “common” sense.

Spectroscopic factor (again)

Leave a comment

In the last posts on SF, my mind have been biased to the mean field model. I started everything from the mean field. Now, lets forget the mean filed, just starting from experimental view point and any complete basis.


Suppose an arbitrary complete single-particle basis, say, a Woods-Saxon numerical basis, 3D spherical harmonic oscillator, or spherical Gaussian basis, and donate it as \phi_i . Using this basis,, we can form the Slater determinate of N-nucleons as

\displaystyle \Phi_i = \frac{1}{\sqrt{N!}}\begin{vmatrix} \phi_1(1) & ... & ... \\ ... & \phi_r(s) & ... \\ ... & ... & \phi_N(N) \end{vmatrix} .

The solution of any N-body Hamiltonian can be calculated by solving the eigen-system with the matrix elements

\displaystyle H_{ij} = \left< \Phi_i |H| \Phi_j \right>

Say, the solution is

\displaystyle  H \Psi_i = E_i \Psi_i


Experimentally, for the A(d,p)B reaction, we are observing the transition probability

\displaystyle T = \left< \chi_p \Psi_i(B)  |V|\chi_d \Psi_0(A) \right>,

where \chi_k is the distorted plane wave of  particle k and V is the interaction during the reaction. The integration summed all internal degree of freedom, particularly,

\displaystyle  \left<\Psi_i(B) |V_{BA} | \Psi_0(A) \right > \approx  \left<\Psi_i(B) | \Psi_0(A) \right > = \sum_i \alpha_i \phi_i = \psi

where \psi is called the quasi-particle wave function and it is expressed in the complete basis of \phi_i .

And then, the quasi-particle wave function will also be integrated.

\displaystyle T = \left< \chi_p \psi |V'|\chi_d \right>

In fact, it is the frame work of many calculations, that require 4 pieces of input,

  1. the optical potential for the incoming channel to calculate the distorted wave of the deuteron,
  2. the optical potential for the outgoing channel to calculate the distorted wave of the proton,
  3. the quasi-particle wave function, and
  4. the interaction.

If the Woods-Saxon basis is being used, then we can have shell model picture that, how the quasi-particle “populates” each Woods-Saxon states, and gives the spectroscopic factor.

The Woods-Saxon basis is somewhat “artificial”, it is an approximation. the nature does not care what the basis is, the transition probability is simple an overlap between two nuclei (given that the interaction is not strong, i.e. in the case of Direct reaction). Since the spectroscopic factor is always in reference for a basis, i.e., the quasi-particle wave function has to be projected into a basis, so, the spectroscopic factor is model dependence. If a weird basis is used, the shell model picture is lost. In this sense, SF is pure artificial based on shell model.

In experiment, we mainly observe the scattered proton, and from the angular distribution, the orbital angular momentum is determined. And from the conservation of angular momentum, we will know which orbital the neutron is being added to nucleus A.

But it is somehow weird. Who said the quasi-particle must be in an orbital angular momentum but not many momenta? If the spin of nucleus A is zero, because of conservation of angular momentum, the spin of nucleus B is equal to that of added neutron. But if the spin of nucleus A is not zero,  the neutron can be in many difference spins at once. In general, the conservation of angular momentum,

\displaystyle \left < \Psi_i(B) | \phi \Psi_0(A) \right >

could be difference than that of

\displaystyle \left < \chi_d \Psi_i(B) | V | \chi_p \Psi_0(A) \right >


We discussed that the spectroscopic factor is artificial. However, if we use a self-consistence basis, i.e. the basis that describes both nucleus A and B very well (is such basis exist?) , then the quasi-particle wave function is being described using the basis of nucleus B. To be explicit, lets say, we have a basis can describe nucleus B as

\Psi_0(A) = 0.8 \Phi_0 + 0.6 \Phi_1

\Psi_i(B) = \sum_{jk} \beta_{jk} \left( \phi_j \otimes  \Psi_k(A) \right)

where the single particle basis capture most of the nucleus A. is the spectroscopic factor meaningful in this sense?


So, is such basis exist? i.e, can it describe a nuclear wave function (almost) perfectly? In ab initial calculation, it can expend the nuclear wave function using some basis and give the overlap. Can the quasi-particle wave function be projected into the “single-particle state” ? How to define the “single-particle state” in ab initial calculation?


The experimental spectroscopic factor is calculated by compare the experimental cross section with that from theory. In the theory, the cross section is calculated by assuming the bound state is from a Woods-Saxon potential without correlation. Thus. the spectroscopic factor is with respect to the Woods-Saxon potential, an approximation.

Some thoughts on the quenching of spectroscopic factor

Leave a comment

Spectroscopic factor plays the central role in unfolding the nuclear structure. In the simplest manner, the total Hamiltonian of the nucleus is transformed into a 1-body effective potential and the many-body residual interaction, i.e.,

\displaystyle H = \sum_i^N \frac{P_i^2}{2m_i} + \sum_{i \neq j}^N V_{ij} = \sum_i^N  \left( \frac{P_i^2}{2m_i} + U \right) + \sum_{i\neq j}^N \left( V_{ij} - U \right) \\ = \sum_{i}^N h_i + H_R = H_0 + H_R

The effective single-particle Hamiltonian has solution:

\displaystyle h_i \phi_{nlj}(i) = \epsilon_{nlj} \phi_{nlj}(i)

where \epsilon_{nlj} is the single-particle energy. The solution for H_0 is

\displaystyle H_0 \Phi_k(N) = W_k \Phi_k(N)

\displaystyle \Phi_k(N)= \frac{1}{\sqrt{N!}}\begin{vmatrix} \phi_{p_1(k)}(1) & \phi_{p_1(k)}(2) & ... & \phi_{p_1(k)}(n) \\ \phi_{p_2(k)}(1) & \phi_{p_2(k)}(2) & ... & ... \\ ... & ... & ... & ... \\ \phi_{p_n(k)}(1) & \phi_{p_n(k)}(2) & ... & \phi_{p_N(k)}(N) \end{vmatrix}

\displaystyle W_k = \sum_i^N \epsilon_{p_i(k)}

where p_i(k) is the set of basis for state k from \phi_{nlj} , and W_k is the eigenenergy.

The residual interaction is minimized by adjusted the mean-field U. Thus, the residual interaction can be treated as a perturbation. This perturbs the nuclear wave function

H \Psi_k(N) = W_k \Psi_k(N),    \Psi_k(N) = \sum_i \theta_{i}(k) \Phi_i(N) .

The normalization requires \sum_i \theta_{i}^2(k) = 1 .


In the Slater determinant \Phi_k, a single-particle wave function for a particular orbital can be pull out.

\displaystyle \Phi_k(N) = \phi_{\mu} \otimes \Phi_{k}(N-1)

where \otimes is anti-symmetric, angular coupling operator. Thus,

\displaystyle \Psi_k(N) = \sum_{\mu i} \theta_{\mu i}(k) \phi_{\mu} \otimes \Phi_i(N-1)

The \theta_{\mu i}^2 (k) is the spectroscopic factor. There are another sum-rule for adding and removing a nucleon. so that the sum is equal to the number of particle in a particle orbital.


I always imagine the quenching is because we did not sum-up the SFs from zero energy to infinity energy (really???), thus, we are always only observing a small fraction of the total wave function. For example,  the total wavefunction would look like this:

\displaystyle \Psi_k(N) =  \phi_{0} \otimes \left(\theta_{00}(k) \Phi_0(N-1) + \theta_{01}(k)\Phi_1(N-1) +....  \right) \\ + \phi_{1} \otimes \left( \theta_{10} \Phi_0(N-1) +...   \right) +...   .

In experiment, we observe the overlap between ground-state to ground-state transition

\displaystyle  \left<  \phi_0 \Phi_0(N-1) | \Phi_0(N-1) \right> = \theta_{00}(0)

for ground-state to 1st excited state transition for the same orbital is

\displaystyle  \left<  \phi_0 \Phi_1(N-1) | \Phi_0(N-1) \right> = \theta_{01}(0)

And since we can only observed limited number of excited states, bounded by either or boht :

  • experimental conditions, say incident energy
  • the excited states that are beyond single-particle threshold.
  • finite sensitivity of momentum

Thus, we cannot recover the full spectroscopic factor. This is what I believe for the moment.


Experimentally, the spectroscopic factor is quenched by 40% to 50%. The “theory” is that, the short-range interaction quench ~25%, the long-range interaction quench ~20%. The long- and short-range interaction correlate the single-particle orbital and reduce the degree of “single-particle”.

Annotation 2020-03-31 010604.png

The short-range interaction is mainly from the “hard-core” of nucleon, i.e., the interaction at range smaller than 1 fm. The long-range interaction is coupling with nearby vibration states of the rest of the nucleus.

For example, from the 19F(d,3He) reaction, the spectroscopic factor for 19F 1s1/2 state is ~0.4, and 0d5/2 is ~0.6.

Annotation 2020-03-31 011328.png

Thus, the wavefunction of 19F is

\left|^{19}\textrm{F}\right> \approx  \sqrt{0.4} \left|1s_{1/2}\right> \otimes \left|^{18}\textrm{O}_{g.s.} \right> + \sqrt{0.6} \left|0d_{5/2} \right> \otimes \left|^{18}\textrm{O}(1.98) \right>

It is worth to note that the above SFs is not re-analysised and the “quenching” is not shown. Many old data had been re-analysised using global optical model and the SF is reduced and show that the sum of SFs is ~ 0.55.

If it is the case for 19F, the wavefunction would become,

\left|^{19}\textrm{F}\right> \approx  \sqrt{0.2} \left|1s_{1/2}\right> \otimes \left|^{18}\textrm{O}_{g.s.} \right> + \sqrt{0.3} \left|0d_{5/2} \right> \otimes \left|^{18}\textrm{O}(1.98) \right> + \sqrt{0.5} \Psi_k

Here I use \Psi_k for the “correlated wavefunction” that the single-particle orbital cannot simply pull out. Nevertheless, if x and y are correlated,

f(x,y) \neq  g(x) h(y)

Am I misunderstood correlation?


My problem is, What does a correlated wave function look like?

In my naive understanding, the Slater determinant \Phi_k is a complete basis for N-nucleon system. A particular single-particle orbital can ALWAYS be pull out from it. If it can not, therefore, the Slater determinant is NOT complete. The consequence is that all theoretical calculation is intrinsically missed the entire CORRELATED SPACE, an opposite of Slater determinant space (of course, due to truncation of vector space, it already missed somethings).

If the theory for correlation is correct, the short-range interaction is always there. Thus, the spectroscopic factor for deuteron 0s1/2 orbital is ~0.8, assuming no long-range correlation. However, we already knew that 96% of deuteron wavefunction is from 0s1/2 and  only 4% is from 1d5/2 due to tensor force. Is it not mean the spectroscopic factor of deuteron 0s1/2 state is 0.96?  Is deuteron is a special case that no media-modification of nuclear force? But, if the short-range correlation is due to the hard core of the nucleon, the media-modification is irrelevant. Sadly, there is no good data such as d(e,e’p) experiment. Another example is 4He(d,p)5He experiment. What is the spectroscopic factor for g.s. to g.s. transition, i.e. the 0p3/2 orbital? is it ~0.6 or ~ 1.0?

Since the experimental spectroscopic factor has model dependency (i.e. the optical potential). Could the quenching is due to incomplete treatment of the short- and long-range correlation during the interaction, that the theoretical cross section is always bigger?

In the very early days, people calibrate their optical potential using elastic scattering for both incoming and out-going channel, and using this to produce the inelastic one. At that time, the spectroscopic factors are close to ~1. But since each optical potential is specialized for each experiment. It is almost impossible to compare the SF from different experiments. Thus, people switch to a global optical potential. Is something wrong with the global optical potential? How is the deviation?


Let me summarize in here.

  1. The unperturbed wave function should be complete, i.e. all function can be expressed as a linear combination of them.
  2. A particular single-particle orbital can be pull out from the Slater determinate \Phi_k.
  3. The residual interaction perturbs the wave function. The short-/long-range correlation should be in the residual interaction by definition or by construction of the mean field.
  4. The normalization of wave function required the sum of all SF to be 1.
  5. Another sum rule of SFs equals to the number of particle.
  6. By mean of the correlation, is that many excited states have to be included due to the residual interaction. No CORRELATED space, as the Slater determinant is complete. (pt. 1)
  7. Above points (1) to (7) are solid mathematical statements, which are very hard to deny.
  8. The logical result for the quenching of the observed SF is mainly due to not possible to sum up all SFs from all energy states for all momentum space.
  9. The 2-body residual interaction can create virtual states. Are they the so called collective states?
  10. But still, collective states must be able to express as the Slater determinant (pt. 1), in which a particular single-particle orbital can be pull out (pt. 2).
  11. May be, even the particular single-particle orbital can be pull out, the rest cannot experimentally observed ? i.e. \Phi_k(N-1) is not experimentally reachable. That go back to previous argument for limitation of experiments (pt. 8).
  12. For some simple systems, say doubly magic +1, deuteron, halo-nucleon, very weakly bounded exited state, resonance state, the sum of SF could be close to 1. Isn’t it?
  13. The theoretical cross section calculation that, the bound state wave function is obtained by pure single-particle orbital. I think it is a right thing to do.
  14. The use of global optical potential may be, could be not a good thing to do. It may be the METHOD to deduce the OP has to be consistence, instead of the OP itself has to be universal. Need more reading from the past.

some thoughts on mass

Leave a comment

Lets forget about the Higg mechanism.

Lets start from deuteron, it contains a proton and neutron, which bounded by strong nuclear force. The kinematics energy between the proton and neutron is about 35 MeV. Since the binding energy of deuteron is 2.2 MeV, thus the strong nuclear force is about -37 MeV. According to Einstein E = mc^2 . This 2.2 MeV is converted to 2.2 MeV/c2 of mass.

In this example, when a system contains internal motion, internal interaction, the “mass” is not just the mass of individual components, but also included the interaction energy.

Thus, what is mass? Only elementary particles had “pure” mass. The mass of all other particles must contains the energy of internal interaction and motion.

So, for the Earth, we know the mass from the orbital motion. But this mass, when the Earth exploded into tiny many pieces, some of the mass was gone because the internal interaction and motion no longer there.

What is the mass of the solar system? We can add all the masses of the sun, the planets, comets, ices, debris, etc, but we also need to add the internal motion and interaction. Although the gravitation force is weak, but still, we need to add this.

What is the mass of the Universe? Again, we can added all visible mass, but we also need to add the dark energy, the dark matter ( I suspect the dark matter is kind of illusion that we don’t fully understood the gravity, can the dark matter is due to gravitation energy, which act as mass?).

Simple model for 4He and NN-interaction

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? 

Stopping power and Bethe-Bloch Formula

Leave a comment

I am so surprised that this topic is not in this blog.  But I am not always posting what I did. Anyway…


The Stopping power is energy loss per length, in unit of MeV/cm.

\displaystyle S = -\frac{dE}{dx}

Since the stopping power is proportional to density, it is often normalized with the density in literature and the unit becomes MeV/(ug/cm^2).

The Bethe-Bloch formula based on classical argument.


In order to calculate the range, we can integrate the stopping power. Consider the particle moved by \Delta x distance, the energy loss is

E(x+\Delta x) = E(x) - S(E(x)) \Delta x

rearrange, gives

\displaystyle  \frac{dE(x)}{dx} = - S(E)

\displaystyle \int_{E_0}^{E} \frac{-1}{S(E)} dE = x(E ; E_0)

This is the relation between particle energy and range.

We can plot S(E(x)) to get the Bragg peak.


In here, I use the physical calculator in LISE++, SRIM, and Bethe-Bloch formula to calculate the stopping power for proton in CD2 target.

The atomic mass of deuteron is 2.014102 u. The density of CD2 target is 0.913 g/cm3.

SvsE.PNG

For the Bethe curve, I adjusted the density to be 0.9 mg/cm3, and simple use the carbon charge. And the excitation potential to be 10^{-5} z, where z = 6 , so that the Bethe curve agree with the others.

It is well known that the Bethe curve fails at low energy.

Assume the proton is 5 MeV. The energy loss against distance isEvsX.PNG

and the stopping power against distance is

SvsX.PNG

We can see the Bragg peak is very sharp and different models gives different stopping range. It is due to the rapid decrease of energy at small energy. In reality, at small energy, the microscopic effect becomes very important and statistical. thus, the curve will be smoothed and the Bragg peak will be broaden.

In very short range, the stopping power increase almost linearly, Because the incident energy is large, so that the stopping power is almost constant around that energy range.

For S(E_0) = h ,

\displaystyle  x = \frac{-1}{h}(E - E_0 )

\displaystyle E = E_0 - hx,  hx << E_0

\displaystyle S(E) = S(E_0 - hx) \approx S(E_0) - \frac{dS(E_0)}{dE}(hx)

 

 

Older Entries