Modeling Magnetic Fields with Helical Solutions to Laplace's Equation

TL;DR

This paper develops a series solution to Laplace's equation in helical coordinates, enabling precise modeling of solenoidal magnetic fields with minimal data and high accuracy, capturing complex helical features.

Contribution

It introduces a novel series solution in helical coordinates and demonstrates its effectiveness for accurate magnetic field modeling with sparse data.

Findings

Achieves better than one part per million accuracy in magnetic field modeling.

Effectively captures helical features from solenoid winding.

Uses minimal parameters and data for high-precision modeling.

Abstract

The series solution to Laplace's equation in a helical coordinate system is derived and refined using symmetry and chirality arguments. These functions and their more commonplace counterparts are used to model solenoidal magnetic fields via linear, multidimensional curve-fitting. A judicious choice of functional forms, a small number of free parameters and sparse input data can lead to highly accurate, fine-grained modeling of solenoidal magnetic fields, including helical features arising from the winding of the solenoid, with overall field accuracy at better than one part per million.

Figures

Tables

Table 1. Table 1 : Table of simulation parameters for different solenoid examples.

Parameter	Solenoid A	Solenoid B	Solenoid C
Length	92 m	9.2 m	9.2 m
Radius	250 mm	1000 mm	250 mm
Pitch	10 cm	7.5 mm	10 cm
Current	6114 A	6114 A	6114 A
Average $B_z$	767 G	1 T	767 G

Equations

$${\nabla }^2\Phi =0.$$

$$\Phi (x,y,z)=\sum _n{c}_n{\psi }_n(x,y,z).$$

$${\chi }^2=\sum _i{\Vert -\nabla \Phi (x_i,y_i,z_i)-\mathbf{B}(x_i,y_i,z_i)\Vert }^2.$$

$$\frac{{\partial }^2\Phi }{\partial r^2}+\frac{1}{r}\frac{\partial \Phi }{\partial r}+\frac{1}{r^2}\frac{{\partial }^2\Phi }{\partial {\theta }^2}+\frac{{\partial }^2\Phi }{\partial z^2}=0.$$

$$\frac{R^{\prime \prime }}{R}+\frac{1}{r}\frac{R^{\prime }}{R}+\frac{1}{r^2}\frac{T^{\prime \prime }}{T}+\frac{Z^{\prime \prime }}{Z}=0.$$

$$\frac{R^{\prime \prime }}{R}+\frac{1}{r}\frac{R^{\prime }}{R}-\frac{n^2}{r^2}+c_z=0.$$

$$\begin{array}{rl}{\Phi }_{n,k}(r,\theta ,z)& =J_n(kr)\cos (n\theta )\cosh (kz)\\ {\Phi }_{n,k}(r,\theta ,z)& =I_n(kr)\cos (n\theta )\cos (kz).\end{array}$$

$$\begin{array}{rl}\mathbf{h}(t)& =(x(t),y(t),z(t))=(R\cos (t),R\sin (t),\overline{p}t)\\ & =(r(t),\theta (t),z(t))=(R,t,\overline{p}t).\end{array}$$

$$x$$

$$y$$

$$z$$

$$\frac{\partial (r,\theta ,z)}{\partial (r,\varphi ,\zeta )}=\left[\begin{array}{ccc}{r}_r& r_{\varphi }& r_{\zeta }\\ {\theta }_r& {\theta }_{\varphi }& {\theta }_{\zeta }\\ z_r& z_{\varphi }& z_{\zeta }\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& \overline{p}& 1\end{array}\right].$$

$$(r(t),\varphi (t),\zeta (t))=(R,t,0).$$

$$(r,\varphi +2\pi k,\zeta )=(r,\varphi ,\zeta +kP).$$

$$\frac{\partial f}{\partial \varphi }=\frac{\partial f}{\partial r}\frac{\partial r}{\partial \varphi }+\frac{\partial f}{\partial \theta }\frac{\partial \theta }{\partial \varphi }+\frac{\partial f}{\partial z}\frac{\partial z}{\partial \varphi },$$

$$\frac{\partial }{\partial r}=\frac{\partial }{\partial r}$$

$$\frac{\partial }{\partial \varphi }=\frac{\partial }{\partial \theta }+\overline{p}\frac{\partial }{\partial z}$$

$$\frac{\partial }{\partial \zeta }=\frac{\partial }{\partial z}.$$

$$\dfrac{\partial^{2}u}{\partial{r}^{2}}+\frac{1}{r}\dfrac{\partial u}{\partial r}+\frac{1}{r^{2}}\bigg{(}\dfrac{\partial^{2}u}{\partial{\phi}^{2}}-2\overline{p}\dfrac{\partial^{2}u}{\partial{\phi}\partial{\zeta}}+\overline{p}^{2}\dfrac{\partial^{2}u}{\partial{\zeta}^{2}}\bigg{)}=0.$$

$$r^2{R}^{\prime \prime }+r{R}^{\prime }+r^2{\omega }_2^{2}R+({\omega }_1-\overline{p}{\omega }_2{)}^{2}R=0.$$

$$r^2{R}^{\prime \prime }+r{R}^{\prime }+r^2{\omega }_2^{2}R+({\omega }_1+\overline{p}{\omega }_2{)}^{2}R=0.$$

$${\Phi }_{{\omega }_1,{\omega }_2}(r,\varphi ,\zeta )$$

$${\Phi }_{n,k}(r,\varphi ,\zeta )$$

$${\Phi }_{n,k}(r,\theta ,z)$$

$${\Phi }_{n,k}(r,\varphi ,\zeta )$$

$${\Phi }_{n,k}(r,\theta ,z)$$

$$R(r_0)\left(A_{n,m}\cos (n\theta -k_{m}z)+B_{n,m}\sin (n\theta -k_{m}z)\right).$$

$$\mathbf{B}(\mathbf{x})=\frac{{\mu }_0}{4\pi }\int \frac{Id\mathbf{h}\times (\mathbf{x}-\mathbf{h})}{|\mathbf{x}-\mathbf{h}{|}^3}.$$

$$\begin{array}{rl}\Phi (r,\theta ,z)=& \sum _{m,n}{I}_n(k_{m}r)[C_n\sin (n\theta )+(1-C_n)\cos (n\theta )]\cdot \\ & [A_{m,n}\cos (k_{m}z)+B_{m,n}\sin (k_{m}z)].\end{array}$$

