Implementation of Maxwell's equations in the reconstruction of the magnetic field in the $g-2$ storage ring

We present a method for implementing the constraints that are implied by Maxwell's equations in fits to measurements of the magnetic field in the muon storage ring of the $g-2$ experiment. The method that we use makes use of toroidal-harmonic solutions of Laplace's equation. We point out that the fitting problem can be approximated well as a linear-algebra problem. We have devised an efficient algorithm for the linear-algebra problem that makes it possible to find a solution for $10^5$ data points and $10^4$ harmonics in less than an hour on a present-day desktop computer. We illustrate our method by applying it to some preliminary measurements of the magnetic field in the $g-2$ storage ring.


I. INTRODUCTION
Precise knowledge of the magnetic field inside the muon storage ring is of paramount importance for the measurements of the muon anomalous magnetic moment (g − 2) that will be made by the Fermilab experiment E989 [1]. The primary method of measurement of the magnetic field, which utilizes a circular array of nuclear-magnetic-resonance (NMR) probes, determines only the magnitude of the magnetic field, not its direction. Ancillary measurements using Hall probes have been used to obtain information about components of the field that are orthogonal to its average direction [2].
As we will explain, checks of these measurements can be obtained by making use of the constraints that are implied by Maxwell's equations. Since the inside of the muon storage ring is a vacuum, Maxwell's equations for the magnetic field B become Because the magnetic field is almost uniform in direction and magnitude inside the storage ring, we can write B as a sum of the uniform field B ẑ and a fluctuation B f : whereẑ is a unit vector in the direction of the average of the B field over the active volume of the storage ring. For the E989 ring,ẑ is very close to the vertical direction. B is chosen to be about the same size as the magnetic field strength. The precise value of B is unimportant, as small changes in B can be absorbed into B f . However, as we shall see, it is crucial to our analysis that B f be small compared with B .
Since the uniform field B ẑ satisfies Maxwell's equations in vacuum, so does the fluctuation B f . The equation ∇ × B f = 0 implies that B f can be written as a gradient of a scalar function V , which satisfies Eq. (3b), which follows from ∇ · B f = 0, is Laplace's equation.
The region of validity of Eq. (3b) can be an issue because the single-valuedness of V relies on Eq. (1b), through Stokes' theorem. For the E989 muon storage ring, there are no currents through the plane of the ring at radii that lie inside the outermost part of the vacuum chamber. Consequently, V is single-valued inside the vacuum chamber.
In this paper, we implement Maxwell's equations in fits to the NMR measurements of the storage-ring magnetic field by expressing the fitting function as a sum of solutions of Maxwell's equations. We do this by writing V as an expansion in a complete set of solutions of Laplace's equation (harmonic functions). In particular, we use toroidal harmonics as the basis functions, since their coordinate geometry is particularly well matched to the geometry of a toroidal storage ring.
Given the large number of data points and harmonics and the underconstrained nature of the fitting problem, the fitting problem, in its original form, is quite formidable. We have developed several simplifications/reorganizations of the fitting algorithm that make the fitting problem computationally tractable.
We have applied our fitting algorithm to preliminary data from the g − 2 trolley runs.
While these data do not reflect the ultimate precision of the NMR measurements or the ultimate quality of the g − 2 storage-ring magnetic field, they serve to illustrate the fitting procedure and its utility.
The remainder of this paper is organized as follows. In Sec. II we describe the toroidalharmonic expansion of the magnetic field. In Sec. III we outline our method for fitting the toroidal harmonic series to the data. Here we describe the algorithmic improvements that we have made. Section IV contains some examples of the application of our fitting procedure to some of the early NMR magnetic-field data. Finally, in Sec. V, we summarize our results.

