Experimental TBMEs

Leave a comment

For a doubly magic + 2 nucleons system, the 2 nucleons are in j_1 and j_2 orbitals. If there is no residual interaction, all possible J-states formed by j_1-orbital and j_2-orbital are degenerated. As the total Hamiltonian is a simple sum of “single-particle” Hamiltonian with the mean field, that nucleons “do” not interact with each others.

Lets denote the single-particle energies for the j_1 an j_2 orbitals are \epsilon_1 and \epsilon_2, respectively. The interaction energy between the two valence nucleons is:

\displaystyle E_{JM}(j_1j_2, j_1j_2) = \left<\Psi_{JM}(j_1J_2)|V_{12}| \Psi_{JM}(j_1J_2) \right>

where H is the nuclear Hamiltonian. And \Psi_{JM}(j_1J_2) is the wave function of two nucleons that coupled to spin J .

Sometimes, there could be 2 configurations couples to same J . for example, (d5/2)(d5/2) and (d5/2)(s1/2) can both couple to spin 2, and 3. For example, see this post. In this case, the wave function is a linear combination of the two configurations. Thus, the Hamiltonian is a bit complicated,

\displaystyle H = \begin{pmatrix} E_1 && V \\ V && E2 \end{pmatrix}

\displaystyle E_1 = \epsilon_1 + \epsilon_2 + E_{JM}(j_1j_2, j_1j_2)
\displaystyle E_2 = \epsilon_3 + \epsilon_4 + E_{JM}(j_3j_4, j_3j_4)
\displaystyle V =E_{JM}(j_1j_2, j_3j_4)

The solution is in this post. And the observed energy levels are e_\pm,

\displaystyle e_\pm = \bar{E} \pm \sqrt{ \Delta E^2 + V^2}

and the states have spectroscopic factors,

\displaystyle \Psi_+ = \alpha \Psi(j_1j_2) + \beta \Psi(j_3j_4)
\displaystyle \Psi_- = -\beta \Psi(j_1j_2) + \alpha \Psi(j_3j_4)

Now, the theoretical background is laid down. Experimentally, lets take the 18O , and the 17O(d,p)18O reaction as an example.

The 17O is a single d5/2 neutron on top of 16O. Adding another neutron on 17O, the new neutron will couple to that d5/2 neutron, any state contains d5/2 will be populated. We restrict ourselves only to the d5/2 and s1/2 states. For the J = 0 state, we pick the ground state and the 5.34 state as the mixing between the (d5/2)(d5/2) and (s1/2)(s1/2).

\displaystyle \Psi_0(j_1j_1) = [\phi_{5/2} \times \phi_{5/2}]_{J=0}
\displaystyle \Psi_0(j_2j_2) = [\phi_{1/2} \times \phi_{1/2}]_{J=0}

The spectroscopic factor of the ground state is 1.22, and the SF of the 5.34 MeV state is 0.16. We have to normalize the SF.

\displaystyle \Psi_(0.00) = \sqrt{\frac{1.22}{1.22+0.16}} \Psi_0(j_1j_1) + \sqrt{\frac{0.16}{1.22+0.16}} \Psi_0(j_2j_2) \\= 0.94 \Psi_0(j_1j_1) + 0.34 \Psi_0(j_2j_2)

\displaystyle \Psi_(5.34)= -0.34\Psi_0(j_1j_1) + 0.94 \Psi_0(j_2j_2)

Next, we have to estimate the single-particle energy. The binding energy for the a d5/2 neutron is S_n(17O) = epsilon_1= 4.14 MeV, And the binding energy for 2 neutrons and their interaction is S_{2n}(17O) = 12.19 MeV. Thus, the neutron interaction energy is S_n(17O) - S_{2n}(18O) = -3.90 MeV .

Thus, we have

e_1 = 0.00 - 12.19 = -12.19 MeV
e_2 = 5.34 - 12.19 = -6.85 MeV
\alpha = 0.94, \beta = 0.34


E_1 = -11.57 MeV = \epsilon_1 + \epsilon_1 + E_0(j_1j_1, j_1j_1)
E_2 = -7.47 MeV = \epsilon_2 + \epsilon_2 + E_0(j_2j_2, j_2j_2)
V = -1.71 MeV = E_0(j_1j_1, j_2j_2)

Thus, the TBME of the d5/2-d5/2 neutrons coupled to J = 0 is -3.29 MeV. The TBME of the d5/2-s1/2 neutrons coupled to J = 0 is -1.71 MeV.

The epsilon_2 is the single-particle energy of the s1/2 neutron. The 1/2+ state next to the 5/2+ ground state of 17O is 0.87 MeV. Thus, \epsilon_2 = 3.27 MeV, then, E_0(j_2j_2, j_2j_2) = -0.92 MeV.

