Nuclear Physics 101

Note on lifetime, decay width, and gamma-ray transition probability

May 2, 2024

The Gamma Transition is on this post, but it seems that the equation connecting transition rate $T$ and reduced transition rate $B(qL)$ missed a factor

The relations between mean lifetime $\tau$ , half-life $T_{1/2}$ , (partial) decay width $\Gamma$ , and transition rate $T$ are:

$\displaystyle \tau = \frac{T_{1/2}}{\ln2} = \frac{\hbar}{\Gamma} = \frac{1}{T}, ~ \hbar = 6.582119569\times 10^{-16} \textrm{eV}\cdot\textrm{s}$

The gamma transition rate $T_\gamma$ in s^-1 is related to the electric reduced transition rate $B(EL)$ in e²fm^2L is:

$\displaystyle T_\gamma = K(L) E_\gamma^{2L+1} B(EL)$

$\displaystyle K(L) = \frac{8\pi \alpha c (L+1) }{ L ((2L+1)!!)^2} \left( \frac{10^6}{\hbar c}\right)^{2L+1}$

$\displaystyle \alpha = \frac{e^2}{4\pi\epsilon_0 \hbar c} \sim \frac{1}{137.06}$

$\displaystyle c = 299792458 \times 10^{15} \textrm{fm/s}$

Here is a table for the first few coefficients $K(L)$

Multiple	Decay rate (1/s), $E_\gamma$ in MeV
E1	$1.59021 \times 10^{15} E_\gamma^3$ B(E1) [e²fm²]
E2	$1.22519 \times 10^{9} E_\gamma^5$ B(E2) [e²fm⁴]
E3	$5.70796 \times 10^{2} E_\gamma^7$ B(E3) [e²fm⁶]
E4	$1.69666 \times 10^{-4} E_\gamma^9$ B(E4) [e²fm⁸]
E5	$3.45706 \times 10^{-11} E_\gamma^{11}$ B(E5) [e²fm¹⁰]

The reduced decay rate and reduced excite rate is related by:

$\displaystyle B(qL; I_1 \rightarrow I_2 ) = \frac{2I_2+1}{2I_1+1} B(qL; l_2 \rightarrow l_1)$

For example, in ¹⁰C, the first 2⁺ (3.354 MeV) to 0⁺ (g.s.) has mean lifetime of 219 fs (E. A. McCutchan et al., PRC 86, 014312 (2012)), the transition rate is then $4.566 \times 10^{12}$ s^-1, and the decay width is 3.01 meV. As this is 100% E2 decay, using the above formula, the B(E2) is 8.781 e²fm⁴. And the excite B(E2; 0⁺ -> 2⁺) is 5 times larger.

Poisson process and exponential decay

March 17, 2024

GoLuckyRyan Basic Leave a comment

How does a Poisson process generate the exponential decay?

The exponential decay is derived by

$\displaystyle \frac{dN}{dt} = - k N \rightarrow N(t) = N_0 \exp\left(-kt\right) = N_0 2^{- t / \tau}, ~ k = \frac{\log(2)}{\tau}$

which means the rate of the loss of the number of nuclei is proportional to the number of nuclei. and $\tau$ is the half-life.

This is a phenomenological and macroscopic equation that tells us nothing about an individual decay.

Can we have a microscopic derivation that, assumes the decay of an individual nucleus following a random distribution?

From this post, we derived the distribution for the differences in a list of sequential random numbers. The key idea is the probability of having a decay at exactly time $t - \delta t$ to $t$ , denote as $f(t, \delta t)$ , is

$\displaystyle f(t, \delta t) = P(t-\delta t, 0)P(\delta t, 1)$

where $P(t,k)$ is the probability of having $k$ decay within (time) length $t$ , which is a Poisson distribution.

$\displaystyle P(t, k) = \frac{t^k}{k!} \exp(-t)$

$\displaystyle f(t, \delta t) = \frac{(t-\delta t)^0}{0!} \exp(-(t-\delta t)) \frac{\delta t}{1!} \exp(-\delta t) = \exp(-t) \delta t$

Now, the total probability of all possible decay within time $T$ , assume not more than 1 decay can happen in $\delta t$ ,

$\displaystyle F(T) = \sum_{\delta t \rightarrow 0}^{t < T} f(t, \delta t) = \int_0^{T} \exp(-t) dt = 1- \exp(-T)$

Thus, the number of nucleus not decay or survived within time $T$ is

$\displaystyle N(t) = N_0 (1 - F(t)) = N_0 \exp(-t) = N_0 2^{-t / \tau}$

Radioactivity and Number of Isotopes