II. USE OF TOROIDAL HARMONICS TO IMPLEMENT MAXWELL'S EQUA-TIONS
The expansion of a vacuum magnetic field in toroidal harmonics was discussed in detail in Ref. [3]. Here, we present the essential features of that expansion.
Toroidal coordinates (ζ, η, φ) are given in terms of Cartesian coordinates (x, y, z) as fol- minor radii of R/ tanh ζ and R/ sinh ζ, respectively, whose center is located at x = y = z = 0. The surface of constant η is a sphere of radius R/ sin η, whose origin is at x = y = 0 and z = R/ tan η. lows: Here, 0 ≤ ζ < ∞, 0 ≤ η ≤ 2π, and 0 ≤ φ ≤ 2π. φ is the usual azimuthal angle about the z axis, with φ = 0 corresponding to the x axis. R is a constant that has a dimension of length. The surface of constant ζ is a torus of major and minor radii of R/ tanh ζ and R/ sinh ζ, respectively. (ζ → ∞ corresponds to a circle of radius R in the x − y plane, and ζ → 0 corresponds to (x → ∞, y → ∞) and to the z axis.) The poloidal angle η specifies, for each circular φ slice of the torus, the position on the circle. The surface of constant η is a sphere of radius R/ sin η, whose origin is at x = y = 0 and z = R/ tan η. The toroidal coordinates are illustrated in Fig. 1.
Although there is some freedom in choosing the value of R, it should not be too different from the radius of the storage ring, so that a torus defined by a constant ζ roughly coincides with the interior of the storage ring. For numerical calculations, R should be adjusted so that the points of interest do not lie exactly on the circle of radius R in the x − y plane; otherwise, one encounters ζ = ∞ in coordinate transformations.
The general solution to Laplace's equation, which is found by using the technique of separation of variables in toroidal coordinates, is a linear combination of the following functions, which are defined for nonnegative integers m and n by cosh ζ − cos η P n m−1/2 (cosh ζ) cos(mη) cos(nφ), cosh ζ − cos η P n m−1/2 (cosh ζ) cos(mη) sin(nφ), cosh ζ − cos η P n m−1/2 (cosh ζ) sin(mη) cos(nφ), (5c) where the toroidal harmonics P n m−1/2 (cosh ζ) and Q n m−1/2 (cosh ζ) are half-integer Legendre functions of the first and second kind, respectively. P n m−1/2 (cosh ζ) and Q n m−1/2 (cosh ζ) can be expressed in terms of the hypergeometric function as [3,4] The P n m−1/2 (cosh ζ) correspond to the solutions that are singular at ζ → ∞, and the Q n m−1/2 (cosh ζ) correspond to the solutions that are regular at ζ → ∞. Therefore, harmonics containing P n m−1/2 (cosh ζ) correspond to contributions from sources in the interior of the ring, and harmonics containing Q n m−1/2 (cosh ζ) correspond to contributions from sources outside the ring.
In practical calculations, it is necessary to work with a finite number of harmonics. This may result in bad convergence properties if the potential V is not single valued 1 . The 1 For example, the Fourier expansion for the function f (θ) = θ for 0 ≤ θ ≤ 2π does not converge uniformly. potential V may not be single valued when the region of space that satisfies Eq. (3b) is multiply connected. For example, the condition V (ζ, η + 2π, φ) = V (ζ, η, φ) can be spoiled if there is a current that flows inside the ring because the line integral C dx · B along a closed curve C does not vanish if C encloses a nonzero net current. For the muon storage ring, the condition V (ζ, η + 2π, φ) = V (ζ, η, φ) holds because, in the absence of the muon beam itself, the interior of the ring is a vacuum. The condition V (ζ, η, φ + 2π) = V (ζ, η, φ) is ensured by the fact that C dx · B = 0, where the path C is on the interior of the ring and encloses the z axis. The vanishing of this line integral follows from the fact that no net current flows through the plane of the ring at radii that lie inside the outermost part of the vacuum chamber.
We set the cutoffs for m and n to be M and N , respectively. In order to eliminate solutions that correspond to currents that pass through the plane of the storage ring at radii that lie inside the outermost part of the vacuum chamber, we exclude the harmonics P n m−1/2 (cosh ζ) from the solution. Then, we obtain CC cos(mη) cos(nφ) + C mn CS cos(mη) sin(nφ) +C mn SC sin(mη) cos(nφ) + C mn SS sin(mη) sin(nφ)] , where, for numerical stability, we divide each toroidal harmonic by Q n m−1/2 (cosh ζ 0 ), so that Q n m−1/2 (cosh ζ)/Q n m−1/2 (cosh ζ 0 ) is always order one when ζ ≈ ζ 0 . There is some arbitrariness in the choice of ζ 0 . An appropriate choice for ζ 0 is one such that R/ tanh ζ and R/ sinh ζ are approximately equal to the major and minor radii, respectively, of the storage ring.
Once we have fixed the (2M + 1) × (2N + 1) constants 2 C mn CC , C mn CS , C mn SC , and C mn SS , we can give an approximate reconstruction of the magnetic field in the storage ring. The magnetic field components are given by the gradient of Eq. (7). The magnetic field can be written in terms of cylindrical basis vectors as 2 Since sin 0 = 0, C m0 CS , C 0n SC , C 0n SS and C m0 SS do not appear in Eq. (7) where [3] Hereafter, we refer to B ρ , B z , and B φ as radial, vertical, and longitudinal fields, respectively.