$$\begin{array}{rl}{B}_r=& \sum _{m,n}{I}_n^{\prime }(k_{m}r)[C_n\sin (n\theta )+(1-C_n)\cos (n\theta )]\cdot \\ & [A_{m,n}\cos (k_{m}z)+B_{m,n}\sin (k_{m}z)]k_m\\ B_{\theta }=& \sum _{m,n}{I}_n(k_{m}r)[C_n\cos (n\theta )-(1-C_n)\sin (n\theta )]\cdot \\ & [A_{m,n}\cos (k_{m}z)+B_{m,n}\sin (k_{m}z)]\frac{n}{r}\\ B_z=& \sum _{m,n}{I}_n(k_{m}r)[C_n\sin (n\theta )+(1-C_n)\cos (n\theta )]\cdot \\ & [-A_{m,n}\sin (k_{m}z)+B_{m,n}\cos (k_{m}z)]k_m,\end{array}$$

$$\begin{array}{rl}\Phi (r,\theta ,z)=\sum _{m,n}{I}_n\left(h_{m}r\right)[& D_{m,n}\cos (n\theta -h_{m}z)+\\ & E_{m,n}\sin (n\theta -h_{m}z)],\end{array}$$

$$\begin{array}{rl}{B}_r=& \sum _{m,n}{I}_n^{\prime }\left(h_{m}r\right)[D_{m,n}\cos (n\theta -h_{m}z)+\\ & E_{m,n}\sin (n\theta -h_{m}z)]h_m\\ B_{\theta }=& \sum _{m,n}{I}_n\left(h_{m}r\right)[-D_{m,n}\sin (n\theta -h_{m}z)+\\ & E_{m,n}\cos (n\theta -h_{m}z)]\frac{n}{r}\\ B_z=& \sum _{m,n}{I}_n\left(h_{m}r\right)[D_{m,n}\sin (n\theta -h_{m}z)-\\ & E_{m,n}\cos (n\theta -h_{m}z)]h_m.\end{array}$$

Full text

High Energy Physics; Magnetic Fields; Numerical Methods

I Introduction

The use of superconducting solenoids is commonplace in scientific research, medical, and industrial applications. The ever-increasing demand for higher resolution (e.g. medical imaging, charged particle reconstruction) requires increased accuracy and precision for the magnetic field Feher et al. (2018). The physical geometry and helical winding of a solenoid generates an axial asymmetry that may cause significant deviations of the magnetic field from that of an ideal solenoid Nouri and Plaster (2013); Gu et al. (2011), and it may be necessary to account for these deviations. The deviations will be especially pertinent for solenoidal coils that have a relatively large pitch; some large-pitch solenoids are currently employed in active particle physics experiments Wittgenstein et al. (1990); Swoboda et al. (2002). A common analysis approach used within the high energy physics community is to combine a set of sparse field measurements (e.g. provided by 2D or 3D Hall or NMR probes) with a series of functions derived from solutions to Maxwell’s equations to reconstruct the magnetic field within a given volume Amapane et al. (2005); Klyukhin et al. (2010); Aleksa et al. (2008). These methods have not attempted to use functions derived with helical symmetry, however, which may be a better choice for modeling real solenoids.

Any magnetic field must obey Maxwell’s equations. For a static field in a region free of current and magnetic materials, the magnetic field ${\bf B}$ can be expressed as ${\bf B}=-\nabla\Phi$, where the scalar field $\Phi$ satisfies Laplace’s equation:

[TABLE]

When treated as a boundary-value problem, Eq. 1 can sometimes be solved via a separation of variables, and in the case of solenoids, this is typically done in a cylindrical coordinate system Jackson (1999). The axial symmetry inherent to the cylindrical coordinate system is broken by the helical winding of the solenoid, however. In this case it is more appropriate to use a helical coordinate system, for which solutions to Laplace’s equation via separation of variables are known to exist Overfelt (2001). This coordinate system naturally encompasses the axial asymmetry of a helically wound solenoid, and thus is a better choice for modeling the associated magnetic field.

In this paper, we first derive the series solution to Laplace’s equation in a helical coordinate system. We then calculate magnetic fields due to solenoids with different radius, pitch, and length properties, and from those calculations produce datasets that are used as inputs for a parametric modeling of the field. Using a judicious choice of terms for the helical series, cylindrical series, and other physically valid functions for magnetic fields, we construct a model of the calculated magnetic fields using only a sparse subset of the data and a linear sum of analytical functions whose coefficients are determined through least squares fitting. We then compare the reconstructed field with a high-granularity dataset from the same calculated field in order to determine the overall precision and accuracy of the model.

This paper is structured as follows. In Section II, we introduce the mathematical formalism and solve Laplace’s equation in a helical coordinate system. In Section III, we establish the exact definitions of the functional forms used for modeling magnetic fields in the following section. In Section IV, we describe the process used to calculate helically-wound solenoids, and then model three specific examples using the functions discussed in the previous sections. For the third example, we investigate the impact of measurement error and data density on the model performance. We discuss the results and implications of these modeling efforts in Section V.

II Mathematical Formalism

Representing a scalar potential function as an infinite linear combination of orthogonal basis functions $\{\psi_{n}\}$ is a standard technique in potential theory, especially in service of solving boundary value problems for Laplace’s equation.

[TABLE]

The convergence of this series is understood in the sense that $\int|\Phi-\Phi_{N}|^{2}dV\rightarrow 0$, as $N\rightarrow\infty$, where $\Phi_{N}=\sum_{n=1}^{N}c_{n}\psi_{n}(x,y,z)$ is the $N^{th}$ partial sum in the series representation of $\Phi$ Evans (2010). For special geometries, various classical families of basis functions exist. For example, we will recall in what follows the various forms of the harmonic functions which separate in cylindrical coordinates, and lead to complete orthogonal bases of harmonic functions for certain function spaces Jackson (1999). For our application, it is not useful to have representations of potentials in an infinite expansion unless we can well represent the target field using only a small number of terms. We are motivated to choose the family of basis functions so that the above convergence happens as fast as possible.

Given field data ${\bf B}(x_{i},y_{i},z_{i})$ for a finite set of discrete points $(x_{i},y_{i},z_{i})$, we propose a method to find a scalar potential function $\Phi(x,y,z)$ which minimizes the sum of squared residuals:

[TABLE]

