using beam dynamics models and simulations to produce $M^{c}(x,y,t)$ for all
azimuthal locations where calorimeters ($c$) are placed. (4) Fold
$M^{c}(x,y,t)$ distributions with acceptance, asymmetry, and phase maps to
obtain phase shift versus time in fill for each calorimeter,
$\varphi^{c}_{pa}(t)$. (5) Generate Monte Carlo data using the full fit
function, including the predicted behavior of the phase,
$\varphi^{c}_{pa}(t)$. Fit these pseudodata using the normal fit function to
determine the difference in the extracted $\omega_{a}^{m}$ compared to the
input. This difference yields the correction factor $C_{pa}$ for each
calorimeter.

###  VIII.1 Measurement of the time-changing muon distribution

The two tracker stations produce the time evolution of the muon distribution
maps over the course of a fill, $M^{T}(x,y,t)$. These distributions are
corrected for resolution and acceptance as described in Sec. II.4. Betatron
oscillations cause the distribution to change rapidly, which makes extraction
of the slower terms due to the damaged resistors more difficult. These
oscillations are removed from the data by randomizing the time information for
each track according to a uniform distribution $U(-T/2,T/2)$, where $T$ is the
time period of the oscillation. The modulations at $\omega_{\text{CBO}}$,
$\omega_{\text{2CBO}}$, $\omega_{\text{VW}}$, $\omega_{y}$, and
$\omega_{a}^{m}$ are all removed. This randomization procedure also has the
effect of removing the CBO decoherence phase-acceptance effect, which mostly
cancels in the sum of all calorimeter data as described above. Figure 11 shows
the smooth nature of the vertical projections of the mean $\langle y\rangle$
and width $\sqrt{\langle y^{2}\rangle}$ versus time after the randomization
procedure. Naively, one might assume that an average of the two tracker maps
will provide the most accurate representation of the muon distribution versus
time and decay coordinates. However, as indicated in Fig. 9, the early-to-late
behavior of the storage ring optical lattice functions varies significantly
versus azimuth, because the damaged resistors are located in only one section
of the storage ring. The correction is, therefore, evaluated independently
starting from each tracker station, and the final result is taken as their
average.

###  VIII.2 Estimation of beam distribution around the storage ring

The trackers measure $M^{T}(x,y,t)$ at two locations around the ring, but the
extraction of $\omega_{a}^{m}$ is performed using calorimeters at 24 azimuthal
locations. The damaged resistors were located on the opposite side of the ring
with respect to the two tracker positions. This asymmetry leads to differences
in the muon transverse distribution measured by the trackers and in the
proximity of calorimeters elsewhere in the ring. Both COSY and GM2RINGSIM
predict similar relative behavior of $\langle y\rangle(t)$ by azimuthal
location, as is shown with a comparison to calorimeter data in Fig. 12. But
the reduction of the vertical width (Fig. 11b)—which cannot be measured by the
calorimeters—is more important to the $C_{pa}$ calculation. The $M^{T}(x,y,t)$
distributions are evolved separately to all azimuthal locations using the
predictions from either COSY or GM2RINGSIM. This provides a complete spatial
and time distribution of the muons $M^{c}(x,y,t)$ in the vicinity of a
calorimeter based on a measurement at a single point in the ring. The time
dependent muon profiles were calculated at an azimuthal angle approximately
$22\text{\,}\mathrm{\SIUnitSymbolDegree}$ upstream of each calorimeter’s front
face, where acceptance is maximal.

Similarly to the treatment for the pitch correction (Sec. V), the vertical
distribution is anchored to the tracker data at its azimuthal location, and
the distribution at other parts of the storage ring is obtained by scaling the
width as a function of azimuth relative to that at the tracker station:

$y_{\text{pred}}(\phi,t)=y_{\text{tkr}}(t)\frac{y_{\text{pred}}^{\text{rms}}(\phi,t)}{y_{\text{pred}}^{\text{rms}}(\phi_{\text{tkr}},t)}.$
(24)

The width ratio is evaluated directly from virtual tracking planes in
GM2RINGSIM and via the ratio of beta functions in COSY
($\sqrt{\beta_{y}(\phi,t)/\beta_{y}(\phi_{\text{tkr}},t)}$). From such
transformations, vertical beam drifts are properly projected around the ring,
and further effects from permanent vertical closed-orbit distortions are
accounted for in the systematic errors of the correction.

The contribution to $C_{pa}$ from the radial motion of the beam is
subdominant, but for a complete treatment the radial beam changes are also
included. A similar procedure to that for the vertical changes is followed,
although there are added complications due to the interplay of muon momentum
and radial position as well as closed-orbit distortions due to azimuthal
variation of the main dipole field [10]. Each entry from the measured
distribution is scaled as

$x_{\text{pred}}(\phi,t)=\frac{x_{\text{pred}}^{\text{rms}}(\phi,t)}{x_{\text{pred}}^{\text{rms}}(\phi_{\text{tkr}},t)}\cdot(x_{\text{tkr}}(t)-\bar{x}_{\text{tkr}}(t))+\bar{x}_{\text{pred}}(\phi,t),$
(25)

where the predicted values for $x_{\text{rms}}(\phi,t)$ and $\bar{x}(\phi,t)$
are taken from tracking planes in GM2RINGSIM and calculated in COSY using
$\beta_{x}(\phi,t)$, $D_{x}(\phi,t)$, and the reconstructed distribution of
muon momenta. Measurements of the main dipole storage field are used to
calculate the closed-orbit distortion included in
$\bar{x}_{\text{pred}}(\phi,t)$. After including both vertical and radial
manipulations, the procedure has been shown to provide excellent agreement
between calculations of $C_{pa}$ using either extrapolated $M^{c}(x,y,t)$ or
beam-tracking planes in a closed-loop simulation test.

###  VIII.3 Phase-acceptance correction: Results and uncertainty evaluation

The muon distribution is combined with the simulated maps to extract the time-
dependent phase for each calorimeter $\varphi_{pa}^{c}(t)$ using the following
weighted sum:

$\varphi^{c}_{pa}(t)=\arctan\cfrac{\sum_{ij}M^{c}(x_{i},y_{j},t)\cdot\varepsilon^{c}(x_{i},y_{j})\cdot
A^{c}(x_{i},y_{j})\cdot\sin(\varphi^{c}_{pa}(x_{i},y_{j}))}{\sum_{ij}M^{c}(x_{i},y_{j},t)\cdot\varepsilon^{c}(x_{i},y_{j})\cdot
A^{c}(x_{i},y_{j})\cdot\cos(\varphi^{c}_{pa}(x_{i},y_{j}))}\,.$ (26)