III. FITTING THE TOROIDAL-HARMONIC SERIES TO DATA
Now we discuss the strategy to constrain the magnetic field by comparing the expressions in Eqs. (7) and (9) with the measured magnetic field.

A. Linearization of the problem
We determine the constants C mn CC , C mn CS , C mn SC , and C mn SS by fitting the expressions in Eqs. (7) and (9) to the measured magnetic field strengths B measured . The most straightforward way to do this is to find a set of constants that minimizes the sum of squares of residuals S, which is defined by where x p is the position at which the measurement of the p'th data point is made, N data is the total number of measurements, and V is given by Eq. (7).
Although minimization of S is a general method with which to fit the data, it is not feasible to carry out in practice, especially when M , N , and N data are large. Each component of B = B ρρ + B zẑ + B φφ at a given position in space is a linear polynomial that depends on N constants ≡ (2M + 1) × (2N + 1) constants. Then, the field strength |B| = B 2 ρ + B 2 z + B 2 φ that appears in Eq. (10) is a square root of a quadratic polynomial that depends on N constants constants, which can contain up to 1 2 (N constants + 1)(N constants + 2) terms. Since the sum of the square roots of different polynomials does not, in general, simplify, the number of terms in S is about 1 2 N 2 constants N data . In the case of our fits to the storage-ring field data, N constants is typically of order 10 3 and N data is typically of order 10 6 . Hence, the computation of S becomes very compute-time intensive.
In order to simplify the fitting problem, we take the following approach. We expand the magnitude of the field strength in powers of B f : Since the transverse components of the magnetic field in the g − 2 storage ring are typically about 100 ppm, the corrections beyond linear order are of order 0.01 ppm, which is of the same order as the nominal error in the NMR-probe measurements. Therefore, we can reasonably ignore the corrections beyond the linear order. In the linear approximation, we have |B| = B z . Therefore, in the linear approximation, the sum of squares of residuals becomes B z is a linear polynomial containing N constants constants. Then, each term in the sum in Eq. (13) is a quadratic polynomial containing N constants constants, which can have up to 1 2 (N constants + 1)(N constants + 2) terms. When we add polynomials to perform the sum in Eq. (13), we can simply add the coefficients numerically. Thus, the number of terms in S L is 1 2 (N constants + 1)(N constants + 2). This is a factor of 1/N data smaller than the number of terms in S. By minimizing S L , we find an approximate solution to the following set of equations, which, in general, is overconstrained: where the measurement index p takes on the values p = 1, 2, . . . , N data . Since B z (x p ) is a linear polynomial that contains N constants constants, Eq. (14a) can be rewritten as where c 1 , c 2 , . . . , c Nconstants denote the constants C mn CC , C mn CS , C mn SC , and C mn SS .
A further advantage of the linear approximation is that it simplifies the minimization problem into a linear-algebra problem. Since S L is a sum of squares of linear polynomials, it is bounded from below, and, if S L has one local minimum, the local minimum coincides with the global minimum. The condition for a local minimum is for all i = 1, 2, . . . , N constants , where c 1 , c 2 , . . . , c Nconstants denote the constants C mn CC , C mn CS , C mn SC , and C mn SS . Since S L is a quadratic polynomial, the left side of Eq. (15) is a linear polynomial. Therefore, Eq. (15) can be rewritten as where The square matrix A and the column matrix b are numerical and do not depend on the c's. A is symmetric and positive semidefinite. Eq. (16) has a unique solution if A is positive definite, i.e., if A has no vanishing eigenvalues. However, as we will explain in the next section, this is not the case in our application. That is, our linear least-squares problem is underconstrained.
There are many well developed methods for solving Eq. (16). One such method is singularvalue decomposition (SVD). 3 In this method, the expression ∂B z (x p )/∂c i on the left side of Eq. (14b) is regarded as a matrix with indices p and i. Then SVD is used to find the pseudoinverse 4 of that matrix. A solution to the linear least squares problem is obtained by applying the pseudoinverse matrix to the right side of Eq. (14b). Again, as we will explain in the next section, the solution of the linear least-squares problem is not unique.
The SVD method is particularly well suited to non-unique linear least-squares problems, such as ours, because it allows one to identify the zero eigenvalues of the matrix A in Eq. (16). 3 See, for example, Ref. [5]. 4 Here, and throughout this paper, we use the term "pseudoinverse" to mean the Moore-Penrose inverse of a matrix [6][7][8][9]. See, Ref. [5] for a discussion of the pseudoinverse in the context of SVD.
However, we have found that the direct application of the SVD method to solve Eq. (16) is impractical for several reasons. First, the appearance of large numbers of different linear combinations of unconstrained c i 's makes it difficult to identify those linear combinations and to sequester them. Second, the computation, when it can be carried out, is very intensive computationally. Finally, the data display fluctuations between measurements at different azimuthal angles that are much larger than differences between the measurements of the various probes at a fixed azimuthal angle. In order to deal with these issues, we make a further rearrangement of the linear least-squares problem. In this re-arrangement, which we describe in detail in Sec. III C, we first analyze the azimuthal dependence of the data, using Fourier decomposition, and then analyze the ρ-z of the Fourier coefficients in order to determine the toroidal-harmonic coefficients.