March 16, 2024

GoLuckyRyan Basic Leave a comment

Suppose the radioisotope has half-life of $\tau$ and initial number of $N_0$ , the number of isotopes at give time is

$\displaystyle N(t) = N_0 2^{-\frac{t}{\tau}}$

The number of decay at a given time, measured with a time period $\Delta t$ is

$\displaystyle D(t, \Delta t) = N(t - \Delta t) - N(t)$

The rate of number of decay at a given time is

$\displaystyle R(t) = \lim_{\Delta t \rightarrow 0} \frac{N(t- \Delta t) - N(t)}{\Delta t} = - N'(t) = \frac{\log2}{\tau} N_0 2^{-\frac{t}{\tau}} = \frac{\log2}{\tau} N(t)$

$\displaystyle N(t) = \frac{\tau}{\log2} R(t)$

This, the rate of number of decay is always proportional to the number of radioactive isotopes, therefore, by measuring the radioactivity, the number of radioactive isotopes can be deduced.

For example, suppose we have a sample of tritium with half live of $3.89 \times 10^8$ sec. When we measure the sample, we got 100 beta decay in 1 min, that means, the radioactivity is 1.67 decay per second. Therefore the number of isotopes is $9.14 \times 10^8$ tritium. Notice that, this number is the number of isotopes at the end of measurement.

For very short half-life, say, few mili-sec, if we collect the number of decay in a few seconds, the number of isotopes changed a lot and the simple ratio $\frac{\tau}{\log2}$ is no longer accurate.

$\displaystyle \frac{D(t, \Delta t)}{N(t)} = \frac{N(t - \Delta t) - N(t)}{N(t)} = 2^{\frac{\Delta t}{\tau}} - 1$

and

$\displaystyle \lim_{\tau \rightarrow \infty} \left( 2^{\frac{\Delta t}{\tau}} - 1 \right) \frac{\tau}{\log 2} = 1$

Transform Hamiltonian in rotating frame

March 15, 2024

GoLuckyRyan Basic Leave a comment

I revisit the NMR theory and technical stuff. Doing a 90 degrees flip and measure the Free Induction Decay (FID) signal. The NMR theory is simple, under a magnetic field, the proton spin is precessing at the Larmor frequency, but the spin aligns with the field (alone the z-axis) and the change of the magnetic field on the x-y plane is very weak. To increase the change of the magnetic field, we can flip the spin 90 degrees, make it precesses on the x-y plane and maximize the change of the magnetic field. With the change of the magnetic field induct electric field and voltage on a coil (on the x-axis for example), an signal can be detected at the Larmor frequency. Since the spin will graduate align itself back to the z-axis and so the signal will decay. The technique to flip the spin is called Nuclear Magnetic Resonance. I made a very clumsy note on that many years ago in this pdf.

In order to flip the spin, we need to apply a rotating field with frequency $\Omega$ , in the rotating frame created by the rotating field, an inducted field is created alone the negative z-axis as the world is rotating around (in an opposite direction), and also, a constant field points along the rotating x-axis. When the induced field in the rotating frame can cancel the external magnetic field, the over all effective magnetic field would be the one along the rotating x-axis. Since nuclear spin was on the z-axis, and since switch on of the rotating field, it feels a net magnetic field along the rotating x-axis and start to precess around it. If the rotating field switch off at the time that the spin precesses on the x-y plan, a 90 degree flip is done.

Mathematically, the Hamiltonian a spin under and external magnetic field is

$\displaystyle H_0 = -\gamma \vec{S} \cdot \vec{B} = -\gamma B S_z = \Omega_I S_z$

where $\gamma = \frac{e g_p}{2 m_p} = g_p \frac{\mu_p}{ \hbar } = 42.577~\textrm{MHz}/\textrm{T}$ is the gyromagnetic ratio of proton, in which $g_p$ is the g-factor of proton and $\mu_p$ is the proton magneton. The Larmor frequency is $\Omega_I = -\gamma B$ .

with a oscillating magnetic field (which is effectively combination of a right-hand and left-hand rotating field), the Hamiltonian is

$\displaystyle H = \Omega_I S_z + 2k \cos(\Omega t) S_x$

The factor $k$ represent the strength of the rotating field. The factor of 2 will comes in handy later. In order to change to the rotating frame, suppose the Lab wavefunction is $\phi$ , the rotating frame wavefunction is $\phi_R = \exp(i\Omega t S_z) \phi$ , the rotating frame Hamiltonian must satisfy