We will use a series representation $\Phi=\sum_{j}a_{j}\psi_{j}^{(cyl)}+\sum_{j}b_{j}\psi_{j}^{(hel)}$ in terms of a family of harmonic functions which we motivate below by the geometry of a solenoid represented as a truncated helical curve. The cylindrical symmetry of a helix will cause us to work in cylindrical coordinates $(r,\theta,z)$ at first, when we recall some facts about separable cylindrical harmonics. Later we will work in a helical coordinate system introduced by Waldron Waldron (1958) and introduce a family of separable harmonic functions suited to the geometry of any particular semi-idealized solenoid. Using this functional form for the scalar potential $\Phi$ will allow us to model features of the magnetic field that are due to the helical nature of the coil and the finite length of the solenoid.

II.1 Cylindrical Harmonics

We will briefly discuss the standard procedure of separation of variables of Laplace’s equation in cylindrical coordinates, in order to illuminate the connections to the following analogous work in the helical coordinate system.

In cylindrical coordinates $(r,\theta,z)$, given by the standard change of variables to Cartesian coordinates $x=r\cos(\theta),y=r\sin(\theta),z=z$, Laplace’s equation takes the following form:

[TABLE]

If we assume that $\Phi$ separates in the cylindrical variables, $\Phi(r,\theta,z)=R(r)T(\theta)Z(z)$, we see that the functions $R(r),T(\theta),Z(z)$ must satisfy the following:

[TABLE]

The quantities ${Z^{\prime\prime}}/{Z}$ and ${T^{\prime\prime}}/{T}$ are constant functions of the variables, so that $Z^{\prime\prime}=c_{z}Z$, $T^{\prime\prime}=c_{\theta}T$ for some real constants $c_{z},c_{\theta}$. The function $T(\theta)$ is $2\pi$-periodic, which forces $c_{\theta}=-n^{2}$ for some integer $n$. Consequently, $R$ must satisfy the following ordinary differential equation, which can be reduced to a Bessel equation of order $n$ using the linear change of variables $s=cr$, where $c=\sqrt{c_{z}}$:

[TABLE]

These relationships lead to the following forms for separable cylindrical harmonics depending on the sign of $c_{z}$. If $c_{z}=k^{2}$, then $Z(z)=A\cosh(kz)+B\sinh(kz)$ and $R(r)=CJ_{n}(kr)+DY_{n}(kr)$, where $J_{n}$ and $Y_{n}$ are the order $n$ Bessel functions of the first and second kind. If $c_{z}=-k^{2}$, then $Z(z)=A\cos(kz)+B\sin(kz)$ and $R(r)=CI_{n}(kr)+DK_{n}(kr)$, where $I_{n}$ and $K_{n}$ are the order $n$ modified Bessel functions of the first and second kind.

Because we are interested in potentials which do not have singularities at $r=0$, we discard the Bessel functions of the second kind. Thus, for each choice of a positive integer $n$ and positive real $k$, we have harmonic functions of the following form (note, for each term, there are three other terms corresponding to replacing $\cos$ and $\cosh$ with $\sin$ and $\sinh$, respectively):

[TABLE]

This family of cylindrical harmonics is indexed by the integer $n$ and the real number $k$. When interested in solving a boundary value problem on a cylinder $r<R$, $|z|<L$, one need only consider the countable complete orthogonal basis corresponding to choosing $k$ to be $k_{m}=m\pi/L$, where $m$ is an integer.

An alternative approach is to assume that $T(\theta)=\exp(\omega_{1}\theta)$ and $Z(z)=\exp(\omega_{2}z)$ are complex exponential functions parametrized by the complex separation constants $\omega_{1},\omega_{2}$, and then determine the solution to the resulting ordinary differential equation for $R(r)$. In the case that the constants $\omega_{1}=in$, and $\omega_{2}$ is real $\omega_{2}=k$ or purely imaginary $\omega_{2}=ik$, this approach results in the same family of harmonic functions as the standard procedure. This alternative method is used by Waldron and Overfelt Waldron (1958); Overfelt (2001) to obtain separable solutions to the wave equation and Laplace’s equation in the helical coordinate system.

II.2 Helical Coordinate Systems

An idealized infinite helix is uniquely defined by the radius $R>0$, pitch111The pitch, $P$, is defined as the advance along the axis of the helix for one complete turn. $P>0$, and handedness ($+1$ or $-1$ corresponding to right-handed and left-handed helices, respectively). We use the convention that a right-handed helix is one in which moving along the helix in the counter-clockwise direction (around the $z$-axis when viewed from above) produces movement in the positive $z$ direction.

When considering an idealized helix in space, it is convenient to work in cylindrical coordinates and orient the helix so that it is contained in the cylinder given by the equation $r=R$, with the $z$-axis serving as the axis of the helix. By rotating the helix (or the coordinate system), we may assume that the helix intersects the $z=0$ plane on the positive $x$ axis, at the point $(r,\theta,z)=(R,0,0)$.

With these conventions, each right-handed helix may be parametrized by the following equations, where $\overline{p}=P/2\pi$ and $t$ is a parameter representing the polar angle $\theta$:

[TABLE]

Any parametrization of a helix in Cartesian coordinates will require that all three coordinates depend on the parameter $t$, and only by orienting the helix appropriately in a cylindrical coordinate system can we parametrize the helix so that only two of the coordinates depend on that parameter. For applications in which there is helical symmetry, it is convenient to choose a coordinate system in which a helix can be parametrized using only one of the coordinates.

The following helical coordinate system $(r,\phi,\zeta)$ was introduced by Waldron Waldron (1958), and has many similarities with a cylindrical coordinate system:

[TABLE]

We need the Jacobian of the change of variables for the above coordinate transformations, from which all of the identities concerning the coordinate directions, partial differential operators, etc. are derived:

[TABLE]

Figure 1 shows the three coordinate directions passing through a point $(r_{0},\phi_{0},\zeta_{0})$. The helical $r$ coordinate direction is identical with it’s $r$ counterpart in cylindrical coordinates. While the $\theta$ coordinate direction would wind around the cylinder in a circle, the $\phi$ coordinate direction winds up the cylinder along a helical path. It is important to note that although $\theta=\phi$ for each point in space, the vectors that point in the $\theta$ and $\phi$ coordinate directions are different at each point. To be more explicit, the vector (in cylindrical coordinates) which points in the $\phi$ coordinate direction is the second column in Eq. 4, which has a component in the $\theta$ direction and a component in the $z$ direction. Since the $\zeta$ coordinate differs from $z$ only by a translation, they are parallel: $\hat{z}=\hat{\zeta}$. The $\zeta$ coordinate can be thought of as the distance in the $z$ direction between the point and the surface $z=P\theta/2\pi$, illustrated by the length of the line segment Fig. 2.

