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.