The sum is over all spatial bins. The time dependence enters via the muon
distribution $M^{c}$. The acceptance, asymmetry, and phase maps for
calorimeter $c$ are represented by $\varepsilon^{c}$, $A^{c}$, and
$\varphi_{pa}^{c}$, respectively. An example evaluation of $\varphi^{c}_{pa}$
for the Run-1d dataset is shown in Fig. 22. As outlined earlier, the phase
increases with time primarily due to the decreasing vertical width of the
beam.

Figure 22: Calculation of $\varphi_{pa}^{c}(t)$ for calorimeter 19 in Run-1d
using data from the tracker station at
$180\text{\,}\mathrm{\SIUnitSymbolDegree}$ and Eq. 26. The fit start time for
the extraction of $\omega_{a}^{m}$ for this dataset is
$50\text{\,}\mathrm{\SIUnitSymbolMicro s}$ after injection in order to
mitigate the effect of $\varphi_{pa}^{c}(t)$.

The size of the mismeasurement of $\omega_{a}^{m}$ as a result of the time-
dependent phase $\varphi_{\text{pa}}^{\text{c}}(t)$ is estimated using
simulated data. A histogram is generated for each calorimeter according to the
full fit function used to extract $\omega_{a}^{m}$, an example of which is
described in Appendix D. The spectrum includes betatron oscillations of the
normalization, asymmetry, and phase and a modified $(g-2)$ phase term that has
a parameterization of $\varphi_{\text{pa}}^{\text{c}}(t)$. The function
parameters for each calorimeter are set according to their data-derived
values. The resulting simulated data are then fit with the full fit function
but excluding the modified phase term. The difference between the values of
$\omega_{a}^{m}$ that were input and extracted from the fit give the value of
the correction $C_{\text{pa}}$ for that calorimeter.

The final values correspond to the average of the corrections for all
calorimeters, obtained starting from measurements from each of the tracker
stations and from performing the procedure using both COSY and GM2RINGSIM.
Given that the phase increases with time as shown in Fig. 22, the measured
$\omega_{a}^{m}$ is larger than the true value. The correction $C_{\text{pa}}$
must reduce the measured value and is, therefore, negative. The $C_{pa}$
corrections are shown in Table 8. The values range from $-184$ ppb for Run-1a
to $-117$ ppb for Run-1c. The condition of the damaged resistors was worse for
the Run-1d dataset, which would generally lead to a larger $C_{pa}$ value. In
order to keep the correction for Run-1d comparable in size to the other
datasets, the fit start time was delayed from 30 to
$50\text{\,}\mathrm{\SIUnitSymbolMicro s}$, which reduces the correction by a
factor of 1.8.

Dataset | Run-1a | Run-1b | Run-1c | Run-1d
---|---|---|---|---
$C_{\text{pa}}$ | -184 | -165 | -117 | -164
Stat. uncertainty | 23 | 20 | 15 | 14
Tracker and CBO | 73 | 43 | 41 | 44
Phase maps | 52 | 49 | 35 | 46
Beam dynamics | 27 | 30 | 22 | 45
Total uncertainty | 96 | 74 | 60 | 80
Table 8: Values of the phase-acceptance correction $C_{\text{pa}}$ (ppb) and
their statistical, systematic, and total uncertainties for each of the Run-1
datasets.

The evaluations of statistical and systematic uncertainty are also detailed in
Table 8. The statistical uncertainty originates from the limited number of
tracks available to form $M^{T}(x,y,t)$ and ranges from 14 to 23 ppb depending
on the size of the dataset.

The sources of systematic uncertainty can be grouped into three main areas
with no single category dominating over the others. First, there is imperfect
knowledge of the tracker’s alignment, resolution, and acceptance corrections,
which all affect the measurement of $M^{T}(x,y,t)$. Uncertainties are
estimated by repeating the evaluation of $C_{\text{pa}}$ while varying the
corrections within a reasonable range based on external constraints.
Additionally, there is a statistical tension in Run-1a between the corrections
extracted independently from each of the two tracking stations. The
uncertainty for this dataset is conservatively inflated to account for this.
As explained previously, the phase-acceptance effect originating from CBO
decoherence will largely cancel when all the calorimeters are summed together.
The degree of this cancellation is dependent on the relative amplitudes of the
CBO-modulated acceptance oscillation around the ring, which is strongly
correlated with the tracker acceptance correction. The uncertainty from the
CBO decoherence cancellation is, therefore, combined linearly with those from
the tracker corrections, resulting in values ranging from 41 to 73 ppb.

Second, uncertainties associated with the estimation of the phase, asymmetry,
and acceptance maps in Eq. 26 are estimated using the GM2RINGSIM simulation.
These uncertainties are dominated by the knowledge of the phase map, which has
a much stronger effect on $C_{\text{pa}}$ than either the asymmetry or
acceptance maps. The possible variation of the map is estimated by comparing
measurements of the fitted phase in the actual $\omega_{a}^{m}$ fits with
predictions of $\varphi_{\text{pa}}^{\text{c}}$ from Eq. 26. The time
dependence of $M^{c}(x,y,t)$ in the latter is collapsed by taking data from
all times in a single time bin. The variation in the final fitted phase as a
function of calorimeter number and as a function of calorimeter $y$ position
matches the prediction to within 20%. This variation is used as an estimated
uncertainty on the phase map, and the procedure is repeated to propagate this
uncertainty to establish a contribution to $\delta C_{\text{pa}}$.
Interactions in material will also affect the detector acceptance and,
therefore, change the phase map (Fig. 29 in Appendix C). The effect of
material interactions in the simulation can be exaggerated to allow for
quantification of an uncertainty due to mismodeling of material effects. The
total uncertainty associated with the map estimation ranges from 46 to 52 ppb.

Last, the procedure utilizes beam dynamics models to extract the calorimeter-
specific $M^{c}(x,y,t)$ distributions from the tracker-measured
$M^{T}(x,y,t)$. Uncertainties are estimated by calculating $C_{\text{pa}}$
while varying the beta functions and magnetic fields within constraints from
measurements, varying the momentum distribution of the stored muons, and from
a comparison between values extracted using either COSY or GM2RINGSIM. This
uncertainty is smaller than those from the tracker and phase categories in
Runs-1a, 1b and 1c (from 22 to 27 ppb), but is comparable in size in Run-1d
($45$ ppb). The larger uncertainty in Run-1d is due to the worsening condition
of the damaged resistors, which increases the azimuthal variation and,
therefore, increases the reliance on the simulations.

Within each group of uncertainties, correlated effects are added linearly, and
then all remaining effects are summed quadratically to give the total. The
total uncertainties (from 60 to 96 ppb) range from 45% to 52% of the
correction value.