For the 4+ state of 18O, it can only be coupled by (d5/2)(d5/2) neutron. Teh excited energy is 3.55 MeV, The TBME of J = 4 is -0.35 MeV.

For the 2+ state, we can repeat the same method. We take the 1.98 MeV and 3.92 MeV, with spectroscopic factors for the d5/2 neutron adding are 0.83 and 0.66 respectively. The 2+ state can be formed by (d5/2)(d5/2) and (d5/2)(s1/2) configuration.

e_1 = 1.98 - 12.19 = -10.21 MeV
e_2 = 5.26 - 12.19 = -8.27 MeV
\alpha = 0.75, \beta = 0.67


E_1 = -10.21 MeV = \epsilon_1 + \epsilon_1 + E_2(j_1j_1, j_1j_1)
E_2 = -8.27 MeV = \epsilon_1 + \epsilon_2 + E_2(j_1j_2, j_1j_2)
V = -0.96 MeV = E_2(j_1j_1, j_1j_2)

Then, the J=2 (d5/2)(d5/2) TBME is -1.06 MeV, the J = 2 TBME for (d5/2)(s1/2) is -1.71 MeV, and the off-diagonal TBME is -0.96 MeV.

At last, the 3+ state can only be formed by (d5/2)(s1/2). The excited state in 18O is 5.375 MeV. 5.375 = \epsilon_1 + \epsilon_2 + E_3(j_1j_2, j_1j_2) , thus, the J=0 TBME for (d5/2)(s1/2) is 0.60 MeV.

In summary,

JTconfigurationTBMEs [MeV]
01diagonald5/2 – d5/2-3.29
21diagonald5/2 – d5/2-1.06
41diagonald5/2 – d5/2-0.35
01diagonals1/2 – s1/2-0.92
01(d5/2)(d5/2) – (s1/2)(s1/2)-1.71
21diagonal(d5/2)(s1/2) – (d5/2)(s1/2)-1.71
21(d5/2)(d5/2) – (d5/2)(s1/2)-0.96
31diagonald5/2 – s1/20.60

The total Hamiltonian is

\displaystyle H = \begin{pmatrix} H_{J=0} (2x2) && 0 && 0 && 0 \\ 0 && H_{J=2} (2x2) && 0 && 0 \\ 0 && 0 && H_{J=3} && 0 \\ 0 && 0 && 0 && H_{J=4} \end{pmatrix}

\displaystyle H_{J=0} = \begin{pmatrix} -3.29 && -1.71  \\ -1.71 && -0.92 \end{pmatrix}

\displaystyle H_{J=2} = \begin{pmatrix} -1.06 && -0.96 \\ -0.96 && -1.71 \end{pmatrix}

\displaystyle H_{J=3} = 0.60

\displaystyle H_{J=4} = -0.35

We can compare these value with our previous estimation. For the (d5/2) (d5/2) configuration.

I state the previous estimation in here, with the single-particle energy already subtracted.

\displaystyle H'_{J=0} = \begin{pmatrix} -3.171 && -1.831  \\ -1.831 && -1.057 \end{pmatrix}

\displaystyle H'_{J=2} = \begin{pmatrix} -0.725 && -0.959 \\ -0.959 && -1.27 \end{pmatrix}

\displaystyle H'_{J=3} =0

\displaystyle H'_{J=4} = -0.302

The 4+ state interaction action energy is -0.35 MeV. In the previous estimation, it is -2/7 V_0 \bar(R) = -0.302 . which is very good agreement.

The 3+ state is actually repulsive. In the previous estimation, we ignored the mutual interaction, because the (d5/2) and (s1/2) states are eigen state and there is no mutual interaction. Possible issue mixing with d3/2 or the interaction is not a delta function.

The J=0 TBME for the (d5/2)(d5/2) state is -3.29 MeV and in previous estimation, it is -3V_0 \bar{R} = -3.17 MeV, good agreement. The J=0 TBME for the (d5/2)(d5/2) – (s1/2)(s1/2) is -1.71 MeV, and the J=0 TBME for (s1/2)(s1/2) is -0.92 MeV, they agree with the previous estimation of -1.83 MeV and -1.057 MeV, respectively.

For the J = 2 TBMEs, the off-diagonal term is -0.92, which is agreed with the -12\sqrt{7}/35 V_0 \bar{R} = -0.96 MeV. However, the diagonal terms are -1.06 and -0.71, which is a bit different from -24/35 V_0 \bar{R} = -0.74 and -6/5 V_0 \bar{R} = 1.27 . Nut we have to notice that the 2+ state is offset in previous calculation.

Binding energy, Mass, Separation energy, and effective interaction energy

Leave a comment

I hate to re-formulate the relations between them, so I write out in here.

The binding energy is negative.