With these coordinates, we may parametrize a right-handed helix with radius $R$ and pitch $P$ as follows:

[TABLE]

There is another difference between the cylindrical and helical coordinate system, relating to non-unique representations of points. In cylindrical coordinates, each point $(r,\theta,z)$ can be equivalently represented as $(r,\theta+2\pi k,z)$ for any integer $k$. In helical coordinates the following more complicated identity holds for each point $(r,\phi,\zeta)$ and each integer $k$:

[TABLE]

This identity captures the fact that $\Delta\phi=2\pi k$ is equivalent to $\Delta\zeta=k{P}$ (i.e., changes in the $\phi$ coordinate which are integer multiples of $2\pi$ are equivalent to changes in the $\zeta$ coordinate which are the same integer multiple of the pitch ${P}$). Informally, one gets to a point $k$ pitch-lengths above by moving in only the $\phi$ coordinate or only in the $\zeta$ coordinate. This identity is a consequence of the nonorthogonality of the $\phi$ and the $\zeta$ coordinate directions, represented by the non-zero off-diagonal term in the Jacobian for the change of variables between cylindrical and helical; relations like this are not possible in orthogonal coordinate systems.

II.3 Helical Harmonics

In Ref. Waldron (1958), Waldron introduces the right-handed helical coordinate system $({r,\phi,\zeta})$ as a method for analytically solving certain classic problems in electromagnetism in the case of helical symmetry. He recasts Laplace’s equation and the wave equation into helical coordinates, and then uses a modified separation of variables technique to write down explicit analytic solutions to the wave equation.

Using the same technique, Overfelt Overfelt (2001) solves Laplace’s equation in all of space and arrives at explicit forms for harmonic functions which separate in the helical coordinate system, hereafter referred to as helical harmonic functions. Here we give an overview of the technique, which parallels the method used for cylindrical coordinates.

As $(r,\theta,z)$ are functions of $(r,\phi,\zeta)$ through the coordinate transformation equations, the chain rule

[TABLE]

together with the appropriate entries of the Jacobian (Eq. 4), allow us to recover the following relation among the differential operators in the coordinate directions:

[TABLE]

Using these relations, we can translate Laplace’s equation from cylindrical coordinates into a helical coordinate system:

[TABLE]

If we look for separable solutions to Laplace’s equation of the form $R(r)P(\phi)Z(\zeta)$ where $P(\phi)=\exp(\omega_{1}\phi)$ and $Z(\zeta)=\exp(\omega_{2}\zeta)$, then $R(r)$ satisfies the following differential equation:

[TABLE]

Equation 5 also results from the choice $P(\phi)=\exp(-\omega_{1}\phi)$ and $Z(\zeta)=\exp(-\omega_{2}\zeta)$. Similarly, by choosing $P(\phi)=\exp(\omega_{1}\phi)$ and $Z(\zeta)=\exp(-\omega_{2}\zeta)$, or $P(\phi)=\exp(-\omega_{1}\phi)$ and $Z(\zeta)=\exp(\omega_{2}\zeta)$, we find that $R(r)$ satifies

[TABLE]

As in the cylindrical case, the ordinary differential equations for $R$ can be transformed into Bessel equations of complex order $i(\omega_{1}-\overline{p}\omega_{2})$ in Eq. 5 and $i(\omega_{1}+\overline{p}\omega_{2})$ in Eq. 6 in the variable $s$ through the linear change of variables $s=\omega_{2}r$. Thus, for arbitrary complex numbers, $\omega_{1}$ and $\omega_{2}$, we obtain the following separable helical harmonics $\Phi_{\omega_{1},\omega_{2}}({r,\phi,\zeta})=R(r)P(\phi)Z(\zeta)$, in terms of Bessel functions of complex order and argument:

[TABLE]

Although the helical harmonics separate in the variables ${r,\phi,\zeta}$, when we add Eqs. 7 and 8, the additive property of the exponential function allows us to recover functions which couple the $\phi$ and $\zeta$ variables inside the argument of trigonometric functions. Converting to cylindrical coordinates gives functions of the form $R(r)F(\theta,z)$ which does not separate in the cylindrical variables $r,\theta,z$.

In particular, choosing $\omega_{1}=in^{*}$ and $\omega_{2}=ik$ in Eqs. 7 and 8, where $n^{*}$ is an integer and $k$ is a real constant, expressing the complex exponentials in terms of trigonometric functions and the Bessel functions in terms of modified Bessel functions, we arrive at the following form of real-valued helical harmonics, which we also represent in cylindrical coordinates and refer to as “left-handed” helical harmonics:

[TABLE]

The same choice of $\omega_{1}=in^{*}$ and $\omega_{2}=ik$ in Eqs. 9 and 10 results in the following “right-handed” helical harmonics:

[TABLE]

Again, we have a family of harmonic functions indexed by a real parameter $k$. For notational ease, we are free to set $n=n^{*}-\overline{p}k$, and since $k$ controls the oscillations of $\Phi$ with respect to $z$, by setting $k$ to be integer multiples of the effective pitch, $k_{m}=m/\overline{p}$, we recover a countable family $\{\Phi_{n,k_{m}}\}$ of helical harmonic functions which respect the periodicity of the solenoidal coil in $z$.

On any cylinder $r=r_{0}$, each right-handed helical harmonic function takes the following form, where $R(r)$ is a modified Bessel function:

[TABLE]

The expression in Eq. 11 is constant along the curves in the cylinder $r=r_{0}$ that satisfy $n\theta-k_{m}z=$ const. That is, it is constant along curves $z=n\theta/k_{m}+$ const, which are all the right-handed helices with pitch $P^{*}=2\pi n/k_{m}$. If we work in a helical coordinate system with pitch $P$, we can make $P^{*}=P$ by choosing $k_{m}=2\pi m/P$ and letting $n=m$. It is also true that a left-handed helical harmonic function is constant along certain left-handed helices, whose pitch is determined as above. This tells us that the family $\{\Phi_{n,k_{m}}\}$ introduced above contains a sequence of terms $\Phi_{n,k_{n}}$ which are all constant along right-handed helices with pitch $P$. In Fig. 3 we show cylindrical and helical harmonic functions restricted to a cylinder $r=r_{0}$.