##  IX  Conclusion

In this paper, we described summaries and conclusions from in-depth studies of
four beam dynamics systematic corrections that are required to adjust the
measured muon precession frequency $\omega_{a}^{m}$ (Ref. [11]) to its true
physical value $\omega_{a}$. The corrections for the vertical betatron
motion—pitch—and the influence of the motional magnetic field on non-magic-
momentum muons—electric field—are well understood and have been documented by
previous generations of $g\\!-\\!2$ experiments. Here, we have refined these
studies and performed detailed uncertainty analyses. The pitch correction
requires knowledge of the vertical stored muon distribution from the in situ
tracker system, which provides detailed time-dependent stored-muon spatial
profiles in two areas of the storage ring. The electric field correction
requires knowledge of the stored muon momentum distribution (alternatively,
radial distribution), which is deduced from studying the time evolution of the
incoming muon bunch. Two distinct analysis methods are used.

Despite an initial scraping procedure meant to remove muons that will not
remain stored in the ring, a small fraction that will eventually strike a
fiducial-defining collimator do survive well into the measurement period. A
careful analysis demonstrated that those muons which are lost do not
significantly alter the ensemble-averaged spin phase, but we apply a small
data-driven correction.

Finally, owing to 2 of 32 damaged high-voltage resistors in the ESQ system
that led to slower-than-designed charging times for two plates, the mean and
rms of the stored muon distribution in Run-1 evolved throughout the first
${\sim}100$ $\text{\,}\mathrm{\SIUnitSymbolMicro s}$ of the measuring period.
Investigation of this effect led to an extensive and new understanding of a
subtle coupling of decay coordinate within the storage volume and acceptance
versus spin phase. The coupling to acceptance and phase is carried out by
detailed GEANT4-based simulations, and the understanding of the beam motion
around the ring through the evolution of the beta functions is simulated in
detail using COSY, whose key results are supported by storage simulations
using GM2RINGSIM. This “phase-acceptance” effect is time-dependent; its
influence on the phase changes rapidly and faster than the muon population
decays. Therefore, in the Run-1d period where the resistor time constants were
longer than during periods 1a, 1b, and 1c, a fit start time in the precession
fits of $50\text{\,}\mathrm{\SIUnitSymbolMicro s}$ is used to limit the
systematic shift to the phase.

###  IX.1 Summary of Run-1 net corrections to $\omega_{a}^{m}$

Table 9 presents a summary overview of key findings discussed in this paper.
The statistical uncertainty from the asymmetry-weighted analyses of the muon
precession data is provided for reference. The four corrections discussed
herein are multiplicative adjustments that are applied to convert the raw
extracted frequency $\omega_{a}^{m}$ from the precession fits to the true
$\omega_{a}$ value that is needed to determine $a_{\mu}$ as in Eq. 1. The
corrections are identified as electric field $C_{e}$, Sec. IV; pitch $C_{p}$,
Sec. V; muon-loss-phase $C_{ml}$, Sec. VII; and phase-acceptance $C_{pa}$,
Sec. VIII. Accordingly,

$a_{\mu}=\left(\frac{m}{q\tilde{B}}\right)\left[\omega_{a}^{m}\cdot(1+C_{e}+C_{p}+C_{ml}+C_{pa})\right],$
(27)

where $\tilde{B}$ represents the muon-weighted average magnetic field that is
discussed in Ref. [10].

| Correction (ppb) | Uncertainty (ppb)
---|---|---
$\omega_{a}^{m}$ statistical | $\cdot\cdot\cdot$ | 434
$C_{e}$ | 489 | 53
$C_{p}$ | 180 | 13
$C_{ml}$ | $-11$ | 5
$C_{pa}$ | $-158$ | 75
$C_{\text{total}}$ | 499 | 93
Table 9: The Run-1 combined beam dynamics corrections to $\omega_{a}^{m}$ from
the four Run-1 datasets. These values are computed with the full correlation
matrix formalism used to provide the measured value for $\omega_{a}$. The
values, thus, reflect fully correlated systematics and weighting from the
statistical uncertainties from each dataset.

##  X Acknowledgments

We thank the Fermilab management and staff for their strong support of this
experiment, as well as the tremendous support from our university and national
laboratory engineers, technicians, and workshops. The Muon $g\\!-\\!2$
Experiment was performed at the Fermi National Accelerator Laboratory, 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. Additional support for the experiment was provided by the
Department of Energy offices of High Energy Physics and Nuclear Physics (USA),
the National Science Foundation (USA), the Istituto Nazionale di Fisica
Nucleare (Italy), the Science and Technology Facilities Council (UK), the
Royal Society (UK), the European Union’s Horizon 2020 research and innovation
programme under the Marie Skłodowska-Curie Grant Agreements No. 690835 and No.
734303, the National Natural Science Foundation of China (Grants No. 11975153
and No. 12075151), MSIP, NRF and IBS-R017-D1 (Republic of Korea), and the
German Research Foundation (DFG) through the Cluster of Excellence PRISMA+
(EXC 2118/1, Project ID No. 39083149).

##  Appendix A Electric Quadrupole Nonlinearity

Consider a 2D projection of the electric field in the ESQ region, as shown in
Fig. 23. Laplace’s equation in two dimensions and Cartesian coordinates is

$\nabla^{2}V(x,y)=\frac{\partial^{2}V}{\partial
x^{2}}+\frac{\partial^{2}V}{\partial y^{2}}=0.$ (28)

A potential corresponding to a pure quadrupole field can be written as
$V=\frac{1}{2}k(x^{2}-y^{2})$ for some constant $k$. The electric field is
$E={\vec{\nabla}}V=k(x\hat{x}-y\hat{y})$ with zero divergence, and it is
linear in both $x$ and $y$. Higher-order terms may appear owing to the plate
geometry, but the symmetry in the 2D Cartesian limit permits only those terms
that are odd in $x$, $y$. In particular there is no sextupole (${\sim}kx^{2}$)
dependence of the field that would be symmetric in displacement.

In the limit of finite curvature, as employed in the storage ring geometry, it
is more appropriate to represent the fields in cylindrical coordinates, where
Laplace’s equation is

$\nabla^{2}V(\rho,y,\phi)=\left[\frac{1}{\rho}\frac{\partial}{\partial\rho}\left(\rho\frac{\partial}{\partial\rho}\right)+\frac{1}{\rho^{2}}\frac{\partial^{2}}{\partial\phi^{2}}+\frac{\partial^{2}}{\partial
y^{2}}\right]V=0$ (29)

with $\rho=R_{0}+x$. Using the assumption—which is not accurate—that the ESQ
plates are continuous around the ring so that there is no dependence on the
angular coordinate $\phi$, the simplest possible solution is