\displaystyle M(^{Z}A) = (A-Z) m_n + Z m_p + BE(^{Z}A)

\displaystyle S_p(^{Z}A) =M(^{Z-1}A-1) + m_p - M(^{Z}A)

\displaystyle S_n(^{Z}A) =M(^{Z}A-1) + m_n - M(^{Z}A)

Separation energy in term of Binding energy

\displaystyle S_p(^{Z}A) =BE(^{Z-1}A-1) - BE(^{Z}A)

\displaystyle S_n(^{Z}A) =BE(^{Z}A-1) - BE(^{Z}A)

Binding energy and Mass in term of separation energy

\displaystyle BE(^{Z}A) = BE(^{Z-1}A) + S_n(^{Z}A+1) - S_p(^{Z}A+1)

\displaystyle BE(^{Z}A) = BE(^{Z+1}A) - S_n(^{Z+1}A+1) + S_p(^{Z+1}A+1)

\displaystyle M(^{Z}A) = M(^{Z-1}A) + S_n(^{Z}A+1) - S_p(^{Z}A+1) + m_p - m_n

Introduce p-n separation energy that remove a proton and a neutron,

\displaystyle S_{pn}(^{Z}A) = M(^{Z-1}A-2) + m_p + m_n - M(^{Z}A)

Also, the 2p and 2n separation energies

\displaystyle S_{2p}(^{Z}A) = M(^{Z-2}A-2) + 2m_p  - M(^{Z}A)

\displaystyle S_{2n}(^{Z}A) = M(^{Z}A-2) + 2m_n - M(^{Z}A)

in term of binding energy

\displaystyle S_{pn}(^{Z}A) = BE(^{Z-1}A-2) - BE(^{Z}A)

\displaystyle S_{2p}(^{Z}A) = BE(^{Z-2}A-2)   - BE(^{Z}A)

\displaystyle S_{2n}(^{Z}A) = BE(^{Z}A-2)  - BE(^{Z}A)

I am not sure how other people think, the separation energy is a more intuitive way to think of interaction energy than using binding energy.

For example, the energy of the d5/2 neutron in 17O is S_n, the energy for the two d5/2 neutrons in 18O is S_{2n}. See this post. Thus, the interaction energy of the two d5/2 neutrons in 18O is S_{2n} - 2 S_n . This is the extra energy caused by the n-n interaction.

Following the same idea, the interaction energy of the pn pair, is S_{pn} - S_p - S_n.

Wigner-Eckart theorem

Leave a comment

This is a theorem that connect a matrix element to a reduced matrix element that pull out the “magnetic sub-state”. We know that all quantum orbital with angular momentum j are degenerate in to 2j+1 sub-states. A complete wave function looks like this \left|\alpha j m \right > , where \alpha is other quantum numbers, and for convenient, we will drop it in later discussion.

A matrix element for a spherical tensor operator with rank k and q-th component is

\displaystyle \left<j m| (T_k)_q | j' m' \right>

And the m dependence can be pull out as

\displaystyle \left<j m| (T_k)_q | j' m' \right>  = (-1)^{j-m} \begin{pmatrix} j & k & j' \\ -m & q & m' \end{pmatrix} \left<j || T_k || j'\right>

Here we used the 3-j symbol, and the term \left<j || T_k || j'\right> is called the reduced matrix element. This stored the essential physical information, and the 3-j symbol represents the coupling between difference m-state, a kind of geometrical effect.

In nuclear physics, we often encounter the matrix element. In fact, the total Hamiltonian is a rank 0 tensor, magnetic quadruple is a rank 2 tensor. And the matrix representation of the Hamiltonian requires the calculation of matrix element.

For the total Hamiltonian, the 3-j symbol becomes

\displaystyle (-1)^{j-m} \begin{pmatrix} j & 0 & j' \\ -m & 0 & m' \end{pmatrix} =\\ (-1) ^{j -j'} \sqrt{\frac{(j-m)!(j'+m)!}{(j+j'+1)(j-j')!(j'-j)!(j-m)!(j'+m)!}},  m=m'

For j = j', it is particularly simple,

\displaystyle (-1)^{j-m} \begin{pmatrix} j & 0 & j \\ -m & 0 & m' \end{pmatrix} = \sqrt{\frac{1}{(2j+1)}}, m=m'

Thus, the Wigner-Eckart theorem becomes

\displaystyle \left<j m| H | j m' \right> = \sqrt{\frac{1}{(2j+1)}} \left<j || H || j\right>, m=m'

The advantage of the theorem is that, we can pick any convenient m-state to calculate the matrix element, and use the theorem to have the reduced matrix element, and then, we can have all other matrix elements with different m-states! for example, we can have the reduced matrix element from the m=0 state,

\displaystyle  \left<j || H || j\right> = \sqrt{2j+1} \left<j 0| H | j 0 \right>