If the field ${\bf B}$ is due to an idealized wire, parametrized by ${\bf h}(t)$ (Eq. 3), carrying a uniform current $I$, then the field is given by the Biot-Savart Law:

[TABLE]

Theorem 1: A field ${\bf B}$ given by the Biot-Savart Law for an idealized infinite right-handed helix with pitch $P$ leads to a field which is cylindrically symmetric along right-handed helical curves with pitch $P$, that is ${\bf B}(r_{0},\theta_{0},z_{0})={\bf B}(r_{0},\theta_{0}+2\pi s,z_{0}+Ps)$ for all $(r_{0},\theta_{0},z_{0})$ and all $s$, when ${\bf B}$ is represented in cylindrical coordinates ${\bf B}=B_{r}\hat{r}+B_{\theta}\hat{\theta}+B_{z}\hat{z}$.

Proof.

For each real $s$, the change of coordinates $C_{s}(r,\theta,z)=(r,\theta+2\pi s,z+Ps)$ is an orientation preserving rigid motion of space under which infinite right-handed helices of pitch $P$ are invariant. $C_{s}$ is a translation of the cylindrical coordinate system along a helix with pitch $P$, and is the composition of a rotation about the $z$-axis by an angle $2\pi s$ and a translation in the $z$ direction with $\Delta z=Ps$. This isometry $C_{s}$ also acts on vectors via the Jacobian $JC_{s}=\frac{\partial C_{s}}{\partial(r,\theta,z)}$, sending a vector ${\bf V}_{0}$ at ${\bf x}_{0}$ represented in cylindrical components, to the vector ${\bf V}_{1}=JC_{s}({\bf x}_{0}){\bf V}_{0}$ at ${\bf x}_{1}=C_{s}({\bf x}_{0})$. Since $C_{s}$ is a translation in cylindrical coordinates, the Jacobian is the identity matrix, and so the vectors ${\bf V}_{0}$ and ${\bf V}_{1}$ have the same cylindrical components. Now, since this rigid motion leaves the geometry of the current carrying device unchanged, since $C_{s}({\bf h}(t))={\bf h}(t+s)$, the resulting magnetic field must be invariant under the transformation $C_{s}$: ${\bf B}(C_{s}({\bf x}_{0}))=JC_{s}({\bf x}_{0}){\bf B}({\bf x}_{0})={\bf B}({\bf x}_{0})$, and the result follows. ∎

In the special case when $s$ is an integer, the translation in $\theta$ is an integer multiple of $2\pi$, so the transformation $C_{s}$ is equivalent to a pure translation in the $z$ direction with $\Delta z=sP$, an integer multiple of the pitch. In this case, the result of Theorem 1 is equivalent to the statement that ${\bf B}$ is $P$-periodic in $z$. Theorem 1 implies the more general statement that in any parallel planes $z=c_{1}$ and $z=c_{2}$, the field ${\bf B}$ is identical, up to a rotation by $\Delta\theta=2\pi(c_{2}-c_{1})/P$.

Theorem 1, together with the preceding properties of the real-valued helical harmonics, leads to our choice of right-handed helical harmonic functions to model a field due to a right-handed helical solenoid, since these are precisely the harmonic functions which are constant along curves with the handedness and pitch of the solenoid.

III Functional Forms

This section details the individual functional forms used during the fitting process. Functions are expressed in either cylindrical or Cartesian coordinates, as these coordinate systems are typically used when measuring magnetic field points in solenoids, and allow for more straightforward comparisons between different forms. Functional forms are broken into three classes: “cylindrical”, “helical”, and “additional.” These classes differentiate between functions derived from series solutions to Laplace’s equation in cylindrical coordinates, series solutions to Laplace’s equation in helical coordinates, and additional functions that help with specific sources of magnetic field asymmetry.

III.1 Cylindrical Functions

Magnetic fields stemming from idealized solenoids exhibit monotonic behavior for all field components as a function of $r$. This motivates the use of modified Bessel functions as a general solution to solenoidal magnetic fields, as modified Bessel functions exhibit monotonic behavior, while regular Bessel functions exhibit oscillatory behavior. From that equation, we can implement a fully real-valued series expression using $2(m\cdot n)+n$ free parameters:

[TABLE]

The coefficients $A_{m,n}$, $B_{m,n}$, and $C_{n}$ are determined by the fitting procedure described in Section IV. The quantity $k_{m}$ is equal to $(m+1)\pi/L$, where $L$ is some effective length scale and $m+1$ is used instead of $m$ in order to avoid the trivial constant term which is better handled as described in Section III.3.

Physically, the scalar potential $\Phi(r,\theta,z)$ cannot be observed or measured directly. Instead, a least-squares fit to the magnetic field components is used to determine the free parameters $A_{m,n},B_{m,n},C_{n}$:

[TABLE]

where $I^{\prime}_{m}=\frac{d}{dr}\left[I_{m}\right]$ is the ordinary derivative of the modified Bessel function.

The expressions in Eq. 14 can be carried out to any number of terms in $m$ and $n$. In practice, the truncation of these series is determined by the required modeling accuracy. In general, more asymmetric features in the $z$-direction necessitate a larger number of $m$ terms, while more axial asymmetry necessitates a larger number of $n$ terms.

III.2 Helical Functions

As expected from the cylindrical functions, there is no direct coupling between the $\theta$ and $z$ variables in Eqs. 13 and 14. In order to account for axial asymmetries induced by the helical coil, one would have to expand the cylindrical series to a very large number of terms. In contrast, the helical harmonic functions link the $\theta$ and $z$ coordinates, and choosing the appropriate pitch parameter will lead to a large reduction in free parameters needed to obtain an accurate modeling of the magnetic field. As described in Section II.3, a right-handed helix leads to surfaces of constant potential via a right-handed series solution. The following examples are based solely on right-handed coils, so only the right-handed terms are used for the helical series expansion. From Eq. 11, the scalar field can be expressed in the following manner:

[TABLE]

where $h_{m}=(m+1)/\overline{p}$. As before, we use $m+1$ instead of $m$ to avoid the trivial constant term. The terms $D_{m,n}$ and $E_{m,n}$ are constants. We obtain the expressions for the magnetic field components:

[TABLE]

III.3 Additional Functions