$V(\rho,y)=\tilde{k}\left[\frac{1}{2}\left(\frac{\rho^{2}}{R_{0}^{2}}-1\right)-\ln\frac{\rho}{R_{0}}-\left(\frac{y}{R_{0}}\right)^{2}\right],$
(30)

where $\tilde{k}$ is a constant. The electric field is given by

$\vec{E}=-\vec{\nabla}V=-\frac{\tilde{k}}{R_{0}^{2}}\left[\left(\rho-\frac{R_{0}^{2}}{\rho}\right)\hat{\rho}-2y\hat{y}\right]$
(31)

and expanding in $x/R_{0}$:

$\vec{E}\approx
k\left[\left(x-\frac{x^{2}}{2R_{0}}+\ldots\right){\hat{\rho}}-y\hat{y}\right],$
(32)

where $k=-2\tilde{k}/R_{0}^{2}$. Evidently, Maxwell’s equations require a term
quadratic in displacement. The quadratic term is equivalent to a sextupolelike
component that will affect the chromaticity, complicate the $E$-field
correction, and possibly drive a third-order resonance. The general conclusion
is that there necessarily exists a symmetric component to the $E$-field that
will contribute to the evaluation of the correction and its systematic
uncertainty. This uncertainty is included in Table 5 as part of the ESQ
calibration.

Figure 23: A 2D representation of the equipotentials and electric field in the
ESQ region. The circle at radius 45 mm represents the storage ring aperture
defined by the collimators.

###  A.1 Path length correction owing to curved geometry

The radial electric field along the trajectory of a muon with equilibrium
radial offset $x_{e}$ and betatron amplitude $x_{\beta}$ is

$E_{r}(s)=k(\eta\delta+x_{\beta})-\frac{1}{2R_{0}}(\eta\delta+x_{\beta})^{2}+\ldots,$
(33)

where $x_{e}=\eta\delta=\eta\frac{\Delta p}{p}$ and the simple solution above
to Laplace’s equation is assumed. The average electric field along the
trajectory is

$\langle
E_{r}(s)\rangle=\frac{k}{L}\int_{0}\left(\eta\delta+x_{\beta}-\frac{1}{2R_{0}}(\eta\delta+x_{\beta})^{2}\right)(1+x_{\beta}/R_{0}){ds},$
(34)

where $L$ is the length of the trajectory. Using $x_{\beta}=x_{\beta
0}\cos\phi(s)$ with $d\phi=\frac{ds}{R_{0}}=\frac{dl}{(R_{0}+x)}$, it can be
shown that $\langle E_{r}\rangle$ is approximately

$\langle E_{r}(s)\rangle\sim
k\left(\eta\delta-\frac{1}{2R_{0}}((\eta\delta)^{2}-\frac{1}{2}x_{\beta
0}^{2})+\ldots\right),$ (35)

where $k=n\left(\frac{v_{s}B}{R_{0}}\right)$. One observes that the
contribution to the average electric field along a trajectory, due to the
curvature of the ESQ plates, scales quadratically with the betatron amplitude.

Substitution of Eq. 35 into Eq. 13 gives the correction to $\omega_{a}$ due to
fractional momentum offset $\delta$ and betatron amplitude $x_{\beta 0}$ as

$\displaystyle C_{e}(\delta,x_{\beta 0})$ $\displaystyle\approx$
$\displaystyle-2\delta\frac{\beta
k}{cB}\left(\eta\delta-\frac{1}{2R_{0}}((\eta\delta)^{2}-\frac{1}{2}x_{\beta
0}^{2})\right)$ $\displaystyle\approx$ $\displaystyle-2\frac{\beta
k}{cB}\left(\eta\delta^{2}-\frac{1}{2R_{0}}(\eta^{2}\delta^{3}-\frac{1}{2}x_{\beta
0}^{2}\delta)\right).$

The next step requires averaging $C_{e}(\delta,x_{\beta 0})$ over the entire
momentum and CBO distribution:

$\displaystyle\langle{C_{e}(x_{e},x_{\beta 0})}\rangle\approx$
$\displaystyle-2\frac{\beta k}{cB}\Bigg{[}\frac{\langle
x_{e}^{2}\rangle}{\eta}$ (36)
$\displaystyle-\frac{1}{2R_{0}}\left(\frac{\langle
x_{e}^{3}\rangle}{\eta}-\frac{1}{2}\langle{x_{\beta
0}^{2}}\rangle\frac{{\langle{x_{e}}\rangle}}{\eta}\right)\Bigg{]},$

where ${\langle{\delta}\rangle}={\langle{x_{e}}\rangle}/\eta$. If
${\langle{x_{e}}\rangle}/\eta=0$, and the betatron amplitude and momentum
offset are assumed to be uncorrelated, and there is no contribution from the
sextupole or path length terms, then

$\langle C_{e}(\delta,x_{\beta 0})\rangle=-2\frac{\beta
k}{cB}\frac{{\langle{\delta^{2}}\rangle}}{\eta}.$ (37)

In fact, momentum offset and betatron amplitude are strongly anticorrelated,
which will be discussed below. If $\langle x_{e}\rangle\neq 0$, that is, the
muon momentum distribution mean is not centered at the magic momentum, then
the fractional change to the $E$-field correction from the betatron amplitude
consideration is