B. Underconstrained nature of the linear least-squares problem
In the linear approximation, the measured magnetic fields constrain only B z directly.
Hence, once B z is determined, the vertical gradients of the radial and longitudinal fields can be constrained. However, aside from the combination neither the field gradients in the radial or the longitudinal directions nor the absolute sizes of B ρ and B φ can be determined from B z . This implies that a comparison of the expression for B z with the measured magnetic field strengths cannot constrain all (2M + 1) × (2N + 1) constants C mn CC , C mn CS , C mn SC , and C mn SS . We will deal with the underconstrained nature of the fits by making use of singular-value decomposition (SVD) to solve Eq. (14b). In this method, the unconstrained modes correspond to zero eigenvalues of the matrix A. We discard the associated eigenvectors from the solution. In practice, this is done by setting to zero any inverse eigenvalues in the SVD pseudoinverse matrix that correspond to eigenvalues that are smaller than a "tolerance" . While this results in a well-defined solution to the fitting problem, it is important to recognize that this is just one of an infinite number of solutions to this underconstrained problem. That means that this solution must be interpreted with care, recognizing that only the gradients with respect to z of the radial field B ρ and the longitudinal field B φ are fully determined by the constraints that are given in Eqs. (18).
C. Further rearrangement of the linear least-squares problem As we have mentioned, even the linearized fitting problem can be formidable to solve computationally, since it may involve on the order of 10 5 data points and 10 4 harmonics.
We can considerably improve the computational efficiency by carrying out the analysis in steps. The first step is to express the φ dependence of the data for each NMR probe in terms of a Fourier series. The second step is to use the Fourier-series coefficients for all of the probes to determine the coefficients of the toroidal harmonics. The first step in this approach relies on two facts: (1) the data from each NMR trolley probe is a function of φ alone (at fixed ρ and z); and (2) the toroidal-harmonic expansion of the potential V depends on φ only through sine and cosine factors in each harmonic, and B z does not involve azimuthal derivatives of V .