The helical and cylindrical harmonic series expansions are constructed in the context of an infinitely long solenoid. Therefore, they are not equipped to handle constant fields, nor the edge effects from a finite solenoid; two types of behavior we expect to observe.

Constant fields can be handled simply by using a first-order Cartesian solution to Laplace’s equation. For the examples used in this paper, the solenoid axis runs along the $z$-axis, and there is a near-constant $B_{z}$ field component. This implies that a $\Phi(x,y,z)=cz$ term in the potential expansion can represent that field appropriately. In practice, by including the $cz$ term we are able to reduce the total number of free parameters needed to accurately model the magnetic field. In particular, this term also replaces the zeroth order Bessel function terms ($I_{0}(r)$), which are the only terms capable of producing non-zero values for the magnetic field at $r=0$.

The final set of terms used in the following fitting procedure are derived directly from the Biot-Savart Law (Eq. 12). An infinitesimal current located at $(x_{0},y_{0},z_{0})$ element can be modeled using the following form:

[TABLE]

where ${\bf F}=(F_{x},F_{y},F_{z})$ are the free parameters associated with an effective infinitesimal current element, and ${\bf r}=(x-x_{0},y-y_{0},z-z_{0})$. These terms are included to account for small asymmetric contributions arising from the initial and terminal endpoints of the solenoid simulations, detailed in the following section. Two of these infinitesimal current element terms are included during fitting, with the spatial free parameters ($x_{0},y_{0},z_{0}$) constrained to resided within $1\,\textrm{cm}$ of the respective endpoints.

IV Calculation and Modeling

IV.1 Calculating a Helical Solenoid

In order to determine the benefit of a helical solution to Laplace’s equation, we simulate magnetic fields due to helically-wound solenoids. Using MATLAB MATLAB (2018), we define three helical coil geometries which would be typical of high-field, large-bore solenoids used in particle physics detectors. The coils consist of flat, ribbon-like cables (Fig. 4) with a constant and homogeneous current and without any additional leads or connections. The cables have a width of 5.25 mm and a height of 2.34 mm. The magnetic fields they generate are calculated through use of the Biot-Savart law, in which a 3D, parallelized numerical integral is computed to obtain the magnetic field at a grid of points inside the coil. When performing the integral, the angular step size, $\Delta\alpha$, between cross-sectional planes of the cable is roughly $5.8\cdot 10^{-4}$ radians. The parameters used for the different fields generated for this study are shown in Table 1.

The geometries of the three different solenoids are chosen to probe different magnetic field features. Solenoid A is long (92 m) compared to its radius (25 cm). It also has a large pitch (10 cm), which emphasizes the helical field features. These properties lead to a field that exhibits ideal helical behavior in the center region of the bore, with little impact from edge effects. Solenoid B is both shorter (9.2 m) and has a larger radius (1 m) than solenoid A. It also has a smaller pitch (7.5 mm), which reduces the helical field features to a negligible degree. This field can be modeled with cylindrical terms only. Solenoid C has geometric properties from both A and B. It has the same radius and pitch as Solenoid A, but with a length equal to Solenoid B. To model this field appropriately, we use a combination of helical and cylindrical terms.

Depictions of the field components for Solenoid C are given in Fig. 5. These magnetic fields have a relatively large, constant field in the $z$-direction, as is characteristic of solenoids. The radial and axial components are orders of magnitude smaller than the $z$-component (as long as one is far from the edges of the solenoid), but they are not negligible. Due to the helical nature of the solenoid, all field components exhibit wiggles in the $z$-direction, with a wavelength given by the pitch of the solenoid.

IV.2 Fitting Magnetic Field Data

The magnetic fields can be generated using as fine a grid as required (or that computing power allows). However, a real-world solenoid cannot be surveyed on an arbitrarily fine grid. Instead, a solenoid is typically surveyed sparsely, and from that survey, the full magnetic field is reconstructed. In this section, we use the harmonic functions derived above to fit a sparse set of magnetic field values, and reconstruct a continuous field using those functions. Judicious choices of the number of parameters and length scales allow us to model a range of magnetic fields with a high degree of accuracy.

The free parameters are determined by fitting $(B_{r},B_{\theta},B_{z})$ from a subset of the grid points. The remaining grid points are used for validation and visualization. The sparse subset of grid points used for fitting is selected by mimicking the process of surveying a real-world solenoid using a mechanical field mapper Feher et al. (2018). This results in a cylindrically symmetric grid with constant spacing in $r,\theta,z$, the exact values of which are detailed in the following subsections.

These grid points and their associated magnetic field components, ${\bf B}(x_{i},y_{i},z_{i})$, are used as independent variables for the fit. The sum of squared residuals, $\chi^{2}$ (Eq. 2), is the objective function which is minimized using least-squares optimization algorithms from the Scipy and lmfit packages for the Python software language Jones et al. (01); More (1978); Branch et al. (1999); Newville et al. (2014), and linear least-squares solver in MATLAB MATLAB (2018).

In addition to either the partial derivatives of the cylindrical (Eq. 13) and/or helical (Eq. 15) series expansions, the additional terms discussed in Section III.3 are included with all fitting functions.

To determine the optimal number of terms in a given series expansion, fits are performed iteratively, increasing the number of $m$ and $n$ terms at each iteration. Once the desired accuracy is reached, no additional terms are added, even if additional terms leads to slight improvements with respect to the quality of fit. The accuracy requirements for the three solenoid examples are discussed in more detail in the following sections.

Additional studies are performed on Solenoid C, in order to investigate the impact of measurement error and grid density on the fit performance. Gaussian noise is introduced to the sparse field component values to emulate measurement error. The density of the sparse grid is varied, reducing the total number of data available during fitting.

IV.2.1 Solenoid A

Solenoid A (Table 1) is intended to bring out the helical features that arise from the winding of the solenoid, while minimizing edge effects. The solenoid is very long with respect to its radius, has a relatively large pitch, and region of interest is restricted to the middle of the solenoid, far from the edges. This region is fit with two different series expansions: the cylindrical expansion (Eq. 14) and the helical expansion (Eq. 16) in order to ascertain which expansion leads to greater accuracy and fewer terms. The two expansions are compared using the objective function (Eq. 2), the per-component residuals, and the overall accuracy using the validation data set.