$\displaystyle\frac{{\langle{\Delta
C_{e}(x_{e},x_{\beta_{0}}}\rangle}}{{\langle{C_{e}}\rangle}}$
$\displaystyle\sim$
$\displaystyle\frac{-\frac{1}{2R_{0}}\left[\frac{{x_{e}^{3}}}{\eta}-{\frac{1}{2}}\frac{{\langle{x^{2}_{\beta_{0}}}\rangle}{\langle{x_{e}}\rangle}}{\eta}\right]}{\frac{x_{e}^{2}}{\eta}}$
$\displaystyle\sim$
$\displaystyle\frac{-{\langle{x_{e}^{3}}\rangle}+{\langle{x_{e}}\rangle}{\langle{x^{2}_{\beta
0}}\rangle}/2}{2{\langle{x_{e}^{2}}\rangle}R_{0}}.$

For the equilibrium radial distributions typical of Run-1 (see Fig. 6b),

$\frac{{\langle{\Delta C_{e}}\rangle}}{{\langle{C_{e}}\rangle}}\sim-1.1\times
10^{-3}+(1.9\times 10^{-6}{\rm mm^{-2}}){\langle{x_{\beta_{0}}^{2}}\rangle}.$

Conservatively estimating that ${\langle{x_{\beta_{0}}^{2}}\rangle}<(45{\rm
mm)^{2}}$, the maximum radius of the storage volume, we find that the change
in the $E$-field correction, due to the nonlinearity associated with the
curvature of the ESQ plates, is less than 0.5%.

###  A.2 Effect from the intrinsic quadrupole nonlinearity

The rectangular cross section of the ESQ geometry introduces a nonlinearity
compared to the ideal case. The electric field along the horizontal axis
($y=0$) can be expressed by

$E_{x}-iE_{y}=(b_{n}-ia_{n})\frac{x^{n}}{r_{0}^{n}},$ (38)

where $r_{0}=0.045$ m and $a_{n}$ and $b_{n}$ are coefficients of the
multipoles, which have been fit to an azimuthal slice of a OPERA-3D [27] field
map. The fit is for a horizontally pure basis of McMillan functions [43].
Figure 24 shows the horizontal electric field in the midplane, with the
evident nonlinear behavior at the extremes. The values from the OPERA map and
from the multipole expansion are superimposed and are evidently in excellent
agreement.

The fitted sextupolelike coefficient can be compared to the estimate discussed
previously. A solution satisfying the Laplacian in the curved system was given
in Eq. 32. The ratio of the coefficients of the quadratic and linear terms is
$r_{hyp}=-\frac{1}{2R_{0}}=-0.0703{\rm m}^{-1}$, but that hypothetical
solution is not unique. Although it satisfies Maxwell’s equations, boundary
conditions were not yet imposed. The ratio based on the fit to the OPERA field
map, that satisfies both Maxwell and the boundary conditions, is
$r_{fit}=\frac{b_{2}}{b_{1}r_{0}}=-2.71281\times 10^{3}/(1.01609\times
10^{6}\times 0.045{\rm m})=-0.0593{\rm m}^{-1}$, within 16% of our guess.

As was shown above, the effect of the sextupolelike component of the ESQ
field, and the asymmetry of the path length about the magic radius, is that
the E-field correction depends on betatron amplitude and nonlinearly on
equilibrium radial offset. The sextupolelike component and the path-length
component contribute with opposite sign, and the amplitude of the
sextupolelike component is about 1/2 of the path-length piece. Based on the
measured equilibrium radial distribution, the contribution of the
sextupolelike component to the $E$-field correction is less than 1%. The
potential effect from higher-order multipoles was explored in simulation.

Figure 24: The electric field along the $x$ axis in the midplane ($y=0$). The
purple points are computed from the multipole expansion using Eq. 38. The
black curve uses the values from the field map.

##  Appendix B $E$-field and Pitch Corrections Explored Numerically

In simulation, $\omega_{a}$ can be determined directly by spin tracking. The
trajectory of a muon is established by numerically integrating the equations
of motion, and its spin by integration of the Thomas-BMT equation along that
trajectory. In this section, the contribution to $\omega_{a}$ due to the
electric field and the pitch correction (vertical betatron motion) is explored
numerically, in particular, to study the nonlinearity that arises from the ESQ
geometry, mechanical alignment precision, and high-voltage uncertainty. Spin
tracking naturally incorporates the $E$-field and pitch effects; they cannot
be trivially separated.

Full spin tracking verifies to high precision the conclusions for the
$E$-field and pitch corrections as described in Eqs. 16 and 17, respectively.
It is, however, computationally intensive. Furthermore, for a trajectory that
is both off momentum and oscillating vertically, spin tracking cannot
distinguish the contributions from the electric field and pitch. In order to
explore the full landscape of possible systematic uncertainties and to
understand the contributions from pitch and $E$-field independently, an
efficient and equivalently precise method was developed.

The anomalous precession frequency is commonly approximated as

$\displaystyle\vec{\omega}_{a}=-\frac{q}{m}\Bigg{[}a_{\mu}\vec{B}-\Big{(}a_{\mu}$
$\displaystyle-\frac{1}{\gamma^{2}-1}\Big{)}\frac{\vec{\beta}\times\vec{E}}{c}$
(39)
$\displaystyle\left.-\,a_{\mu}\left(\frac{\gamma}{\gamma+1}\right)(\vec{\beta}\cdot\vec{B})\vec{\beta}\right],$

where the experiment actually measures the scalar quantity $\omega_{a}^{m}$
and not its vector approximation. However, Eq. 39 can be used to produce spin
dynamics simulations with short run times that allow for the study of a large
number of configurations. First, consider the $\vec{\beta}\times\vec{E}$
electric field term. Define $\Phi_{a}=\int_{0}^{T}\omega_{a}dt$, there $T$ is
the total integrated time that the muon is precessing. In the absence of
$\vec{E}$, $\Phi_{a}$ can be computed by integrating $\vec{B}_{\perp}$ along
the muon trajectory. The contribution from the electric field to $\phi_{a}$ is
given by

$\begin{split}\frac{\Delta\Phi_{E}}{\Phi_{a}}&=\left(a_{\mu}-\frac{m^{2}c^{2}}{p^{2}}\right)\frac{1}{c}\int_{0}^{T}\vec{\beta}\times\vec{E}dt\\\
&\approx-2\frac{\Delta
p}{p}\frac{1}{c}\int_{0}^{T}\vec{\beta}\times\vec{E}dt,\end{split}$ (40)

where $\Delta p=p-p_{0}$. In the tracking code, it is straightforward to
compute the sums in Eq. 40. The electric field correction for the trajectory,
assuming that $\vec{B}$ is parallel to $\vec{\beta}\times\vec{E}$, is

$C_{e}(T)=-2\frac{\Delta
p}{p}\frac{1}{T}\int^{T}\frac{\vec{\beta}\times\vec{E}}{Bc}dt.$ (41)

In a similar manner, the expressions for the pitch correction $C_{p}$
developed in Sec. V can be rearranged to read

$C_{p}(T)=\frac{\Delta\Phi_{p}}{\Phi_{a}}=\frac{1}{T}\int^{T}(1-|{\hat{\bm{\beta}}}\times\hat{\bf
B}|)dt.$ (42)

To test the equivalence of the shift in precession frequency based on Eqs. 41
and 42, both were computed simultaneously along a single trajectory. The muon
is propagated for 1000 $2\pi/\omega_{a}$ in order to determine the spin tune
by integrating the Thomas-BMT equation with sufficient accuracy, and
$\vec{\beta}\times\vec{E}$ is summed along that same trajectory. The
trajectories chosen to test Eq. 41 have momentum offsets but no vertical
motion. To test Eq. 42, the tracked muon has magic momentum and zero radial
betatron oscillation. The agreement is excellent, as shown in Fig. 25, which
shows $C_{e}$ and $C_{p}$ for both spin tracking and the integration method.

Figure 25: Top: contribution to $\omega_{a}$ due to the electric field as
computed by spin tracking (integration of the Thomas-BMT equation) and
integration along trajectories (Eq. 41) for closed orbits at varying radial
offsets. The vertical closed orbit is zero and there is no vertical betatron
motion. A representative, measured radial distribution is superimposed for
reference. Bottom: contribution to $\omega_{a}$ due to the vertical
oscillation computed with spin tracking and integration along the trajectory
(Eq. 42) for muons with magic momentum and zero radial betatron amplitude.
Note that the formula that assumes linearity overestimates the contribution at
large amplitudes.

In the tests of systematics that are described next, the model of the
quadrupole field is based on the sum of OPERA-based field maps for each flat
ESQ plate and, therefore, includes the nonlinearity associated with the
geometry, as is evident in Fig. 24. At large amplitudes, the formulas that
assume linearity of the electric field slightly overestimate the contributions
of the electric field and pitch to $\omega_{a}$.

The as-built ESQ plate alignment is known with a precision better than $\pm 1$
mm, and the high-voltage set points to $\pm 0.01\%$. To explore the
uncertainty from unknown misalignments or voltage errors, 1024 configurations
were analyzed, varying the conditions within the specifications. In each case,
the $E$-field and pitch corrections were computed and compared to the
reference values. From these studies, a systematic uncertainty of
$\sigma(C_{p})=1.3$ ppb is assigned for the contribution from pitching motion
due to misalignment and voltage errors. The distribution of perturbations to
the electric field correction is not Gaussianm and the standard deviation may
underestimate the uncertainty. Indeed, the error distribution has a strong
dependence on the nature of the simulated misalignment and field-error
configuration. The systematic uncertainty in the measurement of the
contribution from the electric field due to field errors and misalignment is
chosen conservatively to be the full width of distribution simulated values,
yielding $\sigma(C_{e})=8.7$ ppb.

##  Appendix C Maps for the Phase-Acceptance Correction

Estimation of the phase-acceptance correction relies on high-precision maps of
detection acceptance, asymmetry, and phase as a function of muon decay
coordinate $(x,y,\phi)$. In this section, we describe the procedure by which
these maps are produced using our GM2RINGSIM GEANT4-based simulation in which
all detector systems (calorimeters and trackers) and vacuum chamber hardware
of the experiment are fully modeled. The simulation first determines the
random time of a muon decay drawn from an exponential distribution. The muon
four-momentum, position, and polarization direction are calculated based on
the time spent in the storage ring. Finally, the muon is placed at the
precalculated position and allowed to decay immediately. The orientation of
the muon spin at the time of decay is determined from propagating the
precession at the $\omega_{a}$ frequency, which is known to sufficient
precision from the E821 experiment.

Figure 26: GM2RINGSIM simulation of detected positron energy as a function of
time, with muon decay disabled.

The simulated time and energy distributions of the decay positron match well
with the experimental measurement as shown in Fig. 26. This 2D intensity
map—with the decaying exponential removed—shows positron energy versus time in
fill from 30 to 80 $\text{\,}\mathrm{\SIUnitSymbolMicro s}$. The clear
modulations are at the frequency $\omega_{a}$. The asymmetry reverses around 1
GeV, as expected from the kinematic boost of the Michel spectrum. The
simulation does not include CBO modulations because we are only interested in
the average detection acceptance effect; the beam dynamic effects that are
oscillating at a different frequency are strongly suppressed. Roughly 9
billion muon-decay events are simulated and analyzed to produce maps for three
$\omega_{a}$ analysis methods with different positron weightings (see Ref.
[11]). In practice, the asymmetry-weighted method provides the most precise
determination of $\omega_{a}$, and, therefore, we will show its maps.

The detection acceptance, asymmetry, and phase maps are produced in the form
of 3D histograms $H^{c}(x,y,E)$, where $x$ and $y$ are the radial and vertical
coordinates, respectively, in the storage ring volume, $E$ is the decay
positron energy, and $c$ is a calorimeter index from 1 to 24. For simplicity,
we can express these maps as $H^{c}_{ijk}$, where $i$ and $j$ are $x$ and $y$
coordinates in 5 mm bins, and $k$ is a 50 MeV $E$ bin.

###  C.1 Acceptance maps

The average acceptance $\bar{\varepsilon}_{ij}$ for a muon decaying in
transverse position bin $ij$ is given by

$\bar{\varepsilon}_{ij}=\sum_{k}f_{k}\varepsilon_{ijk}\sum_{k^{\prime}>k_{th}^{\prime}}R_{ijkk^{\prime}}w_{k^{\prime}},$
(43)

where the first sum is over positron truth energy bins $k$, the second sum is
over positron deposited energy bins $k^{\prime}$ above a threshold
$k_{th}^{\prime}$, $f_{k}$ is the fraction of events with energy $k$,
$\varepsilon_{ijk}$ is the geometric acceptance for muon decays in bin $ij$
and positron energy $k$, $R_{ijkk^{\prime}}$ is the calorimeter response
function, and $w_{k^{\prime}}$ is the analysis weighting factor [11].
$R_{ijkk^{\prime}}$ represents the probability that a positron in position bin
$ij$ and truth energy bin $k$ ends up in the $k^{\prime}$ detected energy bin.
This response function is illustrated in Fig. 27.

Figure 27: The average calorimeter response function showing the deposited
energy versus the truth energy, integrated over the decay volume.

Figure 28 shows the relative acceptance function maps for each calorimeter,
weighted for the asymmetry-method. The variations among calorimeters are
attributed to the material upstream of each, as is shown in Fig. 29. We can
loosely identify five categories: i) behind ESQ plates; ii) behind kicker
plates; iii) downstream of the trackers; iv) next to the larger inflector
vacuum chamber; and v) nothing in front.