Fourier series for the probe data
In order to display the sine and cosine factors in the harmonics explicitly, we write so that Applying Eq. (9b) to Eq. (20), we obtain where Now we define a truncated Fourier series for each NMR probe that has the same form as the truncated Fourier series in Eq. (21): Here, the index q, which runs from 1 to N probe , labels the NMR probe, whose ζ and η coordinates are ζ q and η q , respectively. Note that we retain the index N in the arguments of the coefficients β as a reminder that these coefficients depend on the order of truncation of the Fourier series.
2. Using the probe data to constrain the Fourier coefficients A standard method for determining the coefficients of a Fourier series is to make use of the orthonormality relations for sines and cosines. In our case, application of these orthonormality relations would lead to the following equations for the Fourier coefficients β(q, n; N ): where n = 1, 2, . . . , N . Approximate values of the coefficients can be obtained by replacing B Fourier (q, φ; N ) by the measured magnetic field strengths and by replacing the integrals with finite sums that approximate the integrals. This method, which is analogous to numerical integration methods, has an obvious weakness that the samplings in φ are limited to the points at which the measurements were made. Therefore, accuracies of the Fourier coefficients that are computed by using this method deteriorate as n approaches the sampling rate of the data. A further drawback to this method is that it couples the issue of the numerical accuracy of the Fourier coefficients to the issue of the stability of the Fourier coefficients with increasing N , making it difficult to assess the latter.
The method we employ to determine B Fourier (q, φ; N ) is a least-squares fit. That is, we find the minimum of where the sum is over all values of φ at which the measurements were made, and N φ is the number of measurements. In practice we solve this linear least-squares problem by making use of the SVD method. That is, we find the pseudoinverse of the matrix whose indices are i and n, and apply it to the column vector B measured (ζ q , η q , φ i ) − B e i , whose index is i. (Here e i is a column vector with all unit entries.) We denote by χ min q (N ) the minimum value of χ q (N ): We also define It is obvious that the minimum of χ q (N ) will decrease with increasing N . Of course, this does not mean that we can increase N indefinitely. If N is too large, the Fourier coefficients become very sensitive to small fluctuations in data, which may be comparable to the uncertainties in measurement. Furthermore, the resulting Fourier series B Fourier (q, φ; N ) can exhibit unphysically large fluctuations at values of φ for which there are no data.
In order to deal with these issues, we use the following method to determine a value of N such that the Fourier series best approximates the actual magnetic field. We note that, if the Fourier coefficients are known exactly, then will decrease with decreasing N − N = ∆N for fixed N , and also for increasing N with fixed ∆N , as long as N is large enough. We will apply these criteria to determine the optimal choice of N for the Fourier analysis. When we computeχ (N, N ), we make use of Parseval's theorem, which yields where n max = max(2N + 1, 2N + 1).

Constraining the toroidal-harmonic coefficients
Having determined the coefficients β(q, n; N ) of the Fourier series in φ for each probe, we can now complete the determination of the toroidal coefficients C mn CC , C mn CS , C mn SC , and C mn SS . (Recall that we denote these toroidal-harmonic coefficients by c i .) Then we determine the c i by minimizing the following quantities for each n: where B C (ζ q , η q , n) and B S (ζ q , η q , n) are given by Eq. (22).
Again, we use the SVD method. We find the pseudoinverses of the matrices where the matrices have indices i, which runs from 1 to 2M + 1, and q, which runs from 1 to N probe . We apply the pseudoinverse of the "C" matrix to the column vector β(q, 2n; N ), which has index q; we apply the pseudoinverse of the "S" matrix to the column vector β(q, 2n + 1; N ), which also has index q. This linear least-squares problem is underconstrained, for the reasons that we discussed in Sec. III B. Therefore, the SVD procedure is very convenient for finding a particular solution.
We must carry out this procedure for each of the 2N + 1 Fourier components. However, the matrices in Eq.

IV. APPLICATION OF THE TOROIDAL-HARMONIC METHOD TO DATA
In this section, we demonstrate the toroidal-harmonic method that we have described by applying it to magnetic-field data from the g − 2 experiment. In these applications, we take B , expressed in units of the NMR frequency, to be 61.8 MHz. We take the tolerance in the SVD calculations of the toroidal coefficients to be = 10 −8 . We have found that varying over several orders of magnitude makes very little numerical difference in the gradients of the magnetic field that are well constrained by Maxwell's equations.
We emphasize that the results that we present are for illustrative purposes only. These results are based upon data that were taken early in the commissioning of the muon-storagering magnetic field, and they data do not reflect either the quality of the magnetic field or the quality of the magnetic-field measurements that will ultimately be used in the experiment.

A. Fourier decomposition
First we show results from the Fourier decomposition in φ of the data from each probe.