The sparse set of data selected to fit for Solenoid A consists of 6 radial positions, chosen such that $12.5\,\textrm{mm}<r<150\,\textrm{mm}$, 16 equally-spaced angles from 0 to 2$\pi$, and 60 equally-spaced $z$ positions $-1500\,\textrm{mm}<z<1500\,\textrm{mm}$. The total number of data points is 5760, and each point contains three magnetic field vector components, $B_{r},B_{\theta},B_{z}$. The validation dataset is generated from a selection of 24 equally spaced radial positions from $6.25\,\textrm{mm}<r<150\,\textrm{mm}$, 126 equally spaced angles, and 120 equally spaced Z positions, leading to a total of 362,880 data points. These points were chosen to probe a region far from the edges of the solenoid to minimize edge effects, and far enough from the solenoid coil to avoid discontinuities while still capturing the helical wiggles. The pitch, $P$, chosen for the helical series is 10 cm, identical to the pitch used to generate the solenoid. The length scale, $L$, used for the cylindrical series is 5 cm. This length scale was chosen such that $k_{m}=h_{m}$ for Eqs. 14 and 16, which determines the wavelength of the wiggles in the $z$-direction.

An example of the fits, with associated residuals, are shown in Figs. 6 and 7. The 3D plots on the left-hand side show the values of the field components with both the calculated data (black points) and the fit function (green mesh) for a single 2D plane of measurements. Note that the parametric fit was determined using a full 3D selection of points, not just the 2D plane displayed. In all three components, the wiggles due to the helical coil are significant, especially at the outer radius of the fit. The size of the wiggles from peak to trough are roughly 3-4 G for the $B_{z}$ component and roughly 0.2-0.3 G for the $B_{\theta}$ and $B_{r}$ components. The residuals, expressed as Data-Fit in units of Gauss, are shown on the right-hand sides of the plots. The vertical striping (especially prominent in the upper right plot of Fig. 6, for example) corresponds to each of the radial positions used in the sparse data set, as the residuals are only evaluated at each of the fitted points for these plots. The horizontal striping corresponds to the variation in the residuals as a function of the $z$ step size. While the left-hand plots for the helical and cylindrical fits may look equivalent, the upper and lower bounds on the residuals for the helical fit are between two to three orders of magnitude smaller than the those of the cylindrical fit. The fit to helical harmonic functions is much more accurate than the fit to cylindrical functions.

The quality of the fits as a function of the total numbers of terms $m$ and $n$ is shown in Fig. 8. For both the cylindrical and helical expansions, the $m$ terms multiply the $z$ coordinate inputs and the $n$ terms multiply the $\theta$ coordinate inputs. The helical series requires only a few terms in the expansion before it reaches a point where additional terms do not greatly improve the quality of fit. In contrast, the cylindrical series continually improves as the number of terms increases, until determining the optimal free parameter values becomes too computationally taxing. Even at roughly 10 times the number of free parameters, the cylindrical expansion does not reach the level of fit quality found using the helical expansion.

The residuals derived from the validation dataset are shown in Fig. 9 for both the helical and cylindrical series. Note that the range of the horizontal axis for the cylindrical and helical plots are vastly different; the range for the helical plot is about six orders of magnitude smaller than the cylindrical. At its worst, the range of the residuals for the helical series are at least four orders of magnitude better than the residuals for the cylindrical series, and with roughly 10 times fewer free parameters. As expected, the helical series is the superior choice for this example, and the relative residuals on the overall magnitude of the field are typically much better than 0.1 parts per million.

The set of terms chosen to construct the validation plots (Fig. 9) was $m\leq 2$ and $n\leq 3$. Additional terms only added a one to two percent improvement to the $\chi^{2}$ statistic, as can be seen from Fig. 8. In Section II.3, terms in which $n=m$ (or $n=m+1$, due to the redefinition in Section III.2) are shown to lead to solutions for constant helical surfaces for a given pitch for a right handed helix. The results in Fig. 8 confirm that prediction; there is a large improvement in the residuals for both the $(m,n)\leq(1,2)$ terms and the $(m,n)\leq(2,3)$ terms. The cylindrical series used for validation used $m\leq 1$ and $n\leq 21$, as additional terms produced no improvement to the $\chi^{2}$ metric or the residual validation plots.

IV.2.2 Solenoid B

Solenoid B differs from Solenoid A in three distinct ways: it has a larger radius (1 m), a smaller pitch (7.5 mm), and a shorter length (9.2 m) than Solenoid A. The field has a roughly 1 T $B_{z}$ component, and $B_{r}$ and $B_{\theta}$ components that are 3-5 orders of magnitude smaller. These properties lead to a magnetic field which does not exhibit helically-induced wiggles in the field except very close to the coil. Additionally, the large radius and small pitch lead to a large radius-to-pitch ratio, which presents challenges to the helical harmonic series. A very large ($\gg 10$) radius-to-pitch ratio is difficult to use because it is the argument of the modified Bessel functions used in the helical series expansion. This large argument causes the Bessel functions to grow until they are computationally unwieldy ($>10^{30}$). The shorter length of the solenoid impacts the uniformity of the field components as a function of $z$, such that they vary smoothly over the length of the solenoid. In order to address all these features, the cylindrical series expansion is used instead of the helical series, where the input length scale, $L$, is the total length of the solenoid. The results of this example are used to motivate the parameterization of the fit to Solenoid C, which is a hybrid of Solenoids A and B.

The sparse data set used for fitting consists of the same $\theta$ and $z$ selection used for Solenoid A, along with 6 radial positions between 15 and 90 cm, in 15 cm increments. Examples of the fits to Solenoid B are shown in Fig. 10. As for the to Solenoid A fits, the horizontal and vertical striping in the residual plots is a consequence of the discrete nature of fit data. The oscillations in the residuals, especially prominent in the $B_{z}$ residuals, are due to the shorter wavelength $m$ terms used in the parameterization. If more terms were added to the parameterization, the magnitude of these oscillations would be reduced. To determine the number of terms to include in the series expansion, the fits were conducted iteratively, increasing the number of terms with each iteration. The results of these fit iterations are shown in Fig. 11. Based on these results, the series is expanded to $m\leq 5$ and $n\leq 3$, as additional terms only lead to diminishing (1-2%) improvements to the residuals.