Figure 28: The asymmetry-weighted relative acceptance maps for all
calorimeters. The differences in acceptance derive from the materials just
upstream of each station; see Fig. 29.

###  C.2 Phase and Asymmetry maps

To produce the asymmetry and phase maps, two types of time spectra are
generated: first, a time spectrum for all decay positrons (no matter whether
the $e^{+}$ is detected or not) and, second, a decay $x-y$ grid of time
spectra for each detected positron energy bin. Each time spectrum in the $ijk$
bin is fitted with the five-parameter function
$N(t)=N_{0}e^{-t/\tau}(1+A\cos{(\omega_{a}t+\varphi)}),$ where $\varphi$ is
the fitted phase. The fitted phase for all decay positrons is denoted as
$\varphi_{0}$, and the fitted phase for each decay $xy$ bin and energy bin is
denoted as $\varphi_{ijk}$. The value of each $ijk$ bin in the phase map is
then given by $\varphi_{ijk}-\varphi_{0}$ Similarly, the value of each $ijk$
bin in the asymmetry map is given by $A_{ijk}$.

Figure 29 shows beam-weighted vertical projections of the extracted phase
maps. Again, the variation from detector to detector shows a strong dependence
on the material in front of the calorimeter. The vertical width changes of the
beam couple with the variation in these distributions such that some
calorimeters, such as those behind ESQ plates, experience a much larger phase-
acceptance effect than those behind the tracker stations.

