Creating Synthetic Proton Radiographs by Particle Tracing

Proton radiography is a diagnostic technique often used to interrogate the electric and magnetic fields inside high energy density plasmas. The area of interest is positioned between a bright source of protons and a detector plane. Electric and magnetic fields in the plasma deflect the protons, producing patterns on the detector. Since this represents a non-linear and line-integrated measurement of the fields, the interpretation of these "proton radiographs" is complicated.

The SyntheticProtronRadiography class creates a synthetic proton radiographs given a grid of electric and magnetic field (produced either by simulations or analytical models). After the geometry of the problem has been set up, a particle tracing algorithm is run, pushing the protons through the field region. After all of the protons have reached the detector plane, a synthetic proton radiograph is created by making a 2D histogram in that plane.

In [ ]:
import astropy.constants as const
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np

from mpl_toolkits.mplot3d import Axes3D

from plasmapy.diagnostics import proton_radiography as prad
from plasmapy.plasma.grids import CartesianGrid

To illustrate the use of this package, we'll first create an example CartesianGrid object and fill it with the analytical electric field produced by a sphere of Gaussian potential.

In [ ]:
# Create a Cartesian grid
L = 1 *
grid = CartesianGrid(-L, L, num=100)

# Create a spherical potential with a Gaussian radial distribution
radius = np.linalg.norm(grid.grid, axis=3)
arg = (radius / (L / 3)).to(u.dimensionless_unscaled)
potential = 2e5 * np.exp(-(arg ** 2)) * u.V

# Calculate E from the potential
Ex, Ey, Ez = np.gradient(potential, grid.dax0, grid.dax1, grid.dax2)
Ex = -np.where(radius < L / 2, Ex, 0)
Ey = -np.where(radius < L / 2, Ey, 0)
Ez = -np.where(radius < L / 2, Ez, 0)

# Add those quantities to the grid
grid.add_quantities(E_x=Ex, E_y=Ey, E_z=Ez, phi=potential)

# Plot the E-field
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection="3d")
ax.view_init(30, 30)

# skip some points to make the vector plot intelligable
s = tuple([slice(None, None, 6)] * 3)


ax.set_xlabel("X (mm)")
ax.set_ylabel("Y (mm)")
ax.set_zlabel("Z (mm)")
ax.set_xlim(-1, 1)
ax.set_ylim(-1, 1)
ax.set_zlim(-1, 1)
ax.set_title("Gaussian Potential Electric Field")

Prior to running the particle tracing algorithm, the simulation instance must be instantiated by providing some information about the setup, including the locations of the source and detector relative to the origin of the grid.

The source and detector coordinates are entered as a 3-tuple in one of three coordinate systems: Cartesian ($x$, $y$, $z$), spherical ($r$, $\theta$, $\phi$) or cylindrical ($r$, $\theta$, $z$). All values should be astropy.units.Quantity instances with units of either length or angle. The vector from the source to the detector should pass through the origin to maximize the number of particles that pass through the simulated fields.