Run 52 data
We begin by examining data from trolley Run 52 (Ref. [10]), which are based on measurements from 25 NMR probes (located in the "shimming trolley") at 7842 azimuthal angles, for a total of 196050 data points.
In Fig. 2 We have carried out the computation of the Fourier coefficients using Mathematica code [11].
In the Mathematica code, it is very important to vectorize the matrix operations, which   27)] in ppm that we have obtained for that probe. First, we note that the values of χ min q (300) that we have obtained are much larger than the nominal probe measurement uncertainty of 10 ppb. This suggests that there may be spatial variations in the magnetic field that are on too small a scale to be represented accurately by the fit, temporal variations in the magnetic field, and/or measurement anomalies in the data. We note that the values of χ min q (300) display an interesting spatial distribution: χ min q (300) is larger at the top and bottom of the trolley than in the middle and is fairly uniform across a given horizontal row. This may provide a clue as to the origin of the unexpectedly large values of χ min q (300). In Fig. 4, we show the distributions of residuals between between the Fourier series fit and the data from Run 52 for probes 1 and 17. As can be seen, the widths of the distributions are rather different, as would be expected from the different values of χ min q (300). Figure 5 shows the Fourier series fit for probe 1 and probe 17, along with the data as a function of φ. These figures reveal significant deviations of the data from the fit, especially in the case of probe 17. One such deviation can be seen in the "zoomed-in" figure for probe 17 near 184 • . There, the data appear to be smooth in the region of the deviation.
However, the fit does not contain enough harmonics to follow the data accurately. It is not possible to increase the number of harmonics for this data set because, as we have mentioned, the fits become unstable as N increases beyond about 300. This instability can be traced to nonuniformity in the sampling of the data over the range of φ. For example, the gap in the measurements in the range φ = 201-202 • can cause instabilities once N is sufficiently large for the harmonics to probe the gap because the fit function is completely unconstrained in the range of the gap. We would expect this 1 • gap to produce instabilities when N ∼ 360 • /1 • , which is in agreement with the stability analysis that is based on

Run 3483 data
For comparison, we also show results from a later trolley run (Run 3483) in which the trolley contained 17 NMR probes. The data for this run [10] contain measurements at 9023 azimuthal angles, for a total of 153391 data points. In Fig. 6 Under each probe position, we list the value of χ min q (500) in ppm that we have obtained for that probe from the Run 3483 data.
we choose N = 500 for our fits, which implies that there are 1001 Fourier components in our analysis of the φ dependence for each probe.
In Fig. 7, we show the approximate positions in the ρ-z plane of the 17 NMR probes.
Under each probe position, we list the value of χ min q [Eq. (27)] in ppm that we have obtained for that probe. Again, the values of χ min q (500) that we have obtained are much larger than the nominal probe measurement uncertainty of 10 ppb. However, they are smaller than the values of χ q from Run 52, perhaps reflecting an improvement in the uniformity of the magnetic field as a result of the shimming process. The values of χ min q (500) display a spatial distribution that is similar to the spatial distribution of the values of χ min q (300) from Run 52. The values of χ min q (500) are larger at the top and bottom of the trolley than in the middle. However, the values of χ min q (500) are less uniform across a given horizontal row than are the values of χ min q (300) from Run 52. In Fig. 8, we show the distributions of residuals between between the Fourier series fit and the data from Run 3483 for probes 1 and 12. As was the case for the Run 52 data, the widths  the presence of large "spikes" in the data, for which a single measurement is significantly displaced from the adjacent measurements. These spikes are particularly apparent in the "zoomed-in" figure for probe 12 and, in general, they are much larger in the data for probe 12 than in the data for probe 15. This suggests that these spikes may be the reason for the large value of χ min q (500) for probe 12. It is perhaps significant that the data for Run 52 were taken while the trolley was stationary, while the data for Run 3483 were taken while the trolley was in motion.