The validation set used to determine the per component residuals for this model is identical in $\theta$ and $z$ with respect to Solenoid A, and used 24 equally spaced radial steps between 3.75 cm and 90 cm. The residuals based on the validation data are shown in Fig. 12. The residuals are never larger than $5\cdot 10^{-3}$ G, and are typically much smaller. Overall, the residuals are smaller than 1 part per million. The skewness of the residual plot for $B_{\theta}$, which affects about 1% of the validation dataset, is driven by an unavoidable asymmetry due to first and last turns of the solenoid, which deviate from an ideal cylinder. This can be seen in the $B_{\theta}$ residual plot of Fig. 10, specifically the darker blue regions in the upper and lower left corners corresponding to the residuals at large radius and $|z|>1.0$ m.

IV.2.3 Solenoid C

Solenoid C is identical to Solenoid A in all respects except for length; it is the same length as Solenoid B, (10 times shorter than Solenoid A), which enhances the $z$-dependent, $\theta$-independent features of the magnetic field due to the finite nature of the solenoid. As can be seen in Fig. 5, Solenoid C has both the large scale, smooth variation of Solenoid B in addition to the high frequency wiggles of Solenoid A. To model this simulation appropriately, a combination of terms is used from both the helical and cylindrical series expansions.

The sparse set of data selected to fit for Solenoid C, as well as the set of validation data, are identical to the sets used for Solenoid A. Instead of performing a scan over a very large potential space of terms, the number of terms for Solenoid A ($m\leq 2$ and $n\leq 3$ for helical) is combined with the number of terms for Solenoid B ($m\leq 3$ and $n\leq 5$ for cylindrical), and all parameters are free. Examples of the fit to Solenoid C are shown in Fig. 13. The residual plots indicate that the fields are modeled very well, with only small ($<10^{-4}$ G) variations between fit and data. The residuals based on the validation dataset, shown in Fig. 14, are on the order of $10^{-4}$ G or less for all field components, and the residuals for the total field magnitude are better than one part per million, which indicates that this combination of terms works well for modeling all major features of this magnetic field.

Two further studies were performed on Solenoid C. The first study investigated the effect of varying the density of the sparse data used for fitting, and the second study explores the effect of measurement error on the fit data.

In the first study, the original sparse data for fitting Solenoid C is a cylindrical lattice whose coordinates are one of 6 radial positions, 16 angular positions, and 60 $z$ positions. The total number of data points is 5760. We can specify subgrids by selecting subsets of each of the radial, angular, and $z$ positions.

Fig. 15 shows the residuals for fits using subgrids with 1/3, 1/8, and 1/9 of the density of the original fitting grid. The 1/3 and 1/9 density subgrids are formed by selecting every $3^{rd}$ and every $9^{th}$ point among the $z$ positions, respectively. The 1/8 density subgrid is formed by taking every other point in each of the radial, angular, and $z$ positions. It is clear that the error increases as the density decreases, but even the largest relative error in magnitude is still on the order $10^{-6}$ for the sparsest grid. We also note that the original grid is already a very coarse subgrid of the fine grid used for validation. The details of these grids are shown in Section IV.2.3.

[TABLE]

The second study examines the response of our method to measurement error on the fit data. For modern solenoid mapping, NMR and Hall probes are required to have measurement accuracy on the order of $10^{-5}$ hall_probe. In order to simulate this type of measurement error, Gaussian noise is added independently to each field component in the data before fitting. The Gaussian mean is 0 and the standard deviation in each component is chosen to simulate a relative error of $\varepsilon$. The residuals are calculated by comparing the fit values to the true field values ($\varepsilon=0$).

Fig. 16 shows histograms of the residuals for fits with $\varepsilon=10^{-3},10^{-4},10^{-5},0$. The histogram for $\varepsilon=0$ corresponds to the nominal fit as seen in Fig. 14. The $\varepsilon=10^{-5}$ error leads to a relative residual of better than $10^{-6}$ for the majority of data points, and even the most pessimistic error simulations are typically less than $10^{-5}$. This implies that our method is robust with respect to expected measurement errors for magnetic field surveys.

V Summary and Conclusions

Starting from first principles, we have introduced general solutions to Laplace’s equation in a helical coordinate system for the purposes of modeling realistic solenoidal magnetic fields. Using sparse sets of data generated for different solenoidal geometries, we modeled the magnetic field using expansions of helical and cylindrical harmonic functions and validated the model with a high-granularity magnetic field calculation. Using a small number of free parameters, we were able to obtain accurate magnetic field models, and capture asymmetric behavior caused by the helical coil. The residuals are better than one part per million, and typically better than 0.1 parts per million.

The helical series expansion is good at capturing the wiggles that arise from the helical winding of a solenoid. Depending on the pitch and radius of the solenoid, these oscillations can become significant compared to the average field strength, and ignoring these contributions would lead to an incorrect model of the magnetic field. The cylindrical series are useful for modeling the low frequency components of a finite-length solenoid. Using an effective length equal to the solenoid length, the cylindrical series models the $z$-dependent, axially symmetric features of the magnetic field. For solenoids of both finite length and with helical features, a combination of helical and cylindrical series solutions provides the most accurate model. While not exhaustive, we examined the impact of realistic measurement errors on our model. Even with pessimistic error assumptions, our model still performs well, with residuals better than 10 parts per million, and typically better than 1 part per million. Further studies of magnetic field perturbation due to machining tolerances, external magnetic field sources (e.g. due to additional circuitry), and other sources of error are left for future investigation.

While this method was shown to be useful for representing a magnetic field, it can also be used to calculate a magnetic field. The calculated fields used in this paper were generated using a common numerical integration of the Biot-Savart law. This process, even when parallelized across GPUs, is computationally taxing, and its resource usage scales both as a function of the size of the solenoid and the number of grid points in which one would like to evaluate the field. Conversely, starting with a closed-form solution of a complex magnetic field scales to the number of grid points, and the evaluation of each grid point requires executing only a handful of mathematical terms instead of a series of numeric integrals. While the method for determining the coefficients of this series without the use of fitting is not explored in this paper, it should be possible for an appropriately specified field calculation to do so.

Acknowledgements.

The authors would like to thank Matt Bonakdarpour for helpful discussions. We gratefully acknowledge the support provided by the Department of Energy under award number DE-SC0015910. This document was prepared by members of the Mu2e Collaboration using the resources of the Fermi National Accelerator Laboratory (Fermilab), a U.S. Department of Energy, Office of Science, HEP User Facility. Fermilab is managed by Fermi Research Alliance, LLC (FRA), acting under Contract No. DE-AC02-07CH11359. We thank the reviewers for providing helpful comments on earlier drafts of the manuscript.