Figure 29: The asymmetry-weighted vertical projections of the phase maps by
calorimeter. The calorimeters are ordered sequentially in azimuth with respect
to the injection location. The differences between stations is caused by
material differences that impact the transmission of positron decays enroute
to detectors. They also impact the acceptance (see Fig. 28) and the asymmetry
maps. These indices represent: Q, behind ESQ plates; K, behind kicker plates;
T, behind tracker stations; IV, in the constrained inflector vacuum chamber;
and, F, free space in front.

##  Appendix D Precession Fit Function

The complete fit function for the positron arrival time distribution described
in Ref. [11] is

$\displaystyle F(t)=N_{0}\cdot N_{x}(t)\cdot N_{y}(t)\cdot\Lambda(t)\cdot
e^{-t/\gamma\tau_{\mu}}\cdot$ $\displaystyle\left[1+A_{0}\cdot
A_{x}(t)\cdot\cos{(\omega_{a}^{m}t+\varphi_{0}\cdot\varphi_{x}(t))}\right].$
(44)

Here, $\Lambda(t)$ is the muon loss function described in Sec. VII.1. The
terms $N_{x}$, $N_{y}$, $A_{x}$ and $\varphi_{x}$ describe the interplay
between acceptance and beam dynamics that affect the overall calorimeter rate,
the average asymmetry, and the average phase. They are defined as

$\displaystyle N_{x}(t)={}1$
$\displaystyle{}+{}e^{-1t/\tau_{\text{CBO}}}A_{N,x,1,1}\cos(1\omega_{\text{CBO}}t+\varphi_{N,x,1,1})$
$\displaystyle{}+{}e^{-2t/\tau_{\text{CBO}}}A_{N,x,2,2}\cos(2\omega_{\text{CBO}}t+\varphi_{N,x,2,2}),$
(45) $\displaystyle N_{y}(t)={}1$
$\displaystyle{}+{}e^{-1t/\tau_{y}}A_{N,y,1,1}\cos(1\omega_{y}t+\varphi_{N,y,1,1})$
$\displaystyle{}+{}e^{-2t/\tau_{y}}A_{N,y,2,2}\cos(1\omega_{\text{VW}}t+\varphi_{N,y,2,2}),$
(46) $\displaystyle A_{x}(t)={}1$
$\displaystyle{}+{}e^{-1t/\tau_{\text{CBO}}}A_{A,x,1,1}\cos(1\omega_{\text{CBO}}t+\varphi_{A,x,1,1}),$
(47) $\displaystyle\varphi_{x}(t)={}1$
$\displaystyle{}+{}e^{-1t/\tau_{\text{CBO}}}A_{\varphi,x,1,1}\cos(1\omega_{\text{CBO}}t+\varphi_{\varphi,x,1,1}).$
(48)

The parameters of the form $A_{N,x,i,j}$ and $\varphi_{N,x,i,j}$ correspond to
the effect of the $i$th moment of the radial ($x$) beam distribution at the
$j$th multiple of the relevant fundamental frequency (in this case,
$\omega_{\rm CBO}$) on the rate normalization $N$. Analogous parameters model
the modulation of the average asymmetry $A$ and phase $\varphi$, as well as
the effect of moments of the vertical ($y$) beam distribution.