than, put back the reduced matrix element, and we get.

\displaystyle \left<j m| H | j m \right> = \left<j 0| H | j 0 \right>

We can see, for rank 0 tensor, the diagonal elements are all degenerated, as we expected.

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.

Alpha cluster and alpha separation energy

Leave a comment

The alpha separation energy is the energy to add int a nucleus, so that it will break up into an alpha-particle and the rest of the nucleus. If we define the mass of a nucleus with mass number A and charge number Z as M(Z,A), the alpha-separation energy is,

\displaystyle S_\alpha(Z,A) = M(Z-2,A-4) + M(\alpha) - M(Z,A) .

The following plot the nuclides chart for alpha-separation energy,

The chart can be divided in 2 regions. In the upper region, nuclei have negative alpha-separation energies, i.e. nucleus gives out energy when emitting an alpha particle, thus, they are alpha-emitter. In the lower region, the alpha-separation energies are positive. And we can see that there are some local minimum for the S_\alpha.

According the alpha-decay theory, alpha-particle should be formed inside a nucleus before it tunnels through the nuclear potential and get out. the formation of the alpha-particle is described as the preformation factor.

It seems that, the alpha-separation energy, somehow, relates to the preformation factor.

The alpha cluster is studied for the many light nuclei, particularly on the 8Be, 12C, 16O, 20Ne, 24Mg, 28Si, 32S, 36Ar, 40Ca, and 44Ti. [Ref??]

From the above plot, it seems that the alpha-separation energies have no correlation with the magic number, and also no correlation with the Z=even. In a naive imagination, alpha-cluster could appear at all Z=even, A = 2Z nuclei. A trivial example is the the 8Be, its alpha-separation energy is -0.09 MeV and it will decay or split into 2 alpha particles. In 12C, the Hoyle state at 7.7 MeV is a triple-alpha state, note that S_\alpha(12C) = 7.37 MeV.

When we use the shell model to look at 8Be and 12C, we will found that the alpha-cluster is “not” making any sense. The protons and neutrons occupy the 0s1/2 and 0p3/2 orbitals. Since 9Be is stable and ground state spin is 3/2. A 8Be nucleus can be obtained by removing a 0p3/2 neutron. We may guess that, the 4 nucleons at the 0s1/2 orbital, which is an alpha-particle, somehow, escape from the 8Be, leaving the 4 0p3/2 nucleons behind. And the 4 0p3/2 nucleons “de-excite” back to the 0s1/2 shell and becomes another alpha-particle. But it seems that it does not make sense. Also, How to use shell model to describe the triplet-alpha cluster?

Above is the the alpha-separation energy for light nuclei. We can see that, there are few local minimum around 8Be, 20Ne, 40Ca, and 72Kr. Near 20Ne, the 18F, 19F, and 19Ne are having small alpha-separation energies, relative to the near by nuclei. For 20Ne, it could be understand why the alpha-separation energy is relatively smaller, as the 20Ne can be considered as 16O core with 4 nucleons in sd-shell. And 20Ne is well deformed that, there is a chance all 4 nucleons are in the 1s1/2 and form a quai-alpha particle. But the situation is a bit strange for 19Ne, 19F, and 18F. In order to form an alpha-cluster, or a quasi-alpha-particle, one or two nucleon from the p-shell has to excited to sd-shell. But the situation may be even more complicated. For some heavy alpha emitter, the s-orbital nucleons are tiny fraction compare to the rest of the nuclei, so, how the alpha-particle is formed before the decay is still unknown. This suggests that the involvement of s-orbital is not needed in alpha formation.

S_\alpha(19F) = 4.014 MeV, and there are many states around 4 MeV in 19F, I am wondering, one of these state is alpha cluster? If so, 19F(d,d’) reaction could excite those states and we will observed 15N + alpha. [need to check the data]. Similar experiment could be done on 18F, 19Ne, and 20Ne. If it is the case, that could provides some information on the alpha formation.

[need to check the present theory of alpha formation]

The brachistochrone curve

Leave a comment

I always want to do the calculation by myself, just for fun.

The problem is, given 2 points, A and B, with height difference H, and horizontal distance d, only under gravity and without friction, what is the shortest-time curve to connect the two points?

The problem is visualized in above figure. s(x) is the path length. The velocity depends on the height.

\displaystyle \frac{1}{2}mv^2 = mgh \Rightarrow v=  \sqrt{2gy}

The total travel time from point A to point B is

\displaystyle T = \int \frac{ds}{v(x)}

