Simulating Particle Trajectories

This article discusses numerical simulations of particle trajectories in the "silver" permanent magnet.

Description of Problem

The general problem of modeling relativistic particle trajectories in electromagnetic fields boils down to solving the following system of six differential equations:

[math]\frac{d\vec{x}}{dt}=\vec{v}=\frac{1}{\gamma m_0}\vec{p}[/math]

[math]\frac{d\vec{p}}{dt}=q\left(\vec{E}+\frac{1}{\gamma m_0}\vec{p}\times\vec{B}\right)[/math]

where [math]\vec{x}[/math] is the position of the particle, [math]\vec{p}[/math] is the momentum of the particle, [math]m_0[/math] is the rest energy of the particle, and [math]\gamma[/math] is the usual relativistic factor. These equations are solved subject to an appropriate set of boundary conditions, usually the initial position and momentum of the particle.

In simple cases, such as when the fields are uniform, this set of equations can be solved analytically. However, for the "silver" permanent magnet we plan on using for the pair spectrometer, the fields are anything but uniform. In such a case, the system must be solved numerically.

Description of Code

The code, found here, solves the above system of equations using an embedded 8th order Runge-Kutta Prince-Dormand method with 9th order error estimates with adaptive stepsize control, as implemented in the GNU Scientific Library [1]. The zip file contains the relevant source and a script which should build the executable on any unix system with GSL installed in the usual place.

When executing the code, you must specify an input file on the command line. The input file is fairly simple. The first line of the input file is either a "0" or some integer greater than zero, and tells the program to run in "detailed output" or "scan" mode, respectively. In "detailed output" mode, the next five lines specify the energy of the initial gamma ray, the energy of the electron (the energy of the positron is assumed to be the difference of these two numbers), the angle between the axis of the magnet and the photon beam, the "Y" position of the conversion site above the magnet yoke, and the "X" position of the conversion site relative to the left edge of the magnet. In "scan" mode, the integer on the first line tells the program how many cases to run. Only two additional lines are needed: the energy of the photon and the "X" position of the conversion site. Examples of such input files are provided here (detailed mode) and here (scan mode).

The program output also depends on the program mode. In the "detailed output" mode, everything is written to stdout, and can be directed to a file using the ">" director. The output begins with a header, each line of which begins with the "#" character. The header lists the physical properties of the particle, the conditions read from the input file, and the corresponding initial positions and momenta of the electron and positron. Following the header, the program writes four columns of data: the time (in seconds) and the x, y, and z coordinates (in meters).

In "scan" mode, the program randomly selects an electron energy, angle of the magnet axis relative to the beam, and a conversion site "Y" coordinate, all subject to appropriate reasonable limits and kinematic constraints. It then models the trajectory of the particle until either the particle is sufficiently far from the detector, or a certain maximum time is reached. The program then evaluates the angles of the particles' final trajectories relative to the initial photon beam and assigns a flag as follows:

Limitations of the Code

Results of parameter space scans

"Ideal" orientation of magnet for 5 MeV photons

Pairs produced elsewhere on the beam