B. Full toroidal-harmonic fit
Finally, we discuss the full toroidal-harmonic fit to the data, taking M = 8 poloidal terms for each Fourier harmonic. We have written Mathematica code to carry out the computation of the toroidal-harmonic coefficients. The computation of the coefficients for all of the Fourier harmonics takes approximately 6 minutes on a dual 6-core 2.6 GHz Intel Xeon CPU system and consumes less than 1 GB of memory in total.
In Fig. 10 We stress that, in the toroidal-harmonic fit, only ∂B ρ /∂ z is constrained, and so only the shape of the curve for B ρ is meaningful. As can be seen on comparison of the fit results with the Hall-probe data, there is considerable disagreement between the shapes. The data points at displacements of −3, +3, +4, and +5 cm from the central position were taken in different sessions than were the more central data points, and so some offset may have occurred. However, even the more central data points show much larger displacements as the distance from center is varied than does the fit. Since the uncertainties in the Hallprobe measurements are thought to be considerably smaller than these displacements, this suggests that there may be an, as yet, unknown systematic uncertainty at play in the Hallprobe measurements.
We have further tested this comparison between the toroidal-harmonic fit and the Hallprobe data by applying directly the constraint from ∇×B = 0 in Eq. (18b). We used the values of B z from the Fourier-series fit to the NMR-probe data to compute ∂B z /∂ρ = ∂B ρ /∂z; we used discrete differences to compute ∂B ρ /∂z from the Hall-probe data. The results of this computation are shown in Fig. 11. conjecture that the Hall-probe data may be subject to an unknown systematic error.

V. SUMMARY
In this paper, we have presented a method for incorporating the constraints that are imposed by Maxwell's equations into fits to the NMR-probe data for the magnetic field of the g −2-experiment storage ring. Maxwell's equations are enforced in the fits by making use of a fitting function that consists of a truncated series of toroidal harmonics. The toroidal harmonics are expressed as functions of toroidal coordinates, which are well suited to the the storage-ring geometry.
The NMR probes measure the magnitude of the magnetic field, rather than its individual components. It would be impractical to fit the toroidal harmonic series for the magnitude of the magnetic field to the data because the fitting problem is highly nonlinear and because the transverse magnetic-field components (components that are transverse to the average field direction) are not fully determined from the magnitude of the field in combination with Maxwell's equations. Therefore, we have simplified the fitting procedure by approximating the magnitude of the field by its principal (vertical) component. Given that the transverse components of the magnetic field in the g − 2 storage ring are typically of order 100 ppm, the corrections to this approximation are of order 0.01 ppm, which is comparable to the nominal error in the NMR-probe measurements.
Even this linear fitting problem is difficult computationally. Difficulties arise because of the large number of data points (of order 10 5 ), the large number of toroidal harmonics that is required for an accurate fit (of order 10 4 ), and the fact that the fit does not fully constrain the coefficients of the toroidal harmonics. A further difficulty is that the data display fluctuations between measurements at different azimuthal angles that are much larger than the differences between the measurements of the various probes at a fixed azimuthal angle.
In order to ameliorate these problems, we further rearranged the fitting procedure in to two steps. First, we fit the azimuthal dependence of the measurements of each probe to a truncated Fourier series. Then, we used the Fourier series coefficients for all of the probes to determine the coefficients of the toroidal harmonics. The underconstrained nature of the problem is confined essentially to this second step, which involves a fitting problem whose dimensionality is much lower than the dimensionality of the original linear fitting problem.
We deal with the underconstrained nature of the problem by using the singular-value decomposition (SVD) method, which identifies and isolates the underconstrained modes. This modified fitting procedure results in a very efficient computational algorithm that is easily carried out on a present-day desktop computer.
We have used this method to fit some of the early NMR-probe data that were taken by the g − 2 collaboration. While these data are not representative of the magnetic-field homogeneity or the field-measurement precision that will ultimately be used by the g − 2 collaboration to make the g − 2 measurement, they serve to illustrate some of the features of the fitting procedure.
Even at the level of the fits to the truncated Fourier series, it is apparent that the deviations of the data from the fit are much larger than the nominal precision of the NMRprobe measurements. Such fits may be useful in diagnosing sources of uncertainties in real-world field measurements. The full fits to the truncated toroidal-harmonic series are in tension with Hall-probe measurements of the radial component of the magnetic field, suggesting that there may be systematic errors in the Hall-probe measurements that are not yet understood.
We expect that the toroidal-harmonic fits will continue to be a useful tool in understanding the magnetic-field measurements in the g − 2 storage ring as those measurements evolve.
The toroidal-harmonic fits may be even more important for the planned muon electric-dipole (EDM) measurement. That measurement may be more sensitive than is the g − 2 measure-