## References

  * Bailey  _et al._ [1972] J. Bailey, W. Bartl, G. von Bochmann, R. Brown, F. Farley, M. Giesch, H. Jostlein, S. van der Meer, E. Picasso, and R. Williams, Precise measurement of the anomalous magnetic moment of the muon, Nuovo Cimento A 9, 369 (1972).
  * Bailey  _et al._ [1979] J. Bailey _et al._ (CERN-Mainz-Daresbury Collaboration), Final report on the CERN muon storage ring including the anomalous magnetic moment and the electric dipole moment of the muon, and a direct test of relativistic time dilation, Nucl. Phys. B150, 1 (1979).
  * Bennett  _et al._ [2006] G. Bennett _et al._ (Muon $g-2$ Collaboration), Final report of the Muon E821 anomalous magnetic moment measurement at BNL, Phys. Rev. D 73, 072003 (2006).
  * Aoyama  _et al._ [2020] T. Aoyama _et al._ , The anomalous magnetic moment of the muon in the Standard Model, Phys. Rep. 887, 1 (2020).
  * Phillips  _et al._ [1977] W. D. Phillips, W. E. Cooke, and D. Kleppner, Magnetic moment of the proton in H2O in Bohr magnetons, Metrologia 13, 179 (1977).
  * [6] E. Tiesinga, P. J. Mohr, D. B. Newell, and B. N. Taylor, The 2018 CODATA recommended values of the fundamental physical constants (Web version 8.1), https://physics.nist.gov/cuu/Constants.
  * Liu  _et al._ [1999] W. Liu _et al._ , High Precision Measurements of the Ground State Hyperfine Structure Interval of Muonium and of the Muon Magnetic Moment, Phys. Rev. Lett. 82, 711 (1999).
  * Hanneke  _et al._ [2011] D. Hanneke, S. F. Hoogerheide, and G. Gabrielse, Cavity control of a single-electron quantum cyclotron: Measuring the electron magnetic moment, Phys. Rev. A 83, 052122 (2011).
  * Abi  _et al._ [2021] B. Abi _et al._ (Muon $g-2$ Collaboration), companion Letter, Measurement of the Positive Muon Anomalous Magnetic Moment to 0.46 ppm, Phys. Rev. Lett. 126, 141801 (2021).
  * Albahri  _et al._ [2021a] T. Albahri _et al._ (Muon $g-2$ Collaboration), companion paper, Magnetic field analysis for the Muon $g-2$ Experiment at Fermilab, Phys. Rev. A 103, 042208 (2021a).
  * Albahri  _et al._ [2021b] T. Albahri _et al._ (Muon $g-2$ Collaboration), companion paper, Measurement of the anomalous precession frequency of the muon in the Fermilab Muon $g-2$ Experiment, Phys. Rev. D 103, 072002 (2021b).
  * Bargmann  _et al._ [1959] V. Bargmann, L. Michel, and V. Telegdi, Precession of the Polarization of Particles Moving in a Homogeneous Electromagnetic Field, Phys. Rev. Lett. 2, 435 (1959).
  * Jackson [1999] J. D. Jackson,  _Classical Electrodynamics_ , 3rd ed. (Wiley, New York, 1999).
  * Bennett  _et al._ [2009] G. Bennett _et al._ (Muon $g-2$ Collaboration), An improved limit on the muon electric dipole moment, Phys. Rev. D 80, 052008 (2009).
  * Stratakis  _et al._ [2019] D. Stratakis, B. Drendel, J. P. Morgan, M. J. Syphers, and N. S. Froemming, Commissioning and first results of the Fermilab Muon Campus, Phys. Rev. Accel. Beams 22, 011001 (2019).
  * Danby  _et al._ [2001] G. T. Danby _et al._ , The Brookhaven muon storage ring magnet, Nucl. Instrum. Methods Phys. Res., Sect. A 457, 151 (2001).
  * Prigl  _et al._ [1996] R. Prigl, U. Haeberlen, K. Jungmann, G. zu Putlitz, and P. von Walter, A high precision magnetometer based on pulsed NMR, Nucl. Instrum. Methods Phys. Res., Sect. A 374, 118 (1996).
  * Yamamoto  _et al._ [2002] A. Yamamoto _et al._ , The superconducting inflector for the BNL $g-2$ experiment, Nucl. Instrum. Methods Phys. Res., Sect. A 491, 23 (2002).
  * Semertzidis  _et al._ [2003] Y. K. Semertzidis _et al._ , The Brookhaven muon ($g-2$) storage ring high voltage quadrupoles, Nucl. Instrum. Methods Phys. Res., Sect. A 503, 458 (2003).
  * Khaw  _et al._ [2019] K. S. Khaw _et al._ (Muon $g-2$ Collaboration), Performance of the Muon $g-2$ calorimeter and readout systems measured with test beam data, Nucl. Instrum. Methods Phys. Res., Sect. A 945, 162558 (2019).
  * Kaspar  _et al._ [2017] J. Kaspar _et al._ , Design and performance of SiPM-based readout of PbF2 crystals for high-rate, precision timing applications, J. Instrum. 12, P01009.
  * Anastasi  _et al._ [2019] A. Anastasi _et al._ (Muon $g-2$ Collaboration), The laser-based gain monitoring system of the calorimeters in the Muon $g-2$ experiment at Fermilab, J. Instrum. 14, P11025.
  * Metodiev  _et al._ [2015] E. M. Metodiev _et al._ , Analytical benchmarks for precision particle tracking in electric and magnetic rings, Nucl. Instrum. Methods Phys. Res., Sect. A 797, 311 (2015).
  * Arvanitis and Lyon [2014] T. Arvanitis and A. Lyon, artG4: A generic framework for GEANT4 simulations, J. Phys. Conf. Ser. 513, 022023 (2014).
  * Agostinelli  _et al._ [2003] S. Agostinelli _et al._ (GEANT4 Collaboration), GEANT4: A simulation toolkit, Nucl. Instrum. Methods Phys. Res., Sect. A 506, 250 (2003).
  * Poole  _et al._ [2012] C. M. Poole, I. Cornelius, J. V. Trapp, and C. M. Langton, A CAD interface for GEANT4, Australasian Physical and Engineering Sciences in Medicine 35, 329 (2012).
  * Ope [2004] _OPERA-3D User Guide_ (Vector Fields Limited, Oxford, England, 2004).
  * Tarazona  _et al._ [2019] D. A. Tarazona, M. Berz, and K. Makino, Muon loss rates from betatron resonances at the Muon $g-2$ storage ring at Fermilab, Int. J. Mod. Phys. A 34, 1942008 (2019).
  * Makino and Berz [2006] K. Makino and M. Berz, COSY INFINITY version 9, Nucl. Instrum. Methods Phys. Res., Sect. A 558, 346 (2006).
  * Valetov  _et al._ [2019] E. Valetov, M. Berz, and K. Makino, Computation of the main and fringe fields for the electrostatic quadrupoles of the Muon $g-2$ storage ring, Int. J. Mod. Phys. A 34, 1942041 (2019).
  * [31] INTEGRATED Engineering Software, COULOMB, https://www.integratedsoft.com/Products/coulomb.aspx.
  * Berz [1999] M. Berz,  _Modern Map Methods in Particle Beam Physics_ (Academic Press, New York, 1999).
  * Sagan [2006] D. Sagan, BMAD: A relativistic charged particle simulation library, Nucl. Instrum. Methods Phys. Res., Sect. A 558, 356 (2006).
  * Orlov  _et al._ [2002] Y. Orlov, C. S. Ozben, and Y. Semertzidis, Muon revolution frequency distribution from a partial-time Fourier transform of the $g-2$ signal in the Muon $g-2$ Experiment, Nucl. Instrum. Methods Phys. Res., Sect. A 482, 767 (2002).
  * Combley and Picasso [1974] F. Combley and E. Picasso, The muon ($g-2$) precession experiments: Past, present and future, Phys. Rep. 14, 1 (1974).
  * Farley and Semertzidis [2004] F. J. M. Farley and Y. K. Semertzidis, The 47 years of muon $g-2$, Prog. Part. Nucl. Phys. 52, 1 (2004).
  * Miller  _et al._ [2012] J. P. Miller, E. de Rafael, B. L. Roberts, and D. Stöckinger, Muon ($g-2$): Experiment and theory, Annu. Rev. Nucl. Part. Sci. 62, 237 (2012).
  * Bennett  _et al._ [2007] G. W. Bennett _et al._ (Muon $g-2$ Collaboration), Statistical equations and methods applied to the precision Muon ($g-2$) Experiment at BNL, Nucl. Instrum. Methods Phys. Res., Sect. A 579, 1096 (2007).
  * Farley [1972] F. J. M. Farley, Pitch correction in ($g-2$) experiments, Phys. Lett. B 42, 66 (1972).
  * Farley and Picasso [1990] F. J. M. Farley and E. Picasso, The Muon $g-2$ Experiments, Adv. Ser. Dir. High Energy Phys. 7, 479 (1990).
  * Field and Fiorentini [1974] J. H. Field and G. Fiorentini, Corrections to the $g-2$ frequency in weak focusing storage devices due to betatron oscillations, Nuovo Cimento A 21, 297 (1974).
  * Silenko [2006] A. J. Silenko, Equation of spin motion in storage rings in the cylindrical coordinate system, Phys. Rev. ST Accel. Beams 9, 034002 (2006).
  * McMillan [1975] E. M. McMillan, Multipoles in cylindrical coordinates, Nucl. Instrum. Methods 127, 471 (1975).