In [ ]:
source = (0 *, -10 *, 0 *
detector = (0 *, 100 *, 0 *

sim = prad.SyntheticProtonRadiograph(grid, source, detector, verbose=True)

Note that, since the example grid did not include a B-field, the B-field is assumed to be zero and a warning is printed.

Next, a distribution of nparticles simulated particles of energy particle_energy is created using the create_particles() function. Setting the max_theta parameter eliminates particles with large angles (relative to the source-detector axis) which otherwise would likely not hit the detector. Particles with angles less than $\theta_{max}$ but greater than $\theta_{track}$ in the setup figure above will not cross the grid. These particles are retained, but are coasted directly to the detector plane instead of being pushed through the grid.

The particle keyword sets the type of the particles being traced. The default particle is protons, which is set here explicitly to demonstrate the use of the keyword.

By default, the particle velocities are initialized with random angles (a Monte-Carlo approach) with a uniform flux per unit solid angle. However, particles can also be initialized in other ways by setting the distribution keyword.

In [ ]:
sim.create_particles(1e5, 3 * u.MeV, max_theta=np.pi / 15 * u.rad, particle="p")

The simulation is now ready to run. In brief, the steps of the simulation cycle are as follows:

  1. Particles that will never hit the field grid are ignored (until a later step, when they will be automatically advanced to the detector plane).
  1. Particles are advanced to the time when the first particle enters the simulation volume. This is done in one step to save computation time.
  1. While particles are on the grid, the particle pusher advances them each timestep by executing the following steps:

    A. The fields at each particle's location are interpolated using the interpolators defined in the AbstractGrid subclasses.

    B. The simulation timestep is automatically (and adaptively) calculated based on the proton energy, grid resolution, and field amplitudes. This timestep can be clamped or overridden by setting the dt keyword in the run() function.

    C. An implementation of the Boris particle push algorithm is used to advance the velocities and positions of the particles in the interpolated fields.

  1. After all of the particles have left the grid, all particles are advanced to the detector plane (again saving time). Particles that are headed away from the detector plane at this point are deleted, as those particles will never be detected.

When the simulation runs, a progress meter will show the number of particles currently on the grid. This bar will start at zero, increase as particles enter the grid, then decrease as they leave it. When almost all particles have left the grid, the simulation ends.

In [ ]:

The following plot illustrates that, after the simulation has ended, all particles have been advanced to the detector plane.

In [ ]:
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection="3d")
ax.view_init(30, 150)
ax.set_xlabel("X (cm)")
ax.set_ylabel("Y (cm)")
ax.set_zlabel("Z (cm)")

# Plot the source-to-detector axis
    sim.source[0] * 100,
    sim.source[1] * 100,
    sim.source[2] * 100,
    sim.detector[0] * 100,
    sim.detector[1] * 100,
    sim.detector[2] * 100,

# Plot the simulation field grid volume
ax.scatter(0, 0, 0, color="green", marker="s", linewidth=5, label="Simulated Fields")

# Plot the the proton source and detector plane locations
    sim.source[0] * 100,
    sim.source[1] * 100,
    sim.source[2] * 100,

    sim.detector[0] * 100,
    sim.detector[1] * 100,
    sim.detector[2] * 100,

# Plot the final proton positions of some (not all) of the protons
ind = slice(None, None, 200)
    sim.x[ind, 0] * 100, sim.x[ind, 1] * 100, sim.x[ind, 2] * 100, label="Protons",


A 'synthetic proton radiograph' can now be constructed by creating a 2D histogram of proton positions in the image plane. The synthetic radiograph function takes two keywords:

  • 'size' gives the locations of the lower left and upper right corners of the detector grid in image plane coordinates.

  • 'bins' is the number of histogram bins to be used in the horizontal and vertical directions. Using more bins creates a higher resolution image, but at the cost of more noise.

In [ ]:
# A function to reduce repetative plotting
def plot_radiograph(hax, vax, intensity):
    fig, ax = plt.subplots(figsize=(8, 8))
    plot = ax.pcolormesh(,,
    cb = fig.colorbar(plot)"Intensity")
    ax.set_xlabel("X (cm), Image plane")
    ax.set_ylabel("Z (cm), Image plane")
    ax.set_title("Synthetic Proton Radiograph")

size = np.array([[-1, 1], [-1, 1]]) * 1.5 *
bins = [200, 200]
hax, vax, intensity = sim.synthetic_radiograph(size=size, bins=bins)
plot_radiograph(hax, vax, intensity)

As expected, the outward-pointing electric field in the sphere has deflected the protons out of the central region, leaving a dark shadow.

Kugland et al. 2012 and Bott et al. 2017 define the dimensionless "contrast parameter" that separates different regimes of proton radiography:

\begin{equation} \mu = \frac{l \alpha}{a} \end{equation}

Where $l$ is the distance from the source to the grid, $a$ is the spatial scale of the scattering electromagnetic fields, and $\alpha$ is the particle deflection angle. The value of $\mu$ can fall in one of three regimes:

\begin{align} \mu &\ll 1 \rightarrow \text{ linear}\\ \mu &< \mu_c \rightarrow \text{ nonlinear injective}\\ \mu &> \mu_c \rightarrow \text{ caustic}\\ \end{align}

where $\mu_c \sim 1$ is a characteristic value at which particle paths cross, leading to the formation of bright caustics. Correctly placing a radiograph in the correct regime is necessary to determine which analysis techniques can be applied to it.

The maximum deflection angle can be calculated after the simulation has run by comparing the initial and final velocity vectors of each particle

In [ ]:
max_deflection = sim.max_deflection
print(f"Maximum deflection α = {np.rad2deg(max_deflection):.2f}")

The spatial scale of the field constructed in this example is $\sim$ 1 mm, and $l$ is approximately the distance from the source to the grid origin. Therefore, we can calculate the value of $\mu$

In [ ]:
a = 1 *
l = np.linalg.norm(sim.source * u.m).to(
mu = l * max_deflection.value / a
print(f"a = {a}")
print(f"l = {l:.1f}")
print(f"μ = {mu:.2f}")

which places this example in the non-linear injective regime.


For sake of comparison, here is the result achieved by setting distribution = 'uniform' in the create_particles() function.

In [ ]:
    1e5, 3 * u.MeV, max_theta=np.pi / 15 * u.rad, distribution="uniform"
size = np.array([[-1, 1], [-1, 1]]) * 1.5 *
bins = [200, 200]
hax, vax, intensity = sim.synthetic_radiograph(size=size, bins=bins)
plot_radiograph(hax, vax, intensity)