\displaystyle ds^2 = dx^2 + dy^2 \Rightarrow ds = dx \sqrt{1+ y'(x)^2}

put everything together

\displaystyle T = \int \frac{ds}{v(x)} = \int \sqrt{\frac{1+y'(x)^2}{2g y(x)}} dx = \int \frac{1}{\sqrt{2g}} L(x, y, y') dx

We want to minimize the total travel time by variate the path y(x) . This is exactly the same a the Lagrangian mechanics. The Lagrangian equation is

\displaystyle \frac{\partial L}{ \partial y} - \frac{d}{dx} \left( \frac{\partial L}{\partial y'} \right) = 0

Let’s compute the partial derivatives,

\displaystyle \frac{\partial L}{ \partial y} = - \frac{(1+y'^2)^{1/2}}{2 y^{3/2}}

\displaystyle \frac{\partial L}{ \partial y'} = - \frac{y'}{\sqrt{y}(1+y'^2)^{1/2}} =K

The total derivative is

\displaystyle \frac{dK}{dx} = \frac{\partial K}{ \partial y} \frac{dy}{dx} + \frac{\partial K}{ \partial y'} \frac{dy'}{dx} + \frac{\partial K}{ \partial x}

\displaystyle \frac{\partial K}{ \partial y} = - \frac{y'}{2 y^{3/2}(1+y'^2)^{1/2}}

\displaystyle \frac{\partial K}{ \partial y'} = \frac{1}{y^{1/2}(1+y'^2)^{3/2}}

\displaystyle \frac{dK}{dx} = \frac{2y y'' - y'^2 -y'^4}{2y^{3/2}(1+y'^2)^{3/2}}


\displaystyle - \frac{(1+y'^2)^{1/2}}{2 y^{3/2}} = \frac{2y y'' - y'^2 -y'^4}{2y^{3/2}(1+y'^2)^{3/2}} \Rightarrow  -(1+y'^2)^2 = 2y y'' - y'^2 -y'^4

\displaystyle  2y y'' + y'^2 +  1 = 0 \Rightarrow \frac{1}{y'}\frac{d}{dx}\left( y(1+y'^2) \right) = 0


\displaystyle y(1+y'^2)  = C

From here, I am not sure how to get to the cycloid. For a cycloid drawn from a circle with radius R, a downward cycloid is

\displaystyle x(\theta) = R( \theta - \sin\theta ) \\ y(\theta) = R(\cos\theta -1 )

Lets check is the cycloid fulfill the requirement.

\displaystyle \frac{dy}{dx} = \frac{dy}{d\theta}\frac{d\theta}{dx} = \frac{-\sin\theta}{1-\cos\theta}

\displaystyle y(1+y'^2) = R(\cos\theta -1) \left( 1 + \frac{\sin^2\theta}{(1-\cos\theta)^2} \right) = 2R =C

So, the cycloid is the solution. Below is a cycloid with 1 unit of radius.

We now knew that the cycloid is the solution, but we have to find the constant C, or the radius R, so that the curve pass through the points A and B.

Fixing point A at (0,0), the point B at (d, -H) and it is on the cycloid. Thus,

\displaystyle d = R( \theta_0 - \sin\theta_0 ) \\ -H = R(\cos\theta_0 -1 )

Solve for \theta_0 , combine the equations,

\displaystyle \frac{d}{H} = \frac{ \theta_0 - \sin\theta_0 }{1 - \cos\theta_0 } \Rightarrow d - H \theta_0 = d \cos\theta_0 - H \sin\theta_0

Set d = G \cos\phi, H = G \sin\phi,  G = \sqrt{d^2 + H^2}, \tan\phi = H/d,

\displaystyle  d - H \theta_0 = G \cos(\phi+\theta_0),  0 < \theta_0 \leq 2\pi

There may be no analytical solution. The solution is the intercept between the 2 curves below.

In the above example, d = 10, H = 5, and \theta_0 = 3.50837, R = 2.586 . The brachistochrone curve look like this:

Since we know the solution now, so, what is the minimum travel time from point A to point B?

\displaystyle ds =d\theta R \sqrt{ \left(\frac{dx}{d\theta}\right)^2 + \left(\frac{dy}{d\theta}\right)^2} = d\theta \sqrt{2}\sqrt{1-\cos\theta} = \sqrt{2R} \sqrt{y} d\theta

\displaystyle \frac{ds}{v} = \sqrt{2R} \frac{\sqrt{y}}{\sqrt{2g y}} d\theta = \sqrt{R} d\theta

which is a constant of \theta !! So the travel time is

\displaystyle T = \sqrt{\frac{R}{g}} \theta_0

Also, the motion in term of time is

\displaystyle x(t) = R \left( \sqrt{\frac{g}{R}}t - \sin\left(\sqrt{\frac{g}{R}} t\right) \right) \\ y(t) = R\left(\cos\left(\sqrt{\frac{g}{R}}t\right) -1 \right)

Now, Lets us investigate the travel time for a straight line from point A to point B. The path is

\displaystyle x = -\frac{d}{H} y

The travel time is

\displaystyle t = \int \frac{ds}{\sqrt{2gy}} = \int \frac{\sqrt{H^2+d^2}}{H\sqrt{2gy}} dy = \frac{\sqrt{H^2+d^2}}{H\sqrt{g}} \sqrt{2y}

Thus, the coordinate in term of time is

\displaystyle x(t) = \frac{1}{\sqrt{2}} \frac{d H \sqrt{g}}{H^2+d^2} t^2 \\ y(t) = - \frac{1}{\sqrt{2}} \frac{H \sqrt{g}}{H^2+d^2}  t^2

Next, we also check a curve that, if H > d , then, the path go vertical down by H-d, than do a circular path with radius d. If H < d, then do a circular path with radius H, and a horizontal path to point B.

Lets only study the case when d > H .

The circular path motion is

\displaystyle x(\phi) = H - H \cos\phi \\ y(\phi) = -H \sin\phi,  \phi = [0, \pi/2]

the path length is s(\phi) = H \phi

\displaystyle t_1 = \int \frac{ds}{\sqrt{2gy}} = \int \frac{H}{\sqrt{2g H \sin\phi}}  d\phi = \sqrt{\frac{H}{2g}} \left(2.6221 - 2 F\left( \frac{1}{2} \left( \frac{\pi}{2} - \phi \right) , 2 \right) \right)

where the $latex F(x, m) is the elliptic integral of the 1st kind.

The rest of the path takes time

\displaystyle t_2 = \frac{d-H}{\sqrt{2gH}}

Here is the comparison to all 3 paths

Factorio transport belt problem

Leave a comment

Factorio is a logistic game. very addictive.

Inside, the transport belt is the main way to transport materials from one place to another place. For massive transport, 2 or more belts will be used in parallel.

There is also a splitter, that can divide the belt into 2 belts and distribute the load evenly.

A problem is, suppose we have n-belts and different load, how can we use splitters to evenly distribute the load into m-belts?

To solve the problem, lets look at the mathematics for a splitter. Suppose we have 2-belts-in and 2-belts-out, we can place a splitter in the middle. The load for the in-put is a, b.

If we define the load as a vector, the splitter is a matrix that transform the load vector.

\displaystyle  \begin{pmatrix} \frac{a+b}{2} \\ \frac{a+b}{2}  \end{pmatrix} = \frac{1}{2}\begin{pmatrix} 1  & 1 \\ 1 & 1 \end{pmatrix}\begin{pmatrix} a \\ b \end{pmatrix}

Since the total number of load must be the same before and after a splitter, thus, the sum of the element in each row must be 1 for a splitter matrix.

If we have n-belts, a splitter is placed between the i-th and j-th belt. The splitter matrix is

\displaystyle S_{ij} = \begin{pmatrix} 1 & 0 & ... & ... & ... & ... & 0 \\ ... \\ 0 & 0 & 1/2 & ... & 1/2 & ... & 0  \\ ... \\ 0 & 0 & 1/2 & ... & 1/2 & ... & 0 \\ ... \\ 0 & ... & ... & ... & ... & ... & 1 \end{pmatrix}

For example, we have 3 belts, and if we want to evenly distribute among them, we try putting 2 splitters like this.

The splitter matrix are

\displaystyle S_{12} = \frac{1}{2} \begin{pmatrix} 1 & 1 & 0 \\ 1 & 1 & 0 \\ 0 & 0 & 2 \end{pmatrix}, S_{23} = \frac{1}{2} \begin{pmatrix} 2 & 0 & 0 \\ 0 & 1 & 1  \\ 0 & 1 & 1 \end{pmatrix}

The combined matrix is

\displaystyle S = S_{23}S_{12} =  \frac{1}{4} \begin{pmatrix} 2 & 2 & 0 \\ 1 & 1 & 2 \\ 1 & 1 & 2 \end{pmatrix}

Clearly, it is far from a evenly distributed matrix, that all element are \frac{1}{3}.

Lets try to use this black for k-times, but before that, to compute power of matrix, it is better to diagonalize the matrix,

\displaystyle S = P\cdot Q\cdot P^{-1}, P = \begin{pmatrix} 1 & -2 & -1 \\ 1 & 1 & 1 \\ 1 & 1 & 0 \end{pmatrix}, Q = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1/4 & 0 \\ 0 & 0 & 0 \end{pmatrix}


\displaystyle P^{-1} =\frac{1}{3} \begin{pmatrix} 1 &1 & 1 \\ -1 & -1 & 2 \\ 0 & 3 & -3 \end{pmatrix}


\displaystyle S^k= P \cdot Q^k \cdot P^{-1} = \frac{1}{3} \begin{pmatrix} 1 +  2^{1-2k} & 1+ 2^{1-2k} & 1 - 2^{1-2k} \\ 1 - 4^{-k} & 1 - 4^{-k} & 1 + 2^{1-2k} \\ 1 - 4^{-k} & 1 - 4^{-k} & 1 + 2^{1-2k} \end{pmatrix}

We can see that, k \rightarrow \infty, the matrix element goes to 1/3. Thus, in principle, we need infinity block to evenly distribute the loading for 3 belts. In practice, 3 pairs of splitters, setting the input load vector be {1-x-y, x, y}, the output load vector is { 0.344 - 0.03 y, 0.328 + 0.02 y , 0.328 + 0.02 y }. Good enough.

Can we using a similar idea to bull-force the problem for 4-belts, or n-belts in general?

It turns out is yes, for an infinite block, we can evenly distribute the load.

How about for finite splitters? For 4-belts, lets try this setting

the last splitter connect the 1st belt and the 4-belt. The matrix S = S_{14} \cdot S_{23} \cdot S_{34} \cdot S_{12} has the matrix element to be 1/4, exactly!! In the community of factorio player, they found this setting long time ago. The idea behind is, first mixing the 1-2 belts, the 3-4 belts, at this point the loads of 1-2 are the same and that for 3-4 are the same too. next, mix the 2-3 belts, and 1-4 belts will do the job. Following this idea, the 2^n belts can be done.

For 6-belts, we can separate them into 3 pairs, and mix the 3 pairs. But we haven’t find a finite splitters setting for mixing 3 loads.

To investigate this problem more, lets study the algebra for the splitter matrix. We can easily found these porperties

\displaystyle S_{ij} = S_{ji}

\displaystyle S_{ij} \cdot S_{lm} = S_{lm} \cdot S_{ij} for i,j, l, m are all different

\displaystyle S_{ij} \cdot S_{jk} \neq  S_{kj} \cdot S_{ji}

I guess there are somethings interesting for a ring-setting , for example, for 4-belts, S_{12} \cdot S_{23} \cdot S_{34}. But I cannot figure it out. I guess there is some kinds of group structure for the splitter matrix, could be an isomorphism to 2-body exchange group ?

Goodness of Fit

Leave a comment

In this post, we discussed the Goodness of Fit for a simple regression problem.

In nuclear physics, the data we collected are plotted in histogram, and we have to fit the histogram to find the energies and resolutions. This is another class of problem for Goodness of Fit.

In general, depends on how we treat the data,

  • treat it as a regression problem with error for each data follows Poisson distribution.
  • treat it as a testing hypothesis that what is the chance for the experimental distribution be observed, given that the model distribution is truth?

For the regression problem, there are 2 mains metrics, one is the R-squared, and the another is the chi-squared.

For the testing hypothesis, the experimental data is tested against the null hypothesis, which is the model distribution. A common method is the Pearson’s chi-squared test.

Note that the chi-squared and chi-squared test are two different things, although they are very similar in the algorithm.

Another note is that, please keep in mind that there are many things not discussed in below, there are many details and assumptions behind.

The R-squared is very simple, the meaning is How good the model captures the data. In a standard ANOVA (Analysis of variance, can be imagine as the partition of the SSTotal ), we have

\displaystyle SSTotal = \sum_{i} (Y_i - \bar{Y} )^2 = \left( \sum_i Y_i^2 \right) + \bar{Y}^2

\displaystyle SSR = \sum_{i} (Y_i - \hat{Y_i})^2

The SSRTotal is the sum of square of “residuals” against the mean value. i.e., If the data has no feature, no peak, the SSTotal will be minimum. SSTotal represents the “size” of the feature in the data. Larger the SSTotal, the data is more different from the mean values, and more feature-rich.

The SSR is the sum of square of residuals against the fitted function. If the model captures the data very well, SSR will be minimum. Since the data could have some natural fluctuation, the MSR = SSR/ndf represents the unbiased estimator of the “sample variance”, where ndf is the number of degree of freedom.

n is the size of the data, p is the number of parameter of the model, or the size of the model. The SSTotal is the “size” of feature, SSR is the “size” of the things that did not captured by the model. The basic of ANOVA is the study of the partition of the data.

The difference between SSTotal and SSR is the “size” of the model captured data.

\displaystyle R^2 = \frac{SSTotal - SSR}{SSTotal}  < 1

The idea of R-squared is that, the numerator is the size that the model captured from the data, and the denominator is the size of the feature in the data. If the data is featureless, thus the data is mainly from the natural fluctuation, therefore, there is nothing a model can be captured. so, R-squared will be small.

Chi-squared includes the sample variance,

\displaystyle \chi^2 = \sum_i \frac{(Y_i - \hat{Y_i})^2}{\sigma_i^2}

We can see that the \chi^2 is (roughly) the SSR divided by sample variance \sigma_i^2 . As we discussed before, the SSR is the size of the natural fluctuation — sample variance. If we divided the SSR with the sample variance, the \frac{(Y_i - \hat{Y_i})^2}{\sigma_i^2} \rightarrow 1 when the fitting is good. It is because, if the fitting is not good, the SSR will be larger than the sample variance (as there are something the model did not capture), or, if the fitting is too good, the SSR will be much smaller than the sample variance.

For histogram data, the sample variance of each bin should follow the Poisson statistics that \sigma_i^2 = Y_i .

The reduced chi-squared is

\displaystyle \bar{\chi^2} = \frac{\chi^2}{ndf}

For n data points, they are living in the “n-dimensional” space. Suppose there is a model to generate these data points, and the model has p parameters. Thus, the model “occupied” “p-dimensional” space. The rest “(n-p)-dimensional” space is generated by the natural fluctuation. Thus, the “size” of the \chi^2 should be (n-p) or the degree of freedom. For a good fit

\displaystyle \bar{\chi^2} = 1

Pearson’s chi-squared test, is a testing, asking, does the data distribution can be captured by a model distribution, more technically, what is the probability to observed the experimental distribution, given that the model distribution is truth. One limitation for using this test is that, the number of data for each bin cannot be too small.

The test calculate the following quantity,

\displaystyle X^2 = \sum_i \frac{(Y_i - \hat{Y_i})^2}{\hat{Y_i}}

The X^2 follows the \chi^2-distribution with degree of freedom ndf.

As all testing hypothesis, we need to know the probability of the tail. For \chi^2-distribution, only 1-tail is available, and this is usually called the p-value. Here shows an example of the \chi^2-distribution.

The meaning of the p-value is that,

given the null hypothesis, or the model distribution, is truth, the probability of the experimental data happen.

When the p-value is less than 5%, the null hypothesis should be rejected. In other word, the model distribution cannot represents the experimental distribution. Or, if the model distribution is truth, less than 5% chance to observe the experimental data.

Postdoc seminar

Leave a comment

Haven’t had an hour long presentation for a while, I just gave one. Can see from the web site.


Annotation 2020-08-06 141821.png

Annotation 2020-08-06 141951

Rotation of deformed nucleus

Leave a comment

The total Hamiltonian of a nuclear system must be rotational invariant, because the state should not have prefer direction. Thus, the total spin J of the wave function is a good quantum number.

There is a tricky thing should point out here. The spherical harmonic clearly has prefer direction, it rotates around the “z-axis”. Since this “z-axis” is arbitrary, the wave function (or the spherical harmonic) still have no spacial preference and can point to any direction. i.e., the energy is the same for all direction. The nuclear spin also the same. Nuclear spin is a mini magnetic dipole that is direction. The nuclear spin wave function is non-rotation symmetry but its energy is the same for any direction under a Hamiltonian without magnetic field.

For axial deformed J = j = 0 nucleus, the nucleus is not rotated at the ground state, and can orient to any direction. When this nucleus rotate, it has to rotate perpendicular to the body axis, because it is the only principle axis that exhibit a rotation. And since the rotation should preserved parity (?), the rotation increase the total spin by 2 every step, i.e. J = R + j = R

For axial deform nucleus with j > 0 , the ground state has to rotate so that the total spin J = R + j is constant of motion.In Nilsson model, the intrinsic angular momentum is not a good quantum number, in other words, a Nilsson orbital is a superposition of many difference spin j. And the spin j is precessing around the body (symmetry) axis. Thus, the introduction of rotation R is needed to make the total spin to be constant of motion. ( I am not sure this explanation is correct, it seems very artificial, and I don’t know how to construct the Hamiltonian)

In previous paragraph, there is some exceptional cases for the maximum \Omega orbital, for example, the 3/2[101] is a pure p3/2. For such orbital, the Nilsson orbital only contains one single j. In that case, the state for J = j need not be rotate.

Annotation 2020-07-08 222954.png

In above picture, the left illustrates the relation between difference quantum numbers. on the right, it is more interesting. As we mansion before, the spin j is precessing around the body axis. And the rotation R must change in both direction and magnitude with the spin j to keep the total spin J to be constant. So, the rotation of a deformed nucleus can be very complicated for j > 0 . One interesting and simple case is that the J = 1/2, which is align on the body axis. In that case, the rotation direction R is also rotating with a constant magnitude.  Nevertheless, a Nilsson orbital usually contains many j and each of them coupled to different rotation to form a constant J. Therefore, when a deformed nucleus start to rotate, many rotation are happening at the same time.

Because of this complicated motion ( although it is quite a classical point of view ), the nucleon is also moving in this rotating body, so that the Coriolis force will be there. This will be the next topic.

Older Entries