High-dimensional phase space measurements
1. Predicting halo formation
Linear accelerators (linacs) generate powerful pulsed hadron beams for various applications. Increasing the beam power/intensity by orders of magnitude would benefit particle physics, nuclear physics, condensed matter physics, material science, biology, nuclear engineering, and other fields. Unfortunately, significantly increasing the beam intensity is not possible due to beam loss which occurs when particles hit the walls of the accelerator during their journey. The total lost beam power must be kept below 1 watt (W) per meter to avoid dangerous radiation levels. For a 1 MW beam, this limit corresponds to a fractional loss of \(10^{-6}\).
Simulations struggle to predict the number and location of lost particles. The dynamics in the accelerator are complex, but the situation is perhaps unsurprising for other reasons:
Accelerators are composed of thousands of components distributed over hundreds of meters. Uncertainties in each component’s position, alignment, amplitude, etc., lead to errors in the applied electromagnetic field at each position.
The initial distribution of particles in position-momentum space (phase space) is usually unknown.
These issues make it difficult to predict the low-order moments of the beam distribution, let alone low levels of beam loss. The second issue is critical because of the violent space charge forces that influence the beam evolution, coupling the output to the input phase space distribution [1,2].1
Take the Spallation Neutron Source (SNS) as an example. The SNS accelerates negative hydrogen (H-) ions through a 500-meter linac. Hundreds of beam loss monitors (BLMs) record the total number of lost particles as a function of position. Phase space distribution measurements are limited to low-dimensional projections in the low-energy stages (near the boxed region in Fig. 1). Simulations based on these measurements and the accelerator’s design parameters cannot predict the measured beam loss.
The LEDA (Low Energy Test Accelerator) experiment at Los Alamos National Laboratory attempted to predict halo formation over a much shorter distance [3]. One-dimensional (1D) density profiles of a low-energy proton beam were measured after a series of focusing quadrupole magnets. Simulations failed to reproduce the low-density “halo” surrounding the beam core — see Fig. 2. The authors attributed the failure to an inaccurate phase space distribution used to seed the simulation. (Only low-order moments of the distribution were measured.)
2. High-dimensional phase space measurements
We are continuing experimental work on halo prediction at the SNS Beam Test Facility (BTF), a recently commissioned replica of the front end of the SNS linac [4]. A diagram of the BTF is below.
The beam starts as a slow-moving, continuous stream of H- ions extracted from the plasma source. The low-energy beam transport (LEBT) uses electrostatic focusing to guide the beam to the radiofrequency quadrupole RFQ. The RFQ bunches and accelerates the beam to 2.5 MeV. Each bunch drifts through the medium-energy beam transport (MEBT) and finally through a series of permanent quadrupole magnets arranged in a focusing-off-defocusing-off (FODO) pattern. (If halo formation occured, it would likely be here due to the rapid oscillation of the beam core.)
The apparatus just after the RFQ (blue region) measures the particle density in six-dimensional phase space — three positions (\(x\), \(y\), \(z\)) and three momenta (\(p_x\), \(p_y\), \(p_z\)), with the \(z\) axis parallel to the beam. It works by measuring the charge within a small region \(\mathbf{x} \pm \Delta\), where \(\mathbf{x} = [x, p_x, y, p_y, z, p_z]^T\), scanning \(\mathbf{x}\) through the space, and interpolating the data to obtain an estimate of the distribution function \(f(\mathbf{x})\).
The region \(\mathbf{x} \pm \Delta\) is selected using slits. First, a horizontal-vertical pair of slits selects a square in the \(x\)-\(y\) plane. Next, a second pair of slits select \(p_x\) and \(p_y\) (based on their positions relative to the first slits). The momentum \(p_z\) is selected by a vertical slit after a dipole magnet, which bends on-energy particles 90 degrees (faster particles bend less). The last to go is \(z\): The beam passes through a wire, generating an electron beam with the same temporal structure. The electrons are streaked in the transverse plane using an oscillating electric field such that the image of the beam on a screen gives the \(z\) distribution.
We proceed by scanning the slits in a nested loop and recording the image of the beam on the screen at each setting. Each pixel in the image is the intensity at a point in phase space. This measurement is the first of its kind [5].
Once we reconstruct the initial distribution \(f(\mathbf{x})\), we generate a set of phase space coordinates \(\{\mathbf{x}_1, \dots \mathbf{x}_n\}\) by sampling from the distribution.2 These coordinates are propagated through a computational accelerator model and compared with measurements at the end of the beamline. Here we sacrifice dimensionality for dynamic range to image the beam halo. Instead of measuring the one-dimensional projections \(\{ f(x), f(y) \}\), we measure two-dimensional projections \(\{ f(x, p_x), f(y, p_y) \}\). We can reach a dynamic range of \(10^6\) (the dynamic range in Fig. 2 is between \(10^3\) and \(10^4\)). This high-dynamic-range measurement is the first of its kind [6].
When combined with a detailed model of the accelerator and a well-tested simulation code, we hope these measurements will allow us to predict the halo dynamics in the BTF. Such a prediction would be the first step toward model-based loss prediction in a linear accelerator.
Our 6D measurements are currently quite slow and fail to capture sharp features in the distribution. Dropping one dimension (\(z\)) increases the resolution and dynamic range while capturing most of the significant correlations in the data. The resulting 5D measurement is useful for in-depth examinations of the phase space distribution. I’ve gotten a lot of mileage out of a high-resolution 5D measurement I took in June 2022, which is publically available [7]. The scan generating over 285,082 images over 18 hours. Fig. 6. shows the change in image brightness as the slits move in and out of the beam core.
We cropped and downscaled the images to decrease the size of the data set from around 250 gigabytes to 20 gigabytes. Next, we interpolated the pixel intensities to obtain a \(69 \times 88 \times 69 \times 65 \times 55\) image of the phase space distribution \(f(x, x', y, y', w)\). Here we use the energy \(w\) instead of the momentum \(p_z\). (Also \(x'\) and \(y'\) are basically the same as \(p_x\) and \(p_y\).)
3. High-dimensional data visualization
We have spent significant time analyzing the measured initial phase space distribution [5,8–10]. This task is straightforward in two dimensions: one simply looks at an image representing \(f(x, p_x)\) and extracts features from it. High-dimensional data is exponentially more challenging to visualize and interpret.
Here it is helpful to distinguish projections, slices, and partial projections. A projection is an integral/sum that squashes the distribution onto a lower-dimensional space. For example: \[ f(x, p_x) = \iiiint_{-\infty}^{+\infty} {f(x, p_x, y, p_y, z, p_z) dy dp_y dz dp_z}. \] We must eventually project the distribution onto a two-dimensional space to view it on a page or screen. For an \(n\)-dimensional distribution, there are \(n\) one-dimensional projections and \(n (n - 1) / 2\) two-dimensional projections. The classic corner plot displays all these projections in a single figure.
Corner plots are valuable tools, but it is important to remember that they hide relationships between three or more dimensions. Higher dimensional correlations are only visible in slices. Most people are familiar with slices (i.e., cross sections) from three-dimensional medical CT scans. CT scans map the three-dimensional density of, say, a human brain. To examine the brain’s internal structure, medical professionals scroll through a set of two-dimensional images. Each image is the density on a two-dimensional plane intersecting the three-dimensional distribution. We call this a “planar slice”.
Add one more dimension. How is a slice defined now? Well, we could look at the intersection of a three-dimensional plane with the four-dimensional distribution. This is equivalent to fixing the value of one of the coordinates, leaving a conditional distribution like \(f(x, p_x, p_y \vert y{=}0)\). But there is now more than one dimension that we could slice. We could take another three-dimensional plane, orthogonal to the first plane, and view the intersection of both planes and the four-dimensional distribution. We are left with a two-dimensional surface like \(f(x, p_x \vert y{=}p_y{=}0)\).
Since slices can be projected and projections can be sliced, we must distinguish between a full projection and a partial projection — the projection of a slice. More complex slicing/projection sequences are obviously possible in five- or six-dimensional spaces.
The “slice matrix” plot in Fig. 8 is one way to visualize the effect of slicing a four-dimensional distribution. A two-dimensional projection is on the bottom right. Three-dimensional projections — obtained by slicing only one dimension — are on the bottom and right panels. The four-dimensional distribution — obtained by slicing two dimensions — is shown on the top left.
Such figures render slowly and are limited to four-dimensional data. Complete slicing flexibility is available in an interactive widget shown in Fig. 9. This function takes one or more images or point clouds of any dimensionality as inputs. I use it extensively in my data analysis. It is available in the psdist package.
There are other ways to slice and project the data. For example, we could extend the concept of a radial histogram. A radial histogram takes a distribution \(f(x_1, \dots, x_n)\) to a one-dimensional profile \(f(r_n)\), where \(r_n = \sqrt{x_1^2 + \dots + x_n^2}\), by integrating the density on spheres of radius \(r_n\).3 Instead of defining the spheres in the \(n\)-dimensional space, suppose we defined them in a subspace like the \(x_2\)-\(\dots\)-\(x_n\) plane and observed the \(x_1\) distribution on each sphere. In other words, compute \(f(x_1, r(x_2, \dots, x_n))\). The resulting projection is more difficult to interpret but preserves high-dimensional correlations in the data.
We could also define shells as the density contours of the distribution in the subpace — \(f(x_2, \dots, x_n)\) in the above example. We could label these “contour shell slices”. Fig. 10 illustrates this concept. I applied this idea to our measured distribution to view a specific feature in a compact figure. In Fig. 11, I defined density contours in the four-dimensional transverse phase space by thresholding the projection \(f(x, x', y, y')\) between two density levels. I then computed the energy distribution within each shell (possible thanks to numpy
).
Notice the drastic change in energy distribution — unimodal to bimodal — as we move from the outer to the inner core in the transverse phase space. This energy hollowing represents a five-dimensional correlation; it is also visible (although perhaps harder to see) in Fig. 8. I also plotted the transverse projections of the slice farthest from the core. Pretty bizarre-looking images result, but they make some sense; they select particles with large transverse amplitude. Some intuition is gained by considering the projections of hollow three-dimensional distribution.
4. Ongoing work
We would like to explain the origin of all features we find in the distribution. From direct measurements at lower beam currents, we know the space charge drives the energy hollowing in Fig. 10.
Furthermore, we are confident that the energy hollowing develops in the RFQ rather than the MEBT. I am currently using simulations to study the dynamics that drive the hollowing in the RFQ.
Another (more practical) question is: how do “hidden” features in the initial distribution affect the downstream dynamics? In particular, how do they affect halo formation in the BTF or the SNS linac? I am also studying this problem using simulations.
Improvements to the high-dimensional measurements and the accelerator model will resume when the BTF comes back online this fall.
References
Footnotes
Space charge forces arise from the electric field generated by the beam. The strength of space charge forces scales inversely with the beam energy.↩︎
Sampling from a high-dimensional distribution is not trivial. We currently use a simple method that treats the distribution as a histogram. I am looking into approaches that will scale to high-dimensional distribution functions, such as normalizing flows.↩︎
One could just as well define ellipsoids instead of spheres by computing the covariance matrix \(\mathbf{\Sigma} = \mathbf{x}\mathbf{x}^T\) and setting \(r = \mathbf{x}^T \mathbf{\Sigma}^{-1} \mathbf{x}\).↩︎