$\displaystyle i \frac{d}{dt} \phi_R = i \frac{d}{dt} \left( \exp(i\Omega t S_z) \phi \right) = - \Omega S_z \exp(i\Omega t S_z) \phi + \exp(i\Omega t S_z) \left( i \frac{d}{dt} \phi \right) \\ = - \Omega S_z \phi_R + \exp(i\Omega t S_z) H \exp(-i\Omega t S_z) \phi_R \\ = \left( - \Omega S_z + \exp(i\Omega t S_z) H \exp(-i\Omega t S_z) \right) \phi_R = H_R \phi_R$

Thus,

$\displaystyle H_R = ( \Omega_I - \Omega) S_z + 2k \cos(\Omega t) \exp(i\Omega t S_z) S_x\exp(-i\Omega t S_z) \\=(\Omega_I - \Omega) S_z + 2k \cos(\Omega t) ( \cos(\Omega t) S_x - \sin(\Omega t )S_y) \\=(\Omega_I - \Omega) S_z + k S_x + k \left( \cos(2 \Omega t) S_x - \sin(2\Omega t )S_y \right) \\= (\Omega_I - \Omega) S_z + k S_x$

We can drop the last time dependence terms, because it oscillate at twice of the frequency of the rotating frame, that, in average, the contribution is averaged out. And the Hamiltonian is as expected as our early discussion. an induced field $- \Omega S_z$ along the z-axis, and constant field along the x-axis.

The strength of the oscillating field $k$ has a meaning of rotational frequency. If we apply the field so that $T = \pi /4 /k$ and $\Omega = \Omega_I$ , this will rotate the spin 90 degree and give a maximum reading from the NRM coil. The challenge is that, the actual value of $k = -\gamma B_{osc}$ depends on the effective oscillating magnetic field strength $B_{osc}$ , which depends on the power transfer from the source to the NMR coil, the impedance matching from the coil to the source, also the resonance frequency of the chamber. Some these factors are hard to measure and often, the 90 degree spin flip is achieved by trials and errors until the maximum FID is obtained.

GEANT4 installation on MacOS 13.3

February 15, 2024

GoLuckyRyan Basic Leave a comment

The installation in Ubuntu 22.04 can be found at the end of this post. The trick for the installation is enabling visualization. The one I pick is the Qt OpenGL, and it requires to have Qt installation.

To install Qt, I download the Qt installer from the Qt website. In the installer, I pick the manual components to install. The components are

Qt > Qt 6.6.2 > macOS
Qt > Qt 6.6.2 > Qt Quick 3D
Qt > Qt 6.6.2 > Qt Shader Tools
Qt > Qt 6.6.2 > Qt Quick Timeline
Qt > Qt 6.6.2 > Additional Libraries > Qt 3D
Qt > Qt 6.6.2 > Additional Libraries > Quick: 3D Physics
Qt > Developer and Designer Tools > CMake 3.27.7
Qt > Developer and Designer Tools > Ninja 1.10.2

I bet some of the components are not needed. Once installed (I installed at /Users/ryantang/Qt/) , we need to export

export Qt6_DIR=/Users/ryantang/Qt/6.6.2/macos/lib/
export Qt6OpenGLWidgets_DIR=/Users/ryantang/Qt/6.6.2/macos/lib/cmake/Qt6OpenGLWidgets/

Those two line let the cmake know where to find the Qt.

Now, I downloaded the GEANT4 source code

>cd Downloads
>wget https://gitlab.cern.ch/geant4/geant4/-/archive/v11.2.0/geant4-v11.2.0.tar.gz
>tar xzf geant4-v11.2.0.tar.gz
>mkdir geant4-v11.2.0_build
>cd geant4-v11.2.0_build
>cmake -DGEANT4_INSTALL_DATA=On -DCMAKE_INSTALL_PREFIX=/opt/geant4-v11.2.0 -DGEANT4_USE_OPENGL_X11=ON -DGEANT4_USE_QT=ON -DGEANT4_USE_QT_QT6=ON ../geant4-v11.2.0
>make -j10
>sudo make install
>echo 'source /opt/geant4-v11.2.0/bin/geant4.sh' >> ~/.bashrc

For test, we can use the example/basic/B1, or I use the Clarion2 simulation (https://fsunuc.physics.fsu.edu/wiki/index.php/Clarion2). Here is the screen I get.

Some thoughts on the quenching of spectroscopic factor (II)

January 1, 2024

GoLuckyRyan Basic spectroscopic factor Leave a comment

The nature of the spectroscopic factor is discussed here (from the point of view of mean-field) and here (from an experimental point of view). Here are some updates and additional discussions on the topic.

Up to now, I still believe that the quenching of the SF simply reflects the mismatch between experiment and theory, due to experimental conditions that limited configurations can be observed (so that the sum of SF is less than 1), and also the reaction theory does not take into account the correlation (so that the SF of each state is quenched). If both limitations are removed, there is no quenching. Although the SF is quenched, careful normalization gives us a great deal about the nuclear structure.

Is the spectroscopic factor or the cross-section quenched?

It seems that since the spectroscopic factor and the cross-section are closely related by

$\displaystyle \sigma_{exp} = C^2S \sigma_{th}$

so, the question is meaningless, because when SF is quenched from unity, the cross-section is also quenched from the therotical cross-section. However, since cross-section is an observable while SF is not, people consider that it is more proper to say quenching of cross-section.

Whether SF is an observable or not, an article Spectroscopic Factors: Observability and Measurability by L. D. Blokhintsev stated that “Spectroscopic factors are non-invariant under the unitary transformations of nuclear forces“, and Non-observability of Spectroscopic Factors by B. K. Jennings said that “At the level of microscopic calculations, they can be varied by unitary transformations and we expect hard and soft nucleon-nucleon potentials to generate quite different values for them. At the phenomenological level, the spectroscopic factor can be varied by field redefinitions.“. Thus, the spectroscopic factor is a pure construction that is fixed on a specific basis and interaction. Another article Nonobservable nature of the nuclear shell structure: Meaning, illustrations, and consequences by T. Duguet, H. Hergert, J. D. Holt, and V. Somà is worth reading.

However, I think SF is an observable because there should be a natural basis and fixed interaction. when the basis is fixed, and the 1st order of the nuclear interaction is well-known, the SF is well-defined, just like the coefficient of fractional percentage. Another argument is that there are many articles on extracting the spectroscopic factors, and although the uncertainty can be as large as 50% due to various experimental conditions and analysis methods, the community has a consensus on the quenching. That is, the community reached a consensus on the basis and range of interaction. I think that means, we are reaching the natural basis and 1st order understanding of the interaction.

The situation is like if we have the theoretical framework but lack some fundamental parameters fixed. The theory can produce all kinds of phenomena. Say, neutrino mass, without that, we don’t know the neutrino mass hierarchy. But we cannot say the neutrino mass hierarchy should fall in either the normal or inverted hierarchy and therefore neutrino mass hierarchy is not well-defined. My point is once the nuclear interaction and the basis are fixed, the SF is uniquely determined. So, some people would say, SF is an quasi-observable and that it is not an observable yet, but it will. May be, we should call SF is an extractable/deductible, which is a extractable/deductible quantity with some will-be-fixed parameter sets.

The last comment is that, if SF is not an observable, and is a pure arbitrary construct, all the works on it become meaningless. But, I believe SF is not aether, it reveals the structure and nucleon configuration of a nucleus. Also, the effective single-particle energy (ESPE) is an SF weighted energy of an orbital. If SF is arbitrarily constructed and can be anything, so the ESPE. but the ESPE very much agrees with the energy levels from Woods-Saxon potential which is a 1st order approximation of the nuclear mean field. A similar thing to the asymptotic normalization coefficient (ANC), ANC is a ratio of the wavefunction to the Coulomb wave function outside the nuclear potential, and the wavefunction is not an observable, so ANC is also not. But the wavefunction is surely can be deduced/extracted from the experimental data by fixed some nuclear parameters.

So, if SF is an observable, it is meaningless to ask is SF quenched or cross-section quenched.

If, the theoretical cross-section is perfectly calculated, the SF would be unity and no quenching for either SF or cross-section, and the nucleon configuration is embedded into the coefficient of fractional percentage.

Intersect between 2 line segments

December 30, 2023

GoLuckyRyan Basic Leave a comment

The whole discussion is on two-dimensional only. Suppose the 2 line segments are defined by points $(\vec{a_0}, \vec{a_1})$ and $(\vec{p_0}, \vec{p_1})$ . The 2 segments can either be parallel or non-parallel.

Non-parallel case. There are at least 2 methods to determine, straightforward thinking is to construct the equation of the 2 lines using the 2 segments

$\displaystyle \vec{r} = \vec{a_0} + h (\vec{a_1} - \vec{a_0}) \\ \vec{r} = \vec{p_0} + k (\vec{p_1} - \vec{p_0})$

solve for the parameters $h, k$ . One way is using bull force to solve it. Another way is the area ratio:

$\displaystyle h = \frac{ (\vec{p_0} - \vec{a_0}) \times (\vec{p_1} - \vec{a_0}) } { (\vec{p_0} - \vec{p_1}) \times (\vec{a_1} - \vec{a_0}) }$

The nominator is the area of the triangle $\vec{a_0}, \vec{p_0}, \vec{p_1}$ . The denominator is the area of the whole area $\vec{a_0}, \vec{p_0}, \vec{a_1}, \vec{p_1}$ . The whole area is

$\displaystyle A= (\vec{p_0} - \vec{a_0}) \times (\vec{a_1} - \vec{a_0}) +(\vec{a_1} - \vec{a_0}) \times (\vec{p_1} - \vec{a_0}) = (\vec{p_0} - \vec{p_1}) \times (\vec{a_1} - \vec{a_0})$

The simplification uses the identity $\vec{u}\times \vec{v} + \vec{v} \times \vec{w} = (\vec{u} - \vec{w} ) \times \vec{v}$ .Similarly.

$\displaystyle k = \frac{ (\vec{a_1} - \vec{p_0}) \times (\vec{a_0} - \vec{p_0}) } { (\vec{p_0} - \vec{p_1}) \times (\vec{a_1} - \vec{a_0}) }$

The intersect condition is

$\displaystyle 0 \le h \le 1~ \textrm{and} ~0 \le k \le 1$

In the case of parallel, there are 2 cases: with offset or without offset

In the case of offset, the “whole” area $A = 0$ , geometrically, the cross-product $(\vec{p_0} - \vec{a_0}) \times (\vec{a_1} - \vec{a_0}) < 0$ as it is clockwise. Since the 2 segments are parallel, the perpendicular distances from $\vec{p_0}$ or $\vec{p_1}$ to line $(\vec{a_0}, \vec{a_1})$ are the same. So, the area $A = 0$

When the 2 segments are “co-planar”, the nominator becomes zero too. To check if the two segments overlap or not, simply check x- or y-component is enough.

$\displaystyle A = (a_0.x - a_1.x) (p_0.y - p_1.y) - (a_0.y - a_1.y) (p_0.x - p_1.x)$

$\displaystyle h = \frac{1}{A} \left( (a_0.x - p_1.x) (p_0.y - a_0.y) - (a_0.y - p_1.y) (p_0.x - a_0.x) \right)$

$\displaystyle k = \frac{1}{A} \left( (p_0.x - a_0.x) (a_1.y - p_0.y) - (p_0.y - a_0.y) (a_1.x - p_0.x) \right)$

Solving Lorentz equation and the effect of a radial field in HELIOS

December 22, 2023

GoLuckyRyan Basic HELIOS Leave a comment

Haven’t posted anything for a while, was working on extracting carbon isotopes ESPS, the FSU DAQ, and SOLARIS DAQ for the NUMEN collaboration, went to Hawaii for DNP 2023, and got sick for a week…. Anyway, a friend was having trouble doing a HELIOS simulation with a non-uniform field in GEANT4, and I am curious what is the effect of a non-uniform field, so, I decided to do my own simulation. The core of the simulation is solving the Lorentz equation $\vec{F} = Qe \vec{v} \times \vec{B}$ . This is a 2nd order coupled equation. Well, I already solved it before, so it would be not a huge work.

In this post, I already show that using RK4 method to solve the Lorentz equation. Here, I show the method using the Scipy library in Python and also NDSolve[] in Mathematica.

#Define magnetic field
def Bx(x, y, z):
  return  x/ (1 + math.exp( x**2 + y**2) )
def By(x, y, z):
  return  y/ (1 + math.exp( x**2 + y**2) )
def Bz(x, y, z):
  return -2.5 / (1 + math.exp((math.sqrt(x**2 + y**2) - 0.5)/0.05)) / (1 + math.exp((z**2 - 1**2)/0.1))

from scipy.integrate import odeint

def FieldEquations(w, t, p):
  x , y, z, vx, vy, vz = w

  f = [vx,
       vy,
       vz,
       p * ( vy * Bz(x,y,z) - vz * By(x,y,z)),
       p * ( vz * Bx(x,y,z) - vx * Bz(x,y,z)),
       p * ( vx * By(x,y,z) - vy * Bx(x,y,z))]
  return f

def SolveTrace( helios, thetaCM_deg, phiCM_deg, r0_meter, tMax_ns):
  vx0, vy0, vz0 = helios.velocity_b(thetaCM_deg, phiCM_deg)
  w0 = [r0_meter[0], r0_meter[1], r0_meter[2], vx0, vy0, vz0]

  # ODE solver parameters
  abserr = 1.0e-8
  relerr = 1.0e-6
  stoptime = tMax_ns
  numpoints = 500

  # Create the time samples for the output of the ODE solver.
  t = [stoptime * float(i) / (numpoints - 1) for i in range(numpoints)]

  wsol = odeint(FieldEquations, w0, t, args=(helios.QMb,), atol=abserr, rtol=relerr)

  return [t, wsol]

The code is very simple, we first define the magnetic field components, next, in the FieldEquation, the w is the coordinate vector, and f is the corresponding equation for the coordinate vector, for the x-component, which is

$\displaystyle \frac{dx}{dt} = v_x \\ \frac{dv_x}{dt} = p ( v_y B_z - v_z B_y )$

Similar for the y- and z-components. The variable p is the charge-to-mass ratio. The SolveTrace will use the helios class, which is not shown here, but it basically gives the velocity vector for the initial velocity. The output is the list of the time and the wsol. the wsol is an array of array of x, y, z, vx, vy, and vz.

In Mathematica,

fieldX[x_, y_, z_] := 0.1 If[x == 0, 1, x/Sqrt[x^2 + y^2] Exp[-2 (x^2 + y^2)]]
fieldY[x_, y_, z_] := 0.1 If[y == 0, 1, y/Sqrt[x^2 + y^2] Exp[-2 (x^2 + y^2)]]
fieldZ[x_, y_, z_] := -2.5

G2x[t_, r_, v_] := v[[2]] fieldZ[r[[1]], r[[2]], r[[3]]] - v[[3]] fieldY[r[[1]], r[[2]], r[[3]]];
G2y[t_, r_, v_] := v[[3]] fieldX[r[[1]], r[[2]], r[[3]]] - v[[1]] fieldZ[r[[1]], r[[2]], r[[3]]];
G2z[t_, r_, v_] := v[[1]] fieldY[r[[1]], r[[2]], r[[3]]] - v[[2]] fieldX[r[[1]], r[[2]], r[[3]]]

SolvePath[ChargetoMass_, G2x_Symbol, G2y_Symbol, G2z_Symbol,  r0_, v0_, tMax_] :=
 NDSolve[{
   x''[t] == ChargetoMass G2x[t, { x[t], y[t], z[t]}, { x'[t], y'[t], z'[t]}],
   y''[t] == ChargetoMass G2y[t, { x[t], y[t], z[t]}, { x'[t], y'[t], z'[t]}],
   z''[t] == ChargetoMass G2z[t, { x[t], y[t], z[t]}, { x'[t], y'[t], z'[t]}],
   x[0] == r0[[1]], x'[0] == v0[[1]],
   y[0] == r0[[2]], y'[0] == v0[[2]],
   z[0] == r0[[3]], z'[0] == v0[[3]]
   }, {x[t], y[t], z[t]}, {t, 0, tMax}]

The code defines the field components first, the G2x is the $(\vec{v} \times \vec{B} )_x$ . After that, use the NDSolve to solve for $x(t), y(t), z(t)$ .

A small note on the unit. We are going to solve it in units of meters and nano-seconds. So, the initial position and velocity are in meter and meter/nano-sec respectively. and the charge-to-mass ratio is in SI unit, so, using Coulomb over kg. The charge-to-mass ratio also needs to be multiplied by $10^{-9}$ because the left-hand side is 2nd derivative of time in nano-sec, while the right-hand side needs to be multiplied by $10^{-9}$ to give the correct scale. (think about the original Lorentz equation is in mater and sec)

Now, we are in a position to solve the locus of the charged particle for any field. In particular, I like to find out the effect of a radial component.

I set the radial field to be $B_\rho (\rho, z) = k^2 \exp(-8 \rho^2 )$ , where $k$ is the strength of the radial component. The radial field is strong at small $\rho$ and decay with $\rho$ .

$\displaystyle B_x(x, y, z) = k \frac{x}{\sqrt{x^2+y^2}} \exp(-4 (x^2+y^2)) \\ B_y(x, y, z) = k \frac{y}{\sqrt{x^2+y^2}} \exp(-4 (x^2+y^2)) \\ B_z (x,y, z) = - 2.5$

For the reaction, I pick 15N(d,p) at 10 MeV/u. At $\theta_{CM} = 20 ^\circ$ , The velocity vector is $( 0.0195276, 0, -0.0147025 )$ [meter/ns]. The locus when $k = 0.1$ is

Interestingly, the particle experiences a positive z-acceleration. Imagine the particle only moves in the x-y plane, with a radial component of the field, the particle experiences a positive z-force.

Another effect is the distortion of the x-y projection. Without the radial field, the x-y projection should be a perfect circle. Since the radial field gives a +z-force, and the z-field gives a radially inward force. These two forces combined will twist the path.

Here is the locus for $k = 0.5$

as expected, a stronger radial field, the stronger positive z-force, and the particle change the moving direction more quickly. Also, the x-y project is more distorted.

Split-Pole Kinematics

October 16, 2023

GoLuckyRyan Basic Leave a comment

A spectrograph relates the measured momentum to the excitation energy of the heavy recoil, assuming that the light recoil does not excite. In a two-body transfer reaction, there are 3 degrees of freedom, 2 from the scattering angle, and 1 from the excited energy of the heavy recoil. In a spectrograph, the lab angle of the light recoil is fixed with a small acceptance. Thus, the momentum of the light recoil is correlated with the excitation energy of the heavy recoil.

This post is the normal kinematics of a Split-Pole spectrograph. The given condition is the angle of the Split-Pole, the mass of all 4 particles, and the excitation energy of the heavy recoil. In the classical limit, the conservations of KE and momentum are,

$\displaystyle m_b^2 v_b^2 = m_1^2v_1^2 + m_a^2 v_a^2 - 2 m_1 m_a v_1 v_2\cos\theta_a \\ m_1 v_1^2 + Q c^2= m_a v_a^2 + m_b v_b^2 \\ Q = m_1 + m_2 - m_a - m_b$

The momentum balance is directly from the cosine law.

The job is get the $v_a$ in term of $v_1, \theta_a$ by eliminating $v_b, \theta_b$ . Obviously,

$\displaystyle m_1^2 v_1^2 + m_a^2 v_a^2 - 2 m_1 m_a v_1 v_a \cos\theta_a = m_1 m_b v_1^2 + m_b Q c^2 - m_a m_b v_a^2$

Rearrange,

$\displaystyle m_a(m_a +m_b) v_a^2 + (- 2 m_1 m_a v_1 \cos\theta_a) v_a \\ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + m_1 (m_1 - m_b) v_1^2 - m_b Q c^2 = 0$

Thus, the $v_a$ is

$\displaystyle v_a = \frac{m_1v_1}{m_a+m_b} \left( \cos\theta_a \pm \sqrt{\cos^2\theta_a + H} \right) \\ H = \left(\frac{m_b}{m_a} + 1\right)\left(\frac{m_b}{m_1} - 1 + \frac{m_bQ}{m_1^2}\frac{c^2}{v_1^2}\right)$

When $m_b > m_1$ , the sqrt root is always valid, and $v_a$ always has 2 solutions.

We are interested in the forward center of mass angle, so we take the higher value of $v_a$ .

For a magnetic field $B$ , the radius is

$\displaystyle \rho = \frac{m v_a}{ q e B}$

If the mass is in MeV/c^2, we have

$\displaystyle \rho = \frac{(m c^2 / e) \beta \gamma}{ cqB}, ~~ \beta = \frac{v_a}{c}, \gamma = \frac{1}{\sqrt{1-\beta^2}}$

We can also do it in relativistic way. The conservation of energy and momentum give

$\displaystyle E_1 + m_2 = E_i = \sqrt{m_a^2 + k_a^2} + \sqrt{m_b^2 + k_b^2} \\ k_1^2 + k_a^2 - 2 k_1 k_a \cos\theta_a = k_b^2$

The momentum balance is done by cosine law. Solve for $k_a$ ,

$\displaystyle k_a = \frac{x \pm E_i \sqrt{y}}{z}$

$\displaystyle x = k_1 (E_i^2 - k_1^2 + m_a^2 - m_b^2 ) \cos\theta_a \\ y = (E_i^2 - k_1^2)^2 + (m_a^2-m_b^2)^2 - 2E_i^2 (m_a^2 + m_b^2) + 2k_1^2(\cos(2\theta_a) m_a^2 + m_b^2) \\ z = 2 (E_i^2 - k_1^2 \cos^2\theta_a)$

With the momentum of the light recoil particle $k_a$ in MeV/c, the cyclotron radius is under a magnetic field strength $B$ in T is

$\displaystyle \rho [\textrm{m}] = \frac{k_a [\textrm{MeV/c}] }{ c q B [\textrm{T}]}, ~~ c = 299.792458$

An interesting case is when $k_a = 0$ , i.e. like a head-on collision. In that case, the mass of b is, by solving $x = Ei\sqrt{y}$ ,

$\displaystyle m_b^2 = (E_i \pm m_a)^2 - k_1^2 = (m_1 + m_2 - m_a)^2 + 2T (m_2 - m_a)$

I am not sure what does it means…

The Split-pole spectrograph measures the cyclotron radius for various states of particle b, assuming particle a does not excited. i.e. $\rho(E_x) = \frac{k_a(E_x)}{cqB}$ . Since $c q B$ is a constant, if we can find the formula to convert $k_a$ back to $E_x$ , what would be very nice.

Solve for $m_b$ using the relativistic equation of $k_a$ , or consider the energy and momentum balance

$\displaystyle E_b = E_i - E_a \\ k_b^2 = k_1^2 + k_a^2 - 2 k_1 k_a \cos\theta$

square the $E_b$ , the mass of b is

$\displaystyle m_b^2 = (E_i-E_a)^2 - k_1^2 - k_a^2 + 2 k_1 k_a \cos\theta$

$\displaystyle m_b^2 = E_i^2 - 2 E_i E_a + E_a^2- k_1^2 - k_a^2 + 2 k_1 k_a \cos\theta$

using $E_a^2 = m_a^2 + k_a^2$ to replace $E_a$ , we have

$\displaystyle m_b^2 = E_i^2-k_1^2+m_a^2 + 2 \cos\theta k_1 k_a - 2E_i \sqrt{(k_a^2+m_a^2)}$

Here is a plot for 12C(d,p) at 10 MeV/u at 10 degree. The equation is

$\displaystyle E_x = \sqrt{1.97\times 10^8 + 394 k_a - 28020 \sqrt{m_a^2 + k_a^2}} - m_b$

Compare Non-relativistic and relativistic calculation in HELIOS

October 16, 2023

GoLuckyRyan Basic HELIOS Leave a comment

The classical result for inverse kinematic is, where $m$ is in SI unit,

$\displaystyle z = T_c ( -u \cos\theta + V_{cm} ) \\ E = \frac{1}{2} m(u^2 + V_{cm}^2) - m u V_{cm} \cos\theta$

The relativistic result is where energy is in MeV unit,

$\displaystyle z = \frac{2\pi}{cqB} \left( \gamma \beta E_{cm} - \gamma p \cos\theta\right) \\ E = \gamma E_{cm} - \gamma \beta p \cos\theta$

The classical limit can be reached as

$\displaystyle \gamma \rightarrow 1, ~~ \beta \rightarrow \frac{V_{cm}}{c}$

and the classical cyclotron period is $T_c = 2\pi m / qB$ . We have:

$\displaystyle z = \frac{2\pi}{qB} \left( V_{cm} \frac{E_{cm}}{c^2} - \frac{p}{c} \cos\theta\right) \\ E = E_{cm} - V_{cm} \frac{p}{c} \cos\theta$

We make more approximations,

$\displaystyle E_{cm} \approx mc^2, ~~ p \approx mu c$

Then

$\displaystyle z = \frac{2\pi m}{qB} \left( V_{cm} - u\cos\theta\right) \\ E = E_{cm} - V_{cm} mu \cos\theta$

We kept the $E_{cm}$ in the 2nd equation because we have to subtract $mc^2$ at the end.

$\displaystyle z = T_c \left( V_{cm} - u\cos\theta\right) \\ E = E_{cm} - V_{cm} mu \cos\theta$

We can see that, the z-position assumes the form of the classical case, and by comparison,

$\displaystyle E_{cm} \approx mc^2 + \frac{1}{2} m(u^2 + V_{cm}^2)$

That means, in the classical limit, compared to relativistic, the transformation is simply adding CM frame motion.

The intermediate step in eliminating the $\cos\theta$ from either case is

$\displaystyle E = \frac{1}{2} m(u^2 + V_{cm}^2) - mV_{cm}^2 + \frac{m V_{cm}}{Tc} z\\$

$\displaystyle E = \frac{1}{2} m(u^2 - V_{cm}^2) + \frac{m V_{cm}}{T_c} z\\$

Compared to the final result from a fully relativistic

$\displaystyle E = \frac{1}{\gamma} E_{cm} + \frac{mc\beta}{T_c} z , ~~ T_c = \frac{2\pi}{B} \frac{m}{q}$

It may be misleading that $E_{cm} \approx mc^2 + \frac{1}{2} m(u^2 - V_{cm}^2)$ .

Nuclear Physics 101

Note on lifetime, decay width, and gamma-ray transition probability

Poisson process and exponential decay

Radioactivity and Number of Isotopes

Transform Hamiltonian in rotating frame

GEANT4 installation on MacOS 13.3

Some thoughts on the quenching of spectroscopic factor (II)

Intersect between 2 line segments

Solving Lorentz equation and the effect of a radial field in HELIOS

Split-Pole Kinematics

Compare Non-relativistic and relativistic calculation in HELIOS

Recent Posts

Calendar

Archives

Keywords

# of Visitors

Login