Beam dynamics corrections to the Run-1 measurement of the muon anomalous magnetic moment at Fermilab

This paper presents the beam dynamics systematic corrections and their uncertainties for the Run-1 data set of the Fermilab Muon g-2 Experiment. Two corrections to the measured muon precession frequency $\omega_a^m$ are associated with well-known effects owing to the use of electrostatic quadrupole (ESQ) vertical focusing in the storage ring. An average vertically oriented motional magnetic field is felt by relativistic muons passing transversely through the radial electric field components created by the ESQ system. The correction depends on the stored momentum distribution and the tunes of the ring, which has relatively weak vertical focusing. Vertical betatron motions imply that the muons do not orbit the ring in a plane exactly orthogonal to the vertical magnetic field direction. A correction is necessary to account for an average pitch angle associated with their trajectories. A third small correction is necessary because muons that escape the ring during the storage time are slightly biased in initial spin phase compared to the parent distribution. Finally, because two high-voltage resistors in the ESQ network had longer than designed RC time constants, the vertical and horizontal centroids and envelopes of the stored muon beam drifted slightly, but coherently, during each storage ring fill. This led to the discovery of an important phase-acceptance relationship that requires a correction. The sum of the corrections to $\omega_a^m$ is 0.50 $\pm$ 0.09 ppm; the uncertainty is small compared to the 0.43 ppm statistical precision of $\omega_a^m$.

1 dataset of the Fermilab Muon g−2 Experiment. Two corrections to the measured muon precession frequency ω m a are associated with well-known effects owing to the use of electrostatic quadrupole (ESQ) vertical focusing in the storage ring. An average vertically oriented motional magnetic field is felt by relativistic muons passing transversely through the radial electric field components created by the ESQ system. The correction depends on the stored momentum distribution and the tunes of the ring, which has relatively weak vertical focusing. Vertical betatron motions imply that the muons do not orbit the ring in a plane exactly orthogonal to the vertical magnetic field direction. A correction is necessary to account for an average pitch angle associated with their trajectories. A third small correction is necessary, because muons that escape the ring during the storage time are slightly biased in initial spin phase compared to the parent distribution. Finally, because two highvoltage resistors in the ESQ network had longer than designed RC time constants, the vertical and horizontal centroids and envelopes of the stored muon beam drifted slightly, but coherently, during each storage ring fill. This led to the discovery of an important phase-acceptance relationship that requires a correction. The sum of the corrections to ω m a is 0.50 ± 0.09 ppm; the uncertainty is small compared to the 0. 43  The measurement of the muon magnetic anomaly 1 a µ ≡ (g µ − 2)/2, where g µ is the factor describing the relationship of the muon magnetic moment to its spin, has undergone significant development since the late 1960s 1 aµ is a scalar quantity and thus not technically the "anomalous magnetic moment" that is ubiquitous in the literature.
when the idea of using a magnetic storage ring for the measurement was first introduced [1]. Two storage ring experiments at CERN [1,2] and the Brookhaven National Laboratory (BNL) Experiment (E821) [3] have increasingly refined the technique, leading to a determination of a µ to a precision of 0.54 ppm [3]. These experiments determine a µ by measuring the muon spin precession frequency relative to the momentum vector while a muon beam is confined in a storage ring. In an ideal case with muons orbiting in a horizontal plane with a uniform vertical magnetic field, the anomalous precession frequency ω a is given by the difference between the spin (s) and cyclotron (c) frequencies, ω a = ω s − ω c . The observable ω a is proportional to a µ . A highly precise determination of a µ is motivated by the fact that the standard model (SM) prediction [4] is known to 0.37 ppm and thus a direct comparison with the experiment serves as an excellent and sensitive test of its completeness. The E821 a µ measurement [3] is larger than the SM prediction; the level of significance is 3.7 standard deviations (σ). This difference can be interpreted as a hint of new physics or an indication of unaccounted-for systematic errors in the experiment or theory. The Fermilab Muon g −2 Experiment (E989) is designed to test the validity of the BNL result and to go further by improving on the experimental precision.
The principle of an a µ measurement in a storage ring lies in the expression where the experiment measures the two quantities in the ratio ω a /ω p (T r ), and the fundamental constantsthe proton-to-electron magnetic moment ratio [5], the QED factor µ e (H)/µ e [6], the muon-to-electron mass ratio [6,7], and the electron g e factor [8]-are obtained from external measurements and calculations. Reference [9]-a companion to this paper-reports the first results from the analysis of the Run-1 data collected in 2018: a µ (FNAL) = 116 592 040(54) × 10 −11 (0.46 ppm).
The statistical (0.43 ppm), systematic (0.15 ppm), and external fundamental constant (0.03 ppm) uncertainties are combined in quadrature. The result is in good agreement with the BNL result and, when the two measurements are combined and compared to the SM, the significance of the difference between experiment and theory rises to 4.2 σ. In order to obtain this result, a number of beam dynamics effects had to be accounted for in the analysis. This paper reports the methods used to quantify the impact of these effects and their uncertainties. These effects were part of a larger analysis, which can be summarized briefly as follows.
In Eq. 1, termω p (T r ) represents the magnetic field sampled by the muons and measured using pulsed proton NMR. The magnitude of the field is calibrated in terms of the equivalent precession frequency ω p (T r ) of a proton shielded in a spherical sample of water at a reference temperature T r = 34.7 • C. The detailed analysis of the magnetic field determination for this experiment is reported in a second companion paper [10].
The anomalous precession frequency ω a in the numerator of Eq. 1 is extracted from a fit to the time and energy spectrum of positrons, which encodes a modulation of its intensity at the anomalous precession frequency. Importantly, the frequency extracted by fitting the data is the measured quantity ω m a , and not the desired quantity ω a needed to obtain a µ in Eq. 1. Details of the measured precession frequency ω m a analysis, including statistical and systematic uncertainties, are reported in the final companion paper to this report [11].
The muons in the experiment do not orbit in a perfectly horizontal plane in a homogeneous vertical magnetic field. An electrostatic quadrupole (ESQ) system is required to provide relatively weak vertical focusing. The resulting presence of an electric field and vertical betatron oscillations, as well as muons lost during the measurement period, and a faulty ESQ component that led to a time-dependent optical lattice necessitate "beam dynamics" corrections that must be applied to ω m a in order to obtain ω a , the quantity needed to determine the magnetic anomaly.
Two corrections are a result of a physical reduction in the spin precession frequency of the muons. The electric-field correction C e is a result of muons traveling orthogonally through the radial electric field from the ESQ system and experiencing a motional vertical magnetic field. The pitch correction C p compensates for the vertical betatron motion that causes, on average, a slightly nonorthogonal relation between the muon momentum and the vertically aligned magnetic field. While these were known from previous experiments, two additional corrections must be made to compensate for a dynamic change to the stored muon ensemble's average spin phase during each injection cycle of muons into the ring. The C ml correction is associated with muons lost during storage that have a slightly different spin phase compared to those that remain stored. Finally, because two high-voltage resistors in the ESQ network had longer than designed RC time constants, the vertical and horizontal centroids and envelopes of the stored muon beam drifted slightly, but coherently, during each storage ring fill. This led to the discovery of an important phaseacceptance relationship that requires a correction C pa . The four corrections are applied in a linear combination 2 to ω m a , ω a ≈ ω m a (1 + C e + C p + C ml + C pa ) .
The total shift to the measured ω m a value is (+0.50 ± 0.09) ppm. Although these small corrections are all at 2 Cross terms between corrections are neglected here which is more than sufficient.
the sub-ppm level, the net correction exceeds the statistical uncertainty on ω m a , and thus the corrections must be scrutinized carefully and their uncertainties determined precisely. A thorough discussion is presented in this paper.
Section II includes background information on the formal experimental principle, Run-1 dataset characteristics, injection and storage information, and descriptions of the experimental instrumentation, magnetic field, and the simulation tools. The electric-field correction is evaluated using the stored muon momentum distribution that is obtained via a data-driven procedure described in Sec. III. The C e correction and its uncertainty can be found in Sec. IV. The pitch correction C p is discussed in Sec. V.
The time-dependent nature of the detected oscillation phase which compels the final two corrections was exacerbated during Run-1 due to damage to 2 out of 32 high-voltage resistors in the ESQ charging system. Section VI examines the consequences of this unexpected configuration for the time-dependent electric quadrupole field and the dynamics of the stored beam. Section VII discusses how this slow change of the steering field leads to larger than expected muon losses and evaluates the resulting correction C ml due to a correlation between muon momentum and initial average spin phase. Section VIII decribes the phase-acceptance correction C pa that is introduced because the average spatial distribution of the beam drifts during storage as a result of the same ESQ resistor issue.
Finally, we summarize and conclude in Sec. IX and provide four appendices covering details of electric field and pitch corrections, the generation of phase, asymmetry, and acceptance maps needed for the C pa correction, and a typical fitting function used to determine ω m a .

A. Experimental principle
Polarized positive muons are injected into a 14.2 mdiameter storage ring having a highly uniform ∼1.45 T vertical magnetic field; each individual injection sequence is termed a "fill". For muons of momentum p µ orbiting in a horizontal plane that is perpendicular to a perfectly uniform magnetic field ( B · p µ = 0), the magnitude of the cyclotron frequency is ω c = |q|B/mγ, where q, m, and γ are the muon charge, mass, and Lorentz factor respectively. The torque on the magnetic moment, together with the Thomas precession, rotates the muon spin at the frequency ω s = g µ |q|B/2m + (1 − γ)|q|B/γm [12]. The anomalous precession frequency ω a can then be described as the difference Because of parity violation in the µ + weak decay, positrons are emitted with an energy and an angular distribution that are each highly correlated to the muon spin direction in its rest frame. Ignoring effects from beam dynamics, the detectors that are distributed uniformly around the inside of the storage ring see a positron count rate versus time in fill t and positron energy E with the following functional form: The normalization, time-dilated muon lifetime, asymmetry, anomalous precession frequency, and phase constant at the time of injection are represented in Eq. 4 by N 0 (E), γτ µ , A(E), ω a and ϕ 0 (E), respectively. The energy dependence of N 0 and A has its origin in the Lorentzboosted muon decay (Michel) spectrum; the actual values depend on the detector acceptance. The asymmetry also depends on the muon ensemble average polarization magnitude P µ , which is ≈ 95% in this experiment. In practice, ω a is determined from a fit to a positron time spectrum having a range of positron energies. The detector design is optimized to accept higher-energy positrons for which the figure of merit N A 2 is maximized.
We emphasize that, throughout this paper, the variable ϕ 0 represents the phase constant in the timedependent phase angle of the cosine in Eq. 4. Its value is determined by the fitting procedure and need not, and cannot, be determined with sufficient precision a priori to the fit of the spin precession data sample. The phase constant ϕ 0 is dependent on the stored muon momentum p µ , the decay positron energy E, and the transverse decay coordinates 3 (x, y) inside the storage ring. The N 0 and A terms in Eq. 4 are also functions of the same quantities, but they couple much more weakly to ω m a . The physical interpretation of these dependencies of the phase constant and their implications to systematic uncertainties in the determination of ω a are most important in the determination of the C ml and C pa correction factors discussed in Secs. VII and VIII, respectively.
The relevant observable for the anomalous precession frequency is the oscillation in the quantityβ · S, where S is the spin vector of the muon. Following the Thomas-Bargmann-Michel-Telegdi (Thomas-BMT) equation [13] for spin and the Lorentz equation for momentum, one observes that, in the presence of an electric field E, where S T is the component of S perpendicular toβ and we have ignored a possible term owing to a nonzero muon electric dipole moment, which has been determined to be negligible [14]. This expression corresponds exactly to the muon spin precession in the ring. The electric field term in Eq. 5 vanishes for muons having the "magic" momentum p 0 = 3.094 GeV/c (γ ∼ 29.3). The experiment is therefore designed around injection and storage of muons centered on p 0 .

B. The Run-1 dataset
The muon delivery and experimental apparatus were commissioned from June 2017 to March 2018. The Run-1 data-taking period began on March 26, 2018 and concluded on July 7, 2018. The datasets in this analysis comprise 14.13 M muon fills. Approximately 5000 muons are stored in each fill, at an average of 11.4 fills per second.
Vertical focusing in the storage ring is achieved using a suite of ESQ plates that occupy 43% of the ring circumference. The field index n, responsible for the relatively weak focusing in the vertical direction, is defined by where R 0 is the central orbit radius, v is the muon velocity, B 0 is the magnetic field, and the gradient in the effective vertical electric field is determined from the plate voltages and geometry. Four distinct datasets, referred to hereafter as 1a, 1b, 1c, and 1d, have been separately analyzed corresponding to four different combinations of kicker magnet and ESQ voltages (see Table I). The separate analyses of these datasets yield consistent results for a µ . The 18.3 kV (n = 0.108) and 20.4 kV (n = 0.120) values were chosen to avoid storage ring betatron resonance conditions that would lead to large muon losses. based on the asymmetry-weighted analysis method [11]; the demand ESQ voltage; the measured field index n which is affected by the damaged ESQ resistors (Sec. II C); and the typical kicker strength in kilovolts used to deflect the incoming beam into the storage volume.
a The precession fit start time for Run-1d was delayed to 50 µs, in contrast to 30 µs for the other data groups.
The different field indices lead to differing beam frequencies, since the horizontal and vertical betatron tunes for a uniform set of quadrupole fields that occupy the full azimuth of the storage ring are given by with corresponding betatron frequencies These expressions are sufficiently accurate for our purposes. Calculations of the frequencies for the different field indices are shown in Table II. Coherent radial and vertical betatron oscillations of the centroid and width of the stored beam are driven by a combination of the mismatch between the beam line admittance and the storage ring acceptance, the intrinsic divergence of the incoming beam, and the strength of the storage ring kicker system. The radial coherent betatron motion of the muon ensemble introduces an oscillatory time dependence to the N 0 , A, and ϕ 0 terms in Eq. 4.
Since ω x > ω c /2, the observed frequency at each detector is the aliased frequency ω CBO = ω c − ω x . This coherent betatron oscillation (CBO) occurs inside the storage volume but becomes imprinted on the positron count spectrum N (t) because of the radial dependence of the detector acceptance. The radial mean of the muon distribution is modulated at frequency ω CBO and the radial width (rms) at 2ω CBO and also at ω CBO if the stored beam is not centered in the aperture. A smaller but similar effect exists for the vertical oscillations, where ω y is observed without aliasing, but the width oscillations at 2ω y are aliased to ω V W = ω c −2ω y . Because of the symmetric nature of the vertical detector acceptance, the ω V W effect is stronger than that from ω y . The calculated values of ω CBO and ω V W are also presented in Table II. These effects impact N 0 , A, and ϕ 0 and are accounted for by additional terms in the fitted positron decay time spectrum (see Appendix D and Ref. [11]). The field indices have been chosen such that the modulations do not couple strongly to ω m a . If the ESQ system does not maintain stable voltages within a fill, ω CBO and ω V W vary as a function of time in fill, which was the case in Run-1 but fixed thereafter. Muons from the Fermilab accelerator complex are created as follows. Bunches of 8-GeV protons strike a target in the AP0 building of the Muon Campus [15]. The average number of protons incident on the target per muon fill was 9.84 × 10 11 . Positive 3.1 GeV/c particles are extracted and transported through the 279-m-long M2/M3 FODO 4 beam line. Roughly 80% of the pions decay to muons in this beam line section, which is optimized to produce a muon beam with an average longitudinal polarization of approximately 95%. The particles then enter the 505-m-circumference Delivery Ring, where they are allowed to circulate for four full turns. The accompanying 3.1 GeV/c protons from the target station lag behind the faster muons and are swept out by an in-ring fastkicker magnet. The purified and polarized muon bunch is then directed along the M4 and M5 beam lines and into the storage ring. The net 2463 m path length from the target to the storage ring reduces the initial pion intensity by a factor 1.5 × 10 6 , an important feature that eliminates the hadronically induced flash at injection, which was challenging in the E821 experiment [3].
Sixteen individual bunches of muons are injected every 1.4 s cycle in two sequences of eight with 10 ms separation and 267 ms between the start of each sequence. The average temporal intensity distribution of muons at the entrance to the storage ring is shown as a dashed red line in Fig. 1. The shape of this intensity-time distribution varies slightly for each of the eight fills in a group. Figure 2 shows a plan view of the main components and detectors in the storage ring. Not shown for clarity in the figure is the "C-shaped" superconducting storage ring magnet [16]. It was transported from BNL to Fermilab, rebuilt, and reshimmed by our Collaboration to a roughly threefold improved uniformity. The NMR probes used to measure the absolute and relative magnetic field are all new, but they closely follow designs developed for E821 [17]. The superconducting inflector magnet [18], the ESQ hardware [19], and the vacuum chambers are used "as is" from E821. The ESQ power supply system and its controls, the entire storage ring kicker system, and all detectors and associated electronics were custom designed for E989.
The storage ring is designed to accept muons in a narrow momentum range around p 0 ; however, the incoming beam has a comparatively wide momentum spread of ±1.6% and no upstream dispersive focus. Only a few percent of the particles are stored with a momentum width less than ±0.15%. The muon bunch enters the ring through a nearly field-free corridor-18 mm wide × 56 mm high × 1700 mm long-provided by the superconducting inflector magnet [18]. Three kicker stations (K1, K2, and K3) create a vertical magnetic field that opposes the main storage ring field and, thus, deflects the muons passing through them outward by ∼10 mrad. Ideally, this transient field is turned off before the muons complete a full turn in about 149 ns. A challenge in this experiment is creating such a transient magnetic kick and timing it with respect to each muon bunch. Figure 1 also shows a sample of the measured kicker-induced magnetic field overlaid with the muon bunch. The kicker pulse shape and magnetic field strength are discussed later in the context of the stored muon momentum distribution within the ring. The kick is not uniform across the full temporal extent of the incoming muon bunch, a fact that impacts the total storage efficiency and the momentum distribution of the stored beam. FIG. 1. The solid black trace represents the shape of the kicker magnetic field. The ringing is owing to an unavoidable impedance mismatch. The dashed red trace represents the average injected muon bunch temporal intensity distribution as measured by the T0 detector. The vertical dotted lines are separated by the 149 ns cyclotron period. The overlay between kicker and beam shapes shows that not all entering muons experience the same magnetic deflection. It is clear from the ringing that muons experience multiple kicks, both inward and outward.
Four ESQ stations are symmetrically placed around the ring. Each consists of a long (L) and a short (S) section spanning 26 • and 13 • , respectively. Thus Q1 in Fig. 2 has sections Q1L and Q1S, each with four plates that are connected to power supplies through individual high-voltage resistors. The plates are raised from ground to operating voltage prior to each fill with RC charging time constants of ∼5 µs. They are returned to ground at the end of the 700-µs-long fill. Plates are charged using either one-step or two-step power supplies. Those plates connected to the two-step supplies rise to a preset voltage that is set to be ∼5-7 kV below their final de-signed voltage. After a programmable delay of approximately 7 µs, they are raised to the full set point voltage. This procedure, known as scraping, initially displaces the beam vertically and horizontally with respect to the central closed orbit. When a muon's horizontal and vertical oscillations conspire to exceed a 45 mm-radius with respect to the quadrupole center, it will likely strike a collimator, scatter, and lose energy. Such muons leave the storage volume in a few turns. When the voltages are symmetrized, the muon distribution relaxes back to the nominal center, where muon losses are minimized. This scraping process is designed to be completed by 30 µs after injection, before the nominal measurement start time.
The charging traces of the ESQ plates for the Run-1 configuration are shown in Fig. 3. As noted, 30 of 32 charging profiles follow the nominal and ideal pattern as shown in the solid black and dotted red traces for the one-step and two-step supplies, respectively. However, following the completion of Run-1 data taking, it was discovered that two resistors had dynamically changed their resistances when high voltage was applied. The resulting traces, measured after the data-taking period, are shown in dotted blue and solid orange in Fig. 3. These resistors were connected to the upper and lower plates of the Q1L system. Because they are asymmetric and did not rise to their proper voltage prior to the fit start time, they introduced a perturbation to the stored muon spatial distribution versus time in fill. As will be evident in later sections, this was likely responsible for the larger-than-expected rate of lost muons, and it created a time-dependent phase shift from the correlations of decay position to average phase owing to the detector acceptance. The changing voltages within the measurement time also led to a change in the storage ring tune values, the consequence of which is a dependence of the CBO frequency versus time in fill. The ω CBO (t) variation is well measured and included in the fits that determine the anomalous precession frequency. These damaged resistors were replaced prior to Run-2.
The temporal bunch structure of the injected muon beam does not fill the storage ring uniformly. Observing how this bunch spreads out owing to the finite momentum distribution is central to determining the stored muon momentum distribution (see Sec. III). For Runs-1a, 1b, and 1c, the measuring period to determine ω m a begins 30 µs after injection, an optimization of minimizing systematic uncertainties and maximizing statistical significance. During Run-1d, the damaged resistors deteriorated further, causing greater perturbations to the stability of the muon storage distribution at early times in the fill. A delayed fit start time of 50 µs is used for this dataset, allowing the ESQ plate voltages to more closely reach to their target values. IG. 2. Plan view of the g −2 storage ring vacuum chamber and instrumentation highlighting the important components discussed in this paper. The beam enters through a fieldfree corridor provided by the inflector magnet. Three kicker stations operate in concert to deflect the beam onto a stable orbit during the first turn. Four ESQ stations, each consisting of a short and long section, provide vertical containment. The T0 and IBMS detectors monitor the injected beam intensity, time profile, and spatial distribution.

D. Instrumentation employed to study the muon beam and decay positrons
The detectors employed to measure the incoming and stored muon bunches are the T0, Injected Beam Monitoring System (IBMS), Calorimeter, and Tracker systems. They are orientated as shown in Fig. 2. Their signals are recorded using dead-time-less digitizers and saved fill by fill for offline processing. The injected muon bunch first passes through the T0 and IBMS detectors, located at the entrance to the storage ring. They are used to measure the intensity of the incoming beam and its temporal and spatial profile and to establish the average entrance time of the bunch on a fill by fill basis, which is required for the analysis of the momentum distribution as discussed in Sec. III. The T0 detector is a 1.0-mm-thick plastic scintillator with dual photomultiplier readout. The two IBMS detectors each consist of a horizontal and vertical array of 16 0.5-mm-diameter scintillating fibers.
Twenty-four electromagnetic calorimeters are positioned symmetrically around the inside radius of the storage ring, adjacent to the storage volume, but outside of the vacuum chamber; see ular aluminum window before striking a detector. Each calorimeter station consists of 54 lead-fluoride Cherenkov crystals read out individually on the downstream side by large-area silicon photomultipliers (SiPMs) [20,21]. The signals are continuously digitized at 800 megasamples per second. The precise time alignment of the 1296 crystals and the system gain stability are enabled using a laser system as described in Ref. [22]. Reconstruction of positron showers in the calorimeter crystals yields energy, time of hit, and impact position. The coincidence of signals in three consecutive calorimeters, each depositing an energy typical of a minimum ionizing muon of about 170 MeV, is used to identify muon losses and measure the loss rate versus time in fill.
The straw-tracker systems are located at approximately 180 • and 270 • with respect to beam injection. They reside within the vacuum chamber in the scallop region just upstream of a calorimeter but outside of the muon storage volume. Both stations consist of eight modules, each made of 128 5-mm-diameter straws oriented at ±7.5 • with respect to the vertical. They are used to track the decay positrons with the intention of tracing the decay trajectory back to its point of tangency inside the storage volume, a good proxy for the decay position of the parent muon. High-quality tracks are selected to construct distributions of the beam position at the tracker locations. These measurements are corrected for momentum-dependent detector resolution and the nonuniform acceptance of the detector. The magnitudes of the corrections are estimated using the GM2RINGSIM package described below. These trackers are critical to all beam dynamics topics discussed in this paper. They provide the muon profile versus time, which is used to determine key storage properties, such as the betatron horizontal and vertical frequencies, their time dependence owing to the damaged resistor influence, and the vertical distribution for the pitch correction.

E. Magnetic field considerations
The highly uniform magnetic field amplitude B = | B| within the storage volume allows for the extraction of the spatial field structure with NMR probe measurements [10]. The magnet iron poles and iron shims are precisely aligned in order to minimize the rms variations of the multipole terms as a function of the azimuth. The reduction of the field nonuniformities improves several systematic uncertainties in the measurement of the magnetic field, including the extraction of the frequency from the NMR probe measurements and the uncertainty due to limited knowledge of the probe positions. Priority is given to the reduction of the lowest-order multipoles that couple to the moments of the muon beam. Surface correction coils on the surfaces of the pole pieces are utilized to reduce the azimuthally averaged field multipole strengths. Typically, the azimuthally averaged multipoles are reduced to below ∼1 ppm, as shown in Table III.
The main measurement uncertainties in the mapping of the field, detailed in Ref. [10], stem from motional effects of the mapping device and the probe positions. Because of the overall field uniformity, despite the significant accuracy demands, the requirements for positioning of the probes are not detrimentally restrictive and can be met in practice. Specifically, a laser alignment survey of the mapper transiting around the ring and a sensitivity analysis determine the impact of the probe position uncertainties on the multipole strengths. This effect generates uncertainties for the azimuthally averaged dipole (12 ppb), normal quadrupole (27 ppb), and skew quadrupole (4 ppb). The azimuthal variations in the magnetic field and muon distributions are incorporated in the muon weighting analysis (< 20 ppb) as well as the beam dynamics simulations described in Sec. II F. The multipole strengths, tracked between field maps by stationary NMR probes outside the storage volume, are highly stable over time; the ranges observed in the Run-1 dataset are summarized in Table III. Additional experimental uncertainties from magnetic field tracking (56 ppb) and the presence of fast transients (99 ppb) are carefully quantified, resulting in total systematic uncertainties on the determination of the muonweighted magnetic field of approximately 114 ppb for the Run-1 dataset. The mean, standard deviation (SD), and peak-to-peak (Pk-Pk) range correspond to the variations of these azimuthally averaged multipole strengths over time and are significantly smaller than the variation in azimuth. The Pk-Pk of the dipole denote lower and upper deviations relative to its mean. The uncertainty on the muon-weighted dipole term is approximately 114 ppb for the Run-1 dataset [10].

F. Beam dynamics simulation tools
Many of the results discussed in this paper incorporate comparison with, or results from, beam dynamics simulations. Several compact simulation packages were used to rapidly estimate effects, but it is the work of three sophisticated and complementary methods that drive the results. To determine critical information, such as momentum-time beam correlations, all three are typically used. In all cases, it is critical that the simulation programs are cross-checked against a set of benchmarks showing that they can evaluate analytically calculated effects with high precision [23]. It is also imperative that they are first tuned to match measurements of the incoming beam properties, the storage distribution within the ring, the CBO amplitudes, frequencies, and their time dependence, and the stored momentum distribution. We describe their essential features below.

GEANT4-based GM2RINGSIM
The GM2RINGSIM program [24] is a GEANT4 [25] based model of the storage ring and the final focus beam line used to steer the beam into the ring. The model includes all of the active detectors and most of the passive components installed in the storage ring. The geometry is constructed from a mixture of GEANT4 native solids and CADMESH objects [26]. Calorimeters and straw-tracker modules are fully described, as are complex solids, such as the vacuum chambers and their inner structures, and the ESQ and kicker plates. The geometry for the strawtracker modules includes coordinates as determined in alignment surveys.
Four sets of realistic time-and space-dependent electric and magnetic fields are implemented. These include a pure dipole magnetic field in the storage region and a radially dependent fringe field that extends toward the center of the ring. Radial magnetic field maps and additional multipole perturbations in the storage region, as described in Sec. II E, can optionally be included. The inflector magnetic field is implemented as a map. The time dependence of the kicker magnetic field is taken from direct magnetometer measurements made at the center of the plates (see Fig. 1). The spatial field map within the kicker region is obtained using the finite-element magnetics modeling package OPERA [27]. The strength and timing of the kicker with respect to the injected beam can be adjusted independently for each of the three kicker plates. The fields associated with the ESQ plates are implemented as a multipole expansion. They are dynamically evolved through the scraping periods at the beginning of each muon fill, and they can accommodate the perturbations and independent time constants caused by the damaged resistors during Run-1.
Simulated datasets are typically generated using two types of "particle guns." The beam gun imports muon distributions at the end of the final focus beam line as determined from G4BEAMLINE simulations of the Muon Campus beam and injects them into the storage ring. The beam gun allows for the injection of a mixture of particles, facilitating studies of proton and positron contamination. The gas gun omits the injection process completely and instead fills the storage region phase space of the ring with muons that then decay at that location. This distribution of muons is matched to reflect the measured vertical and radial offsets of the beam within the storage ring as well as the measured coherent betatron oscillation amplitudes. The gas gun is particularly useful for positron acceptance and reconstruction studies. In both cases, muon spin is appropriately evolved during time in fill, and proper spin-dependent muon decays are employed.

COSY INFINITY
The COSY-based model [28] is a data-driven computational representation of the storage ring. Dedicated packages in COSY INFINITY [29] for the design and analysis of particle optical systems provide the framework for beam physics studies and symplectic tracking simulations in the storage ring. The beam dynamics of the injected muon beam is recreated with high fidelity by representing the magnetic and electric guide fields in the storage region based on measurements of the beam.
The magnetic field inhomogeneities of the storage ring are determined from the magnitude of the magnetic field as measured by the NMR probes (see Sec. II E). Because of the high uniformity of the magnetic field along the vertical direction in the storage region, a reliable extraction of magnetic multipole strengths from the experimental data is performed and implemented in the model as a series of magnetic multipole lattice elements. Each ESQ station is modeled as an optical element su-perimposed on the magnetic field. The nonlinear action of the ESQ on the beam's motion is captured by accounting for the high-order coefficients of the electrostatic potential's transverse Taylor expansion. These coefficients are calculated by recursively iterating the horizontal midplane coefficients-modeled with conformal mapping methods [30]-to satisfy Laplace's equation in curvilinear optical coordinates. The effective field boundary and fringe fields of the ESQ are calculated [30] using COULOMB's [31] boundary element method field solver.
With the aforementioned electric and magnetic guide fields implemented in the COSY-based model, the orbital and spin equations of motion are well defined and integrated to produce high-order transfer maps via differential algebra methods [32]. Transfer maps are computed and combined to recreate either azimuthal segments or the entire storage ring.
Beam-tracking simulations are performed by preparing transfer maps of the storage ring. The muon beam, represented as an array of orbital and spin coordinates around the closed orbit, is transformed turn by turn with transfer maps that encapsulate the time-evolving guide fields as they vary throughout the beam fill. Symplecticity is enforced during beam tracking with high-order transfer maps to account for the truncation of components beyond the map order, in this way controlling energy conservation and error propagation. To account for beam collimation, special routines were developed to efficiently remove muons beyond the collimator apertures and only at the azimuthal locations where they are inserted. This tool, together with the symplectic enforcement during tracking, allows for reliably studying muon loss rates in the storage ring.
The beam conditions after the action of the injection kickers are obtained by preparing high-order transfer maps of the azimuthal segments where the three kickers are placed. Using these maps combined with the transfer maps of the other components of the storage ring, a beam distribution obtained from simulations of the Muon Campus beam lines and the inflector is transferred into the ring's storage region. The structure and time-dependent strength of the kicker magnetic fields are constructed from magnetometer measurements (see Fig. 1). Alternatively, the initial beam distribution after injection is also prepared by calculating with a non-negative leastsquares solver the probability density functions of the beam, based on data from the straw-tracking detectors and the fast-rotation analysis (see Sec. III).
The COSY-based model has been extensively used to calculate lattice configurations (e.g., periodic Twiss parameters, betatron tunes, closed orbits, and dispersion functions) of the storage ring with and without the damaged ESQ resistors. It provided a reconstruction of the special electric field behavior of the ESQ electric fields during Run-1 for the assessment of beam dynamics systematic uncertainties. Further, it has provided a numerical model to link tracker measurements with direct highvoltage (HV) probe measurements of the damaged ESQ resistors, and it has been used to choose optimal configurations of the storage ring to minimize muon losses.

BMAD
BMAD refers to a subroutine library [33] for simulating the dynamics of relativistic beams of charged particles and an associated format for defining beam line elements. So defined, the full complement of the analysis tools of the library can be used to investigate the particle dynamics. Particle-tracking methods include Taylor maps, symplectic integration, and Runge-Kutta integration through field maps. Taylor maps can be predefined or constructed by tracking.
The BMAD formatted representation of the g−2 experiment is comprised of three distinct branches: i) the M5 beam line; ii) the injection channel and inflector; and iii) the storage ring, including a static magnetic field, timedependent quadrupole electric fields and kicker magnetic fields, and collimator apertures. The beam lines are assembled as a sequence of elements with fixed length. The electromagnetic fields in each element are defined by field maps, multipole expansions, or analytic expressions. Time dependence for pulsed kickers and ramped ESQ plates requires custom code. Specialized custom routines are used to incorporate arbitrary kicker pulse shapes or ESQ voltage time dependence. The curvilinear coordinate system can be represented by beginning with a full three-dimensional map and then extracting an azimuthal slice of the map (in x, y at fixed φ) or with a fitted multipole expansion of McMillan functions. The curvature necessarily introduces nonlinearities that are not faithfully included in a two-dimensional Cartesian expansion. The BMAD code allows specification of the quadrupole electric field in terms of field maps or multipoles. The main magnet is represented as a map or analytic function with uniform field. A uniform radial component can be specified. Measured field errors are incorporated analytically. The azimuthal dependence of the error field is expanded as a solution to Laplace's equation in cylindrical coordinates in order to ensure consistency with Maxwell equations. The magnetic field through the hole in the back-leg iron and cryostat and main magnet fringe field is based on a three-dimensional OPERA [27] map. A distinct map is computed for the field in the inflector. The fringe and inflector maps are superimposed as appropriate.

III. DETERMINATION OF THE STORED MUON MOMENTUM DISTRIBUTION
The contribution to ω m a from the electric field depends on the momentum distribution of the stored muon beam or, equivalently, the equilibrium radial distribution. Since the azimuthal speed of the stored muons is nearly uniform over the 9 cm aperture, to a good approximation a muon's rotation frequency is inversely proportional to its equilibrium radius x e . Because the muons are stored over a range of x e , a beam bunched at t = 0 will steadily debunch, as the higher-frequency muons at smaller radii advance with respect to the lower-frequency muons at larger radii.
A technique [34] based on Fourier transformation yields a frequency spectrum that can be converted to radius and momentum. An alternative method extracts the radial distribution by a direct fit to the debunching signal of the muon beam [1,35]. In both cases, the input data are provided by the calorimeters, which measure the time dependence of the intensity of the decay positron distribution. The positron counts from the 24 calorimeters are merged together with a time offset of T /24 per calorimeter where T is the cyclotron period, approximately 149 ns. The upper panel of Fig. 4 shows the intensity 4-5 µs after injection; the lower panel expands the time range to 4-14 µs. The individual turns around the ring, referred to as "fast rotation," are distinct. The slower modulations in the upper envelope are caused by muon decay and the ω m a precession frequency. By the nominal precession fit start time of 30 µs, the rapid cyclotron frequency modulation has largely dephased and eventually disappears. To isolate this fast rotation signal, the calorimeter time spectra are divided by a fit 5 to the envelope modulation using Eq. 4, which effectively removes all the significant physics signals except the fast rotation itself.

A. Momentum-time correlation
The Fourier and debunching methods, discussed below in Secs. III B and III C, both assume that the momentum distribution in the captured muon pulse is uncorrelated with the longitudinal position in the pulse train (see Fig. 1). If the momentum and time are uncorrelated, the distribution is maximally bunched when it enters the ring. The peak intensity of the fast rotation signal only decreases with time. The behavior is symmetric with respect to time reversal. One can imagine that the peak intensity likewise decreases moving backward in time.
In reality, a momentum-time correlation is introduced by the injection kicker. The efficiency with which muons are captured in the ring depends on momentum and the amplitude of the kicker. If the kicker field is low, acceptance is higher for high-momentum muons; if the kicker field is high, low-momentum muons are favored. Since the magnetic field of the kicker varies over the time duration of the incoming muon pulse, so do the momenta of muons that are stored in the ring. Simulations are used to characterize the momentum-time correlation of the captured muons and to estimate the systematic uncertainty in the measurement of the momentum distribution from the fast rotation analysis.

B. Frequency domain fast rotation analysis: Fourier method
For a stored muon beam in which the cyclotron frequencies and injection times are independent, the frequency distribution of the ensemble can be extracted by the cosine transform of the fast rotation signal S(t) [34]: where t 0 is an effective time of symmetry for the ensemble. The optimal t 0 is determined by imposing that the transformŜ(ω) must vanish in the unphysical frequency region outside the range that can be stored. The problem is complicated by the fact that the first few microseconds of the fast rotation signal are contami-nated by beam positrons and by muons that will not be stored. The analysis is based on the intensity signal that begins at about 25 turns, or t s = 4.1 µs after injection, by which time the beam positrons have been lost and the performance of the calorimeter SiPMs has largely recovered from the intense flash of particles at injection. The available cosine transform is, therefore, missing the contribution before the start time t s , which introduces a background modulation B(ω) to the frequency spectrum. As shown in [34], this background can be estimated using an inverse cosine transform and a model for the frequency distributionŜ 0 with the result where the limits of integration correspond to the physical range of frequencies that can be stored. An ansatz for S 0 (ω) is hypothesized and the parameters determined by a fit to the background. IfŜ 0 (ω) = Aδ(ω − ω B ), where A and ω B are fit parameters for amplitude and frequency, respectively, then Empirically, this model is found to give a good fit to the background for start times before 5 µs. More sophisticated functional forms give good background fits for start times up to t s ∼ 25 µs in Monte Carlo simulation, and, in practice, the fitted spectrum is nearly independent of start time for 4 µs < t s < 25 µs. An example of the frequency transformŜ of Eq. 9 is shown in Fig. 5, with the optimal background fit using Eq. 11. The background fit is subtracted, and the corrected transform is taken as the distribution of cyclotron frequencies.

C. Time domain fast rotation analysis: debunching method
An alternative approach, pioneered by the CERN II collaboration, is based on a simple model of the beam's debunching at early times in the measurement period [1,35]. Consider first the contribution to the fast rotation signal from a narrow bin in time and momentum space. The signal is initially rectangular but grows increasingly trapezoidal as it revolves around the ring, due to the uniform momentum spread within the bin. Mathematically, all the essential physics is captured in the propagator function β ijk , which describes this segment's contribution to the overall signal in the detector at time t j . Indices i and k identify the segment's equilibrium radius and position in time within the injected muon bunch, respectively. Using superposition, an ensemble of segments with joint distribution f i I k is given by The background is a consequence of the missing data before the start time of 4 µs. A cardinal sine (sinc) background fit has been applied to the black points and is subtracted from the whole frequency range as a correction.
Here, S j is the calorimeter signal in time bin j, f i describes the radial distribution, and I k describes the time profile of the injected beam. A least-squares fit to the signal S j determines f i and I k .
In the BNL experiment [3], many of the calorimeters were live on the first turn following injection. It was, therefore, simple to make a very good guess of the stored beam's initial time profile I k . As noted above in Sec. III B, detector signals in the Fermilab experiment cannot be used before 4 µs after injection. Therefore, the original CERN method was replaced with a pair of fits, for f i and I k , which are iterated until the results are stable. In the first pair of fits, the time profile is taken from the calorimeter signal at 4 µs, and the momentum distribution is determined by the fit. Then, that momentum distribution is used to update the injected time profile in a second fit. The determination of the momentum and time distributions are computationally identical. Between 50 and 100 iterations of this double fit are generally required for convergence.

D. The radial distributions for Run-1
An example of the radial distribution extracted by both the Fourier method and the debunching analysis is shown in the top panel in Fig. 6. The agreement is sufficiently good that either can be used to extract the electric field correction. The radial distributions for the Run-1a, 1b, 1c, and 1d periods as determined by the Fourier method are shown in the bottom panel in Fig. 6. The slightly different distributions are attributed primarily to the different kicker strengths, as indicated in Table I. The means of the radial distributions-in all cases-do not fall on the magic radius. This fact will enter in the calculation of the electric field correction in Sec. IV. The radial offsets and widths as determined by the Fourier method are given in Table IV. In both plots, the equilibrium radius is defined such that a magic-momentum muon is at 0 mm.

IV. ELECTRIC FIELD CORRECTION Ce
The electric field term in Eq. 5 produces a rest frame magnetic field that affects the measured anomalous precession frequency ω m a . For the simple case where we neglect the vertical betatron motion and ω a = ω s − ω c , Eq. 5 can be simplified to where ω m a is a scalar frequency and the subscripts r and y denote the radial and vertical components respectively. The electric field term vanishes at the magic momentum p 0 = (mc)/ √ a µ , or when E r = 0 which is the case at R 0 as a result of the design of the ESQ system. In practice, the stored muon distribution has a finite momentum spread and is not centered, as discussed in Sec. III and shown in Fig. 6. The mean radial electric field experienced by a muon oscillating about an equilibrium radius x e in an ideal electric quadrupole is where κ is the electric field gradient. Defining p = (p 0 + ∆p) and using where the approximation that the dispersion is a constant value of R 0 /(1−n) is sufficient, we average over the entire muon distribution to obtain the electric field correction where . The correction must be applied to the measured ω m a to obtain the anomalous precession frequency ω a used to determine a µ . The use of an analytic expression for the electric field correction has been previously considered [19,36,37]. The precision goals of our experiment require careful consideration of this correction. An extensive discussion of electric quadrupole nonlinearity factors is given in Appendix A and numerical tests to justify the validity of the use of Eq. 16 are given in Appendix B.
A. Measuring the muon radial distribution and calculating Ce The fast rotation analysis using the Fourier method (see Sec. III B) yields the distribution of cyclotron frequencies f , which are converted to equilibrium radii R by the relation R(2πf ) = v, assuming fixed muon velocity v. The radial offsets x e relative to the magic radius R 0 = 7112 mm are, therefore, x e = v/(2πf ) − R 0 . This conversion yields the distribution of equilibrium radial offsets, as in Fig. 6. The electric field correction depends on the mean and width of this distribution, via x 2 e = σ 2 xe + x e 2 . The recovered radial mean and width exhibit some variation when the fast rotation analysis is repeated over a range of positron energy bins. This is a consequence of variations in calorimeter acceptance with positron energy and radial decay position, supported by simulation with GM2RINGSIM. The final C e is, therefore, weighted over positron energy bins according to the statistical power of ω m a in each bin. For an asymmetry-weighted analysis, 6 this is the average of C e over positron energy bins between 1 and 3.1 GeV, weighted by N (E)A(E) 2 , where N (E) is the number of positron counts in energy bin E and A(E) is the fitted asymmetry of the ω m a modulation in that bin. The resulting corrections C e are tabulated with their dominant uncertainties in Table V, as discussed below.
The statistical uncertainties in x e , σ xe , and C e are estimated by repeating the fast rotation analysis over an ensemble of pseudodata. Each positron count N i in the measured time spectrum is shifted by ± √ N i at random, assuming uncorrelated Poisson statistics. The analysis is repeated over about 1000 randomly altered signals, and the ensemble standard deviation of each recovered quantity is taken as its statistical uncertainty. For subsets of the data with different total positron counts N , the statistical uncertainty was confirmed to scale as √ N . The fast rotation Fourier analysis method relies on several chosen parameters, including the start and end times of the cosine transform, the frequency bin spacing, the frequency distribution model used in the background fit, and the set of frequency bins included in that fit. In each case, a scan is performed over the range of appropriate choices, and the standard deviation in the result is taken as the corresponding systematic uncertainty. The total uncertainty attributed to these parameters is the average of their linear sum and quadrature sum, accounting for probable correlations.
Furthermore, Eq. 9 does not extract the frequency distribution perfectly as intended when there is any systematic relationship between cyclotron frequency and injection time. Referred to here as momentum-time correlation, this relationship is expected as discussed in Sec. III A. The corresponding uncertainty in C e has been estimated as 52 ppb, which is the average discrepancy between truth and reconstruction using GM2RINGSIM and BMAD simulations.
The expression for C e in Eq. 16 is derived under the assumption of continuous ESQ plates and ideal alignment of the quadrupole field relative to the target muon orbit. The effects of discrete ESQ plates, position misalignments, and voltage errors have been studied using The field index n has been measured using the value of ω CBO fitted during the production of the fast rotation signal and the relationship ω CBO = ω c −ω x = (1− √ 1 − n)ω c from Sec. II A. However, the CBO frequency is not constant throughout the fill, rather approaching a stable value with two exponential time constants based on the beam scraping procedure and damaged ESQ resistors. This time dependence has been measured using the tracking detectors (see Sec. VI A), and the uncertainty in n has been estimated using the rms spread in ω CBO over the ω m a measurement period.

V. PITCH CORRECTION Cp
The ESQ system used to vertically confine the muon beam creates vertical betatron oscillations, that is, periodic up-down pitching of the vector β. The a µβ × B term in Eq. 5 affects the value of ω m a . The magnitude of the term is reduced whenβ and B are not perpendicular, as is the case here. The vertical betatron frequencies ω y listed in Table II are an order of magnitude larger than the muon spin rotation frequency ω a , which avoids depolarizing spin resonances. These topics were first recognized and discussed in Refs. [39,40], where the pitch correction C p is derived as The vertical oscillation amplitude A can be extracted from measurements by the tracker detector system. The validity of Eq. 17 has subsequently been confirmed and explored further [23,41,42]. Appendix B describes our numerical simulations that also reaffirm Eq. 17 and furthermore permit modeling of the uncertainty on C p owing to the use of flat electrodes, ESQ plate misalignment, and voltage errors.
A. Measuring the muon vertical distribution and calculating Cp The vertical distribution of decay positrons is measured by the straw trackers, as shown for a subset of the data in Fig. 7a. Throughout Run-1, the temperature in the experimental hall slowly increased from 24.5 • C to 28.5 • C and exhibited typical diurnal fluctuations of roughly 1 • C. These changes produced a slowly varying radial component to the magnetic field, which caused the vertical mean of the beam to change with time. To account for this effect, these data were subdivided into shorter running periods, and a weighted average of their corresponding pitch corrections was computed. When determining the vertical muon distribution, time cuts are applied, which restricted the tracker data to the same time interval used for the measurement of the anomalous spin-precession frequency. The tracker measurements yield a good estimate of the true vertical distribution of the muon beam at the location of the tracker stations. However, it is necessary to take into account azimuthal variations around the storage ring from the discrete ESQ sections. The effective vertical distribution seen by the calorimeters that measure ω m a must also be determined. In order to address azimuthal variations in the vertical distribution of the muon beam, the vertical beta functions β y are evaluated as a function of azimuthal coordinate φ using GM2RINGSIM and COSY. By taking the ratio of β y (φ) with the value at each tracker station, the scale factor which relates the vertical width at each tracker station to any other azimuthal coordinate is obtained. The vertical distribution from a single tracker station is then averaged over azimuth by stretching or shrinking the width ratio in each azimuthal slice using β y (φ). Not all decay positrons hit a calorimeter and enter into the determination of ω m a . To account for this, each calorimeter's acceptance has been estimated as a function of transverse and azimuthal decay position using GM2RINGSIM (see Fig. 28). In each azimuthal slice of the storage ring, the transverse acceptance function is integrated over the radial dimension to produce an effective vertical acceptance function. Then, during the azimuthal averaging procedure, the vertical distribution in each azimuthal slice is masked by the corresponding calorimeter acceptance function after stretching by the vertical width ratio. The nominal results use the acceptance functions for all calorimeters combined, treating the acceptance per calorimeter as a cross-check. Furthermore, the calorimeter acceptance functions are subdivided by positron energy bin. Therefore, the pitch correction is evaluated using calorimeter acceptance from each positron energy bin, and the results are averaged according to the statistical power of ω m a in each energy bin. When considering calorimeter acceptance, the pitch correction can no longer be evaluated using y 2 as in Eq. 17. This is because the simple relation between the pitch angle and y 2 along an oscillation breaks down when the vertical positions are not evenly weighted. Instead, the right-hand side of Eq. 17 is used with A 2 from the distribution of oscillation amplitudes, which accurately describes the distribution of measured pitch angles when vertical acceptance is present. The amplitude distribution may be recovered from the trackers' vertical decay distribution by defining the fit function where N y j is the expected number of decays in the jth position bin, N A i is a fit parameter for the number of muons in the ith amplitude bin, and P (y j |A i ) is the calculable constant probability that a muon from the ith amplitude bin decays to the jth position bin. The amplitude distribution is extracted from a fit to tracker measurements of the vertical counts in each bin, N y j , after correcting for the intrinsic tracker acceptance and resolution. The fit result and corresponding amplitude distribution are shown in Figs. 7a and 7b, respectively. The amplitude distribution is stretched and averaged over azimuth as described above, and each amplitude bin A i is weighted by the average of the calorimeter acceptance over all position bins y j using P (y j |A i ). The result of these corrections is shown by the solid blue line in Fig. 7b. Finally, the pitch correction is calculated using Eq. 17 with A 2 .
Since the corrections have a small effect, the statistical uncertainty of C p can be estimated using the statistical uncertainty of the measured vertical width σ y , propagated to C p ∝ σ 2 y . The dominant systematic uncertainty is from the straw-tracker alignment and reconstruction. There are smaller contributions from the acceptance and resolution corrections. Other possible sources of uncertainty come from the estimation procedure for C p described above, as well as possible errors in alignment and calibration of the ESQ system that are described in detail in Appendix B. The resulting values for C p are summarized in Table VI. The corrections vary from 166 to 199 ppb with the main driver behind the range being the different ESQ settings for the different datasets. Differences with the same n value arise due to different radial magnetic fields. The statistical uncertainty is negligible, and the systematic uncertainty is well under control at the 12-14 ppb level.

VI. DYNAMIC EFFECTS OWING TO TIME-CHANGING FIELDS
A. Changes to the betatron frequencies The slower voltage increase on the Q1L upper and lower ESQ plates as a result of the damaged resistors in Run-1 introduces time dependencies of the storage ring lattice parameters (see the dotted blue and solid orange traces in Fig. 3). Electric normal quadrupole and skew dipole terms are largely proportional to the sum and dif-ference in high voltage between the top and bottom electrodes, respectively, making the beta functions, radial dispersion function, and the closed orbits time-dependent during the measurement period. Figure 8 illustrates the difference in the electrostatic potential in the Q1L ESQ section versus time in fill with respect to the nominal case for all other ESQ sections, which had normal resistors and stable voltages by 30 µs. The beta functions are a consequence of the focusing gradient configuration in the ring, which depends on the normal quadrupole terms at each ESQ section. The solid black line in Fig. 8 illustrates the added quadrupole field at Q1L that introduces a time dependence to the beta functions. The vertical closed orbit is distorted by skew dipole terms from the guide fields, the time dependence of which corresponds to the value of the dashed blue line in Fig. 8.    stored muon distribution at different times within the fill and are, therefore, used to measure the betatron oscillation parameters as well as any slow drift of the beam position or width. The betatron frequencies can be determined to high precision: In fact, it was analysis of tracker data that first drew attention to a possible time dependence of the electric quadrupole field. The measured betatron frequencies are necessary to verify the tune, and the measured betatron amplitudes are used in the optimization of kick strength and inflector deflection angle.
After the end of scraping, with stable voltages during the measuring period, ω CBO and ω VW should not change. However, in Run-1, they continued to evolve after 30 µs by 1.5% in Runs-1a, 1b, and 1c and by 3.0% in Run-1d. Figure 10 shows tracker measurements of ω CBO for Runs-1a and 1d. The fit function contains two exponen-tial terms that describe ω CBO (t) well. The first of these is a fast (∼7 µs) term, which relates to the changing field during scraping. The second term has a longer time constant of ∼60 µs in Runs-1a, 1b, and 1c and ∼80 µs in Run-1d. This term arises from the damaged resistors and the change in their average value during Run-1d due to their further deterioration.  The fit function and parameters are noted in the figure, together with their uncertainties. The difference between the two datasets is attributed to worsening performance of the damaged ESQ resistors. The time dependence of the frequency is included in the fits of the positron data from the calorimeters to accurately incorporate the CBO acceptance dependence.
The fits to the precession data require an accurate model of ω CBO (t) to obtain a good χ 2 and stable fit results when changing the fit start time. We note that the correlation between ω m a and ω CBO in the fit is small (∼2.4%). Varying the time dependence of the CBO frequency within allowable bands determined from the tracker measurements yields a less than 10 ppb shift to the fitted ω m a frequency. The effect from the VW is even smaller than that from the CBO, and it has a negligible effect on the precession fits.

B. Changes to the muon spatial distribution
Because the damaged resistors were connected to the same Q1L plates, and because the circuitry at station Q1L had RC characteristics that differed from the nominal charging conditions (see Fig. 3), we observe a physical effect on the mean and width of the muon distribution. The asymmetrical charging of the top and bottom plates in Q1L created an unbalanced quadrupole component of the ESQ field. This led to a time-dependent focusing gradient, resulting in submillimeter drifts in the beam widths and radial closed-orbit distortions from the timedependent optical lattice. The vertical closed orbit also manifested an in-fill temporal evolution owing to an introduced skew dipole field at Q1L. Figure 11 displays the time-varying vertical mean and rms measured by the 180°t racker for the Run-1d dataset, which is by far the worst case. In the range 30 − 300 µs, the vertical mean shifts by approximately 0.5 mm. Note that the absolute vertical scale is not known to better than ∼1 mm owing to alignment uncertainty and the local radial magnetic field. The vertical rms versus time in fill is unstable at the nominal fit start time of 30 µs but by 50 µs has flattened out, which is when the measurement period for this dataset starts. As noted previously, β y , which is proportional to y 2 , varies around the ring, so one does not expect the magnitude of the width change to be constant versus azimuth.
The tracker station measurements are limited to the 180°and 270°locations in the storage ring. However, the 24 calorimeter stations are positioned at regular intervals covering the complete azimuth. The segmented calorimeters provide a measure of vertical mean versus time in fill at each location. The COSY model of the vertical closed-orbit evolution during the fill predicts the vertical mean change around the ring. Similarly, GM2RINGSIM can determine the vertical mean change by using internal virtual tracking planes at different azimuthal locations. Figure 12 shows the calorimeter measurements of the vertical mean change from 40 to 300 µs and the corresponding predictions of the two simulation programs. The azimuthal variation in the data supports the implementation of the damaged resistors in the simulation and the projection of the dynamics at all azimuthal locations needed for the phase-acceptance correction discussed in Sec. VII B.

VII. MUON LOSS CORRECTION C ml
Several driving mechanisms can lead to loss of muons during storage. For example, the scattering of particles off of the residual gas in the vacuum chamber, noise from residual high-frequency electromagnetic fields in the sys- tem, the sampling of nonlinear fields near the aperture, and nonlinear resonances are potential mechanisms. In Run-1, the muon loss rate was higher than expected owing to a combination of factors including the damaged resistors and the nonoptimized kick. However, measurements still show integrated loss rates of less than a percent during the fill. In general, a muon will be scattered out of the storage region after it strikes one of a set of circular collimated apertures that limit the transverse phase space admittance and its momentum dependence. These collimators have an aperture of radius r 0 = 45 mm and are centered on the ideal orbit. Figure 13 shows where on the collimators these muons strike first, a consequence of circular apertures and normal amplitude distributions.
Monte Carlo beam line simulations and simple analytical calculations both predict that a correlation exists between the injected muon average spin phase and the particle momentum. If such a spin-momentum correlation exists, muons that permanently escape from the storage volume during data taking can potentially bias ω m a by inducing slow drifts in the phase. The ϕ 0 fit parameter in Eq. 4 depends on the average initial spin orientation of the muons that produce the detected decay positrons. In this context, an ensemble of muons is said to have a spinmomentum correlation if dϕ 0 /d p = 0, where p is the mean value of the muon momentum for the ensemble of decaying muons that produces the positron spectra being used to measure ω m a . The population of stored muons is depleted only by decays, while the population of muons that will be lost is depleted at a faster rate due to decays and losses. The stored and lost muon populations have different momentum distributions, and so the different rate of depletion creates a time-dependent average of the muon momenta: d p /dt = 0. The spin-momentum correlation will combine with the time-dependent muon losses to produce a time-dependent phase: The three subsections that follow address, in turn, the following topics: (1) the data-driven determination of the absolute rate of muon losses during a fill, (2) the dataand simulation-driven determination of dϕ 0 /d p at the fit start time, and (3) the data-and simulation-driven determination of the d p /dt during a fill. With these rates and correlations in hand, one can evaluate the impact on the muon loss correction factor C ml .
FIG. 13. The intensity distribution of where muons first strike the r0 = 45 mm radius collimators (black circle). For any particular muon, this will occur when its horizontal and vertical betatron oscillations conspire such that the transverse displacement is at r0 and, at the same time, it is at the azimuthal location in the storage ring of a collimator. Some muons, which will eventually be lost, can survive hundreds to thousands of turns before this condition is met.

A. Muon loss rate determination
Muons that exit the storage ring during the 30−650 µs measuring period deplete the population faster than the expected time-dilated decay e −t/γτµ . The shape of the muon loss function L(t) is accurately measured from the rate of coincident signals in three consecutive calorimeter stations, where each station records an energy deposit of ≈ 170 MeV and the time between stations is ≈ 6.4 ns, corresponding to the energy deposit of a minimum ionizing muon and its time of flight between stations, respectively. The muon precession fit includes a term that multiplies the overall normalization N 0 such that The scale parameter K loss can be accurately extracted from the precession analysis to provide the absolute scale of the muon loss function [11]; this is necessary to estimate the phase-related systematic error. We note that the loss function is a property of the beam and should, therefore, be the same for all calorimeters and energy bins in the precession analyses. An important measure of the rigor of this method is that K loss is independent of which calorimeters are used or which energy bins are selected in the precession fits.
FIG. 14. The integrated fractional muon loss versus time in fill from 30 µs following muon injection. The four curves are from the different run groups. The smaller loss fraction curves are from the n = 0.120 datasets (1b and 1c) and the greater loss fractions are from the n = 0.108 datasets (1a and 1d). The absolute scale is determined from the K loss parameter following final precession frequency fits. The uncertainty bands on the curves come from two different precession frequency analyses and whether a small empirical slow correction term to ensure stability of K loss versus energy is included. Figure 14 shows the accumulated loss fraction (f loss ) for the four datasets in Run-1, defined for t > t s as This gives the fraction of muons that have been lost from the storage ring with respect to the number present at the fit start time. Although all curves rise steeply at early times and gradually at later times, it is clear there are two distinct groups, which are associated with the two different tune n values (see Table I). We note that the loss fraction in this figure is approximately 8 times larger than the loss fraction measured for Run-2, when the damaged resistors had been replaced. This is strong evidence that the dynamic beam motion during storage led to a high degree of scraping on the collimators at early times, until the beam relaxed to its nominal central value when the voltages had stabilized.

B. Phase-momentum correlation determination
Nonzero spin-momentum correlations dϕ 0 /d p are generated in the Muon Campus beam line and muon storage ring. These arise primarily during the circulation of the muons around the Delivery Ring (DR), which is composed of FODO cells and bending dipole magnets. In particular, a dipole bending magnet will change the angle between the muon spin and momentum by ∆ϕ ≈ a µ γη b , where η b is the angle at which the muon bends through the dipole field and the momentum dependence is embed-ded in the γ factor. For k full revolutions of the DR, the angle between the spin and momentum will advance by ∆ϕ ≈ 2πa µ γk, where the sign of phase is defined in the sense that the spin angle precesses with the functional form cos(ωt + ϕ). For a hypothetical muon distribution entering the DR with no phase-momentum correlation, the four turns around the DR will imprint a change in ϕ of 8.6 mrad per 1% of ∆p/p 0 on the overall distribution.
A complete end-to-end simulation has been performed to determine the muon distribution phase space at the exit of the inflector from muons born in all distinct target and beam line regions. The simulation tracks spin and kinematic variables from the production in the target to the delivered beam. A plot of the average phasemomentum correlation from this simulation is shown as the blue band in Fig. 15.
It is possible that this spin-momentum distribution is further perturbed during the storage ring kick, because, as illustrated in Fig. 1, the kick does not apply an equal impulse to muons that are distributed longitudinally throughout the incoming bunch. Although this effect appears to introduce a negligible spin-momentum correlation in the simulation, it was possible to perform a direct measurement of the correlation that exists during the measuring period and compare it to the simulation as shown in Fig. 15.
Three special runs were made with the magnetic field of the storage ring set at nominal (1.45 T), reduced (-0.68%), and increased (+0.67%) values. At each setting, an approximately ±0.15% momentum slice of the broad incoming beam is stored, with its central value corresponding to the momentum of the adjusted field settings. The statistical precision on a few hours of beam is sufficient to determine the average spin phase at injection, the precession frequency, and the time-dilated muon lifetime. The values are obtained from fitting the positron versus time plots to Eq. 3. The change to the precession frequency is proportional to the magnetic field values and readily serves to determine the actual field (and by proxy, momentum) setting. The black points in Fig. 15 show the results of these direct measurements. The fitted slope of (−10.0 ± 1.6) mrad/(%∆p/p 0 ) agrees well with the simulation. The error quoted is from the fit alone.

C. Lost-muon momentum correlation determination
A set of special measurements was made to determine the behavior of the muon losses as a function of time in fill and momentum. The two damaged resistors were reinserted into their Run-1 locations during a short period of systematic tests at the beginning of Run-3. One 8-h period of data collection was acquired in otherwise nominal conditions, to reestablish the Run-1 time dependence of the CBO frequency (see Sec. VI A) and to provide data that could be used in a COSY simulation of the storage ring behavior under these conditions; see below. The FIG. 15. Phase-momentum correlation from an end-to-end simulation (blue band) and from a data-driven approach (black). The simulation gives the result at the entrance to the storage ring. The three data points are obtained by fits to muon precession frequency data at nominal, reduced, and increased central magnetic field values. The phase reported for these data represented muons stored and fit during the measuring period. The phase dependence on momentum from the data is −10.0 ± 1.6 mrad/%∆p/p0. Delivery Ring momentum collimators were used to bias the incoming muon momentum distribution. The collimators can be driven separately on both the high-and low-momentum sides and can traverse the entire horizontal width of the beam. Figure 16 shows horizontal stored-muon distributions F i (x) from a subset of these special runs. The dashed line corresponds to the nominal, full-acceptance distribution used for normal data taking. The three colored distributions are from runs where the low, high, or both momentum collimators were used to bias the stored momentum distribution. The fractions in the legend indicate the intensity with respect to the nominal case. For example, a collimator would be moved until the muon storage rate was reduced to that fraction. The muon loss rate versus time in fill was measured for each collimator setting, providing eight distributions (not all shown in the figure) that were used in the analysis below.
An example of the two extremes-1/5-low (blue) and 1/5-high (red)-correspond to the relative muon loss rates shown in Fig. 17. The low-momentum muons are lost disproportionately early in the fill and the highmomentum muons are lost more often at later times; the other distributions provide intermediate cases.
The eight distributions were used to parameterize an analytical loss-rate function, whose form was motivated by simple simulation phase space studies. Various models were used to assess the reliability of the conclusions from input assumptions. The solid gray curve in Fig. 16 represents a loss rate probability function determined from the integrals of the muon loss versus time distributions from 30 to 70 µs. This is simply meant to be illustrative, as the function evolves in shape throughout the fill. What is needed is a time-dependent muon loss probability function l(x, t), which can be applied to the nominal momentum distribution to yield the time dependence of the average momentum of the stored muons. That information is readily translated into a time-dependent spin phase of the stored muons, ϕ(t). The parameters of the function l(x, t) are determined by fitting the eight distributions over increasingly long time ranges from the fit start time t s = 30 µs. At each time t, the fit is per-formed using eight equations of the form where i runs from 1 to 8 for each of the special runs and R min and R max represent the minimum and maximum, respectively, of the possible radii of the stored muons; F i (x) is the measured intensity of the fast rotation distribution for that run as a function of the equilibrium radius, normalized to 1; and γτ µ ≈ 64.4 µs is the timedilated muon lifetime. The empirical loss function l(x, t), which depends on time and radius, is determined in the fit. The measured integrated triples spectrum L i (t ) is integrated from the fit start time to time t. The extra term e t /γτµ is included in order to follow the convention of the muon loss term in the decay positron fit function, as expressed in Eq. 20. The normalization by H i , the total number of positrons measured in that dataset, ensures that the eight special runs can be correctly compared. An analytic form of l(x, t) is assumed in order to perform the fit. From simulations, l(x, t) is expected to peak near the edges of the storage distribution, as muons that have a high or low equilibrium radius are more likely to be lost. Several forms of l(x, t) with this qualitative behavior were compared, including a sum of two Gaussians and a piecewise sum of two parabolas. The same analytic form was used throughout the fill, and the fit parameters were allowed to vary with time to account for the changing behavior of the lost muons.
At each time, the remaining stored distribution F curr (x, t) is calculated using the equation where F 0 (x) is the fast rotation distribution (see Sec. III) of the full physics dataset, normalized to 1, and represents the radial distribution of the stored muons at the fit start. The second term represents the total number of muons that have been lost up to that time, scaled by the fractional loss correction f loss (t), which emerges from the decay positron fit, as seen in Fig. 14.
The average radius of the stored distribution is then extracted at each time from F curr (x, t). This average radius can be converted into momentum units using Eq. 15. The average ∆p/p 0 of the stored distribution is converted to ϕ using the measured phase-momentum correlation of (−10 ± 1.6) mrad/(%∆p/p 0 ). The ϕ(t) determined by this method for the four Run-1 datasets, using three different analytical forms of l(x), is shown in Fig. 18; ϕ(t) behaves very similarly, regardless of the form of l(x) used.

D. Value of the muon loss correction
The phase angle ϕ(t) is parameterized using a highdegree polynomial function, which is used to generate a  set of points based on the five-parameter decay positron fit function Eq. 4 with the empirical ϕ(t) inserted. A fit is then performed on the generated points using ϕ(t) = ϕ 0 , as is assumed in the physics analysis, with all parameters allowed to float. The bias to ω m a is extracted by comparing the input ω a value to the value extracted from fitting with a constant ϕ. This analysis is repeated for all four datasets.
The correction to ω m a from this effect, C ml , is given in Table VII. The muon loss-induced phase change artificially increases the measured value of ω m a , so C ml is negative. It varies from −17 to −3 ppb depending on the dataset, with Run-1a and Run-1d having larger shifts because of their higher muon loss rates. Three sources of uncertainty are common to all run groups, and they are added linearly. The measured phase-momentum correlation of (−10 ± 1.6) mrad/(%∆p/p 0 ) contributes a 16% uncertainty (1 − 3 ppb). The choice of analytical form of l(x, t) contributes a 0 − 2 ppb uncertainty. The use of different f loss (t) functions from different analyses contributes 1 − 2 ppb. To quantify the latter two sources of uncertainty, the full analysis described in the previous section was repeated for every combination of l(x, t) and f loss function. This procedure yields σ(C ml ) = 2 − 6 ppb, which is small on the scale of other Run-1 uncertainties.

VIII. CORRECTION FOR MUON DISTRIBUTION TIME DEPENDENCE Cpa
The phase of the ω a oscillation at the moment of a muon's decay is related to the orientation of the muon spin vector relative to its momentum at injection into the storage ring. The sign of the phase is defined by the convention in Eq. 4 that the positron intensity spectrum modulation is described by a term proportional to A cos(ω a t + ϕ). As discussed in Sec. VII B, the average phase of the incoming polarized muon beam is determined by upstream beam line components and the number of turns in the DR; its value does not affect the extraction of ω m a . Ultimately, the observed phase is what a calorimeter detects, an integration of decays from all locations that produce a positron signal in an energy bin E. The phase of the fitted distribution corresponds to the injection phase ϕ 0 plus any average orientation of the muon spin with respect to its momentum that maximizes the anomalous precession signal.
In Sec. VI A, we described how the calorimeter acceptance is dependent on the transverse decay coordinate. An azimuthally averaged transverse distribution is shown in Fig. 19. This static image does not convey the radial and vertical oscillations of the mean and width at frequencies associated with the betatron oscillations. They require modification to the normalization, asymmetry, and phase terms in Eq. 4 to fit the spectrum to obtain ω m a . The form of these modifications is discussed in Refs. [3,11], and an example form used in our Run-1 ω a analysis is shown in Appendix D. The phase for a given (x, y) decay coordinate depends on the orientation of the muon's spin that maximizes the acceptance. The 24 calorimeters are finite sized, placed to the inside of the muon trajectory, and positioned at uniformly spaced azimuthal locations. Therefore, the spin orientation of a muon that maximizes acceptance into this system is not parallel to its momentum but rotated slightly radially inward. This rotation, captured by an effective phase shift ϕ pa , is a function of transverse decay coordinate (x, y) because of acceptance effects. Figure 20 is a "phase map" averaged over azimuth and weighted by the asymmetry method used to extract ω m a from the positron intensity time spectra. The map varies more strongly in the vertical direction and less so radially. The procedure to create phase maps using GM2RINGSIM is described in Appendix C. Briefly, muon decays are generated over all (x, y, φ) coordinates with full muon spin precession versus time in fill properly included. Individual positron intensity spectra from decays originating in a matrix of (x, y) transverse bins are fit to determine ω a and ϕ xy . The procedure naturally includes the effects of acceptance, and, indeed, acceptance and asymmetry maps are also obtained for each calorimeter.
The spectrum summed across all detectors and all decay coordinates has a modulation frequency ω a and a net, spatially averaged phase (ϕ 0 +φ pa ). As long as the muon distribution remains constant throughout the measurement period, the net phase is constant. The dependence of the decay-coordinate phase ϕ pa (x, y) on detector acceptance was understood well enough from the E821 experiment to help guide our voltage stability specifications in the development of the new ESQ power-supply network. However, as explained in Sec. VI, the two damaged resistors in Run-1 spoiled the voltage stability requirement on the Q1L upper and lower electrodes, which, in turn, led to a time dependence of the mean and width of the muon distribution. In this section, we evaluate how this effect results in a time dependence of the average coordinate-dependent phase contribution that is not included in the fits to extract ω m a from the positron spectrum. A correction factor C pa is needed to correct ω m a to the true ω a needed to determine a µ . Figure 21 illustrates how the time dependence of the muon distribution vertical and horizontal means and widths can lead to an average phase shiftφ pa → ϕ pa (t). The exaggerated Gaussian profiles are representative of how such a muon distribution would evolve from early (dotted red line) to late (dashed blue line) times in a fill. The average phaseφ pa at any time is calculated by taking a weighted average of the phase values in the map projection, with the weights from the beam intensity (assuming uniform asymmetry and acceptance). The largest effect is from the reduction in the vertical width, whereφ pa is evidently different early compared to late in the fill. In contrast, a small shift of the vertical mean both gains and loses phase nearly symmetrically. The radial phase projection is nearly linear, and, therefore, a mean shift will cause a variation in the phase, whereas a width change is relatively balanced. The radial effects are relatively small, and it is the coherent reduction in the vertical width (see Fig. 11b) that dominates the correction to ω m a and must be evaluated. An additional effect that contributes to the phaseacceptance correction has its origin in CBO decoherence. Early in a fill, the betatron oscillations are largely in phase and the beam moves coherently back and forth in the radial direction. When the beam has fully decohered, all calorimeters sample the full radial distribution at once, and so the relative acceptance of the full aperture determines the average phase. Conversely, early in the fill when the oscillations are still coherent, the calorimeters sample muon decays from only a subset of the radial distribution at any particular time. In this scenario, the relative acceptance between different radial positions is not important. As a result of this difference between the early and later stages of CBO decoherence, the average phase for an individual calorimeter drifts from early-to-late times in the fill with an impact as large as O(100) ppb. However, detectors on opposite sides of the ring see the CBO oscillations out of phase with one another, so when considering the sum of the calorimeter data there is a strong cancellation with no net contribution to the C pa correction at a statistically significant level. We do assign an uncertainty owing to the imperfect cancellation of detector pairs that are π radians out of CBO phase. We note that a similar reduction is realized when comparing individual calorimeter positron decay spectra fits to those of the sum. The terms that describe the CBO oscillation amplitudes are 6−7 times smaller for the summed spectra compared to fits to individual calorimeters, and the reduction would be more complete if the detector acceptances were identical.
The strategy to evaluate the net phase shift to each of the four datasets in Run-1 involves the following steps: (1) Create high-fidelity maps of acceptance, asymmetry, and phase for each calorimeter (see Appendix C). namics 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, ϕ c pa (t). (5) Generate Monte Carlo data using the full fit function, including the predicted behavior of the phase, ϕ c pa (t). Fit these pseudodata using the normal fit function to determine the difference in the extracted ω m a compared to the input. This difference yields the correction factor C pa for each calorimeter.
A. 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 D. 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 ω CBO , ω 2CBO , ω VW , ω y , and ω m a are all removed. This randomization procedure also has the effect of removing the CBO decoherence phaseacceptance 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 y and width y 2 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.

B. 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 ω m a 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 y (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°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 pred (φ, t) = y tkr (t) y rms pred (φ, t) y rms pred (φ tkr , t) .
The width ratio is evaluated directly from virtual tracking planes in GM2RINGSIM and via the ratio of beta functions in COSY ( β y (φ, t)/β y (φ 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 closedorbit distortions due to azimuthal variation of the main dipole field [10]. Each entry from the measured distribution is scaled as where the predicted values for x rms (φ, t) andx(φ, t) are taken from tracking planes in GM2RINGSIM and calculated in COSY using β x (φ, t), D x (φ, t), and the reconstructed distribution of muon momenta. Measurements of the main dipole storage field are used to calculate the closed-orbit distortion included inx pred (φ, 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.

C. 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 ϕ c pa (t) using the following weighted sum: .
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 ε c , A c , and ϕ c pa , respectively. An example evaluation of ϕ 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.
The size of the mismeasurement of ω m a as a result of the time-dependent phase ϕ c pa (t) is estimated using simulated data. A histogram is generated for each calorimeter according to the full fit function used to extract ω m a , 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 ϕ c pa (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 ω m a that were input and extracted from the fit give the value of the correction C 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 ω m a is larger than the true value. The correction C pa must reduce the measured value and is, therefore, negative. The C pa corrections are shown in Table VIII. 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 µs, which reduces the correction by a factor of 1.8.
The evaluations of statistical and systematic uncertainty are also detailed in Table VIII. 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 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 CBOmodulated 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 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 ω m a fits with predictions of ϕ c pa 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 δC 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 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 ω m a (Ref. [11]) to its true physical value ω 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 storedmuon 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 ∼100 µ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 µs is used to limit the systematic shift to the phase.
A. Summary of Run-1 net corrections to ω m a Table IX 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 ω m a from the precession fits to the true ω a value that is needed to determine a µ 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, whereB represents the muon-weighted average magnetic field that is discussed in Ref. [10]. to ω m a from the four Run-1 datasets. These values are computed with the full correlation matrix formalism used to provide the measured value for ωa. The values, thus, reflect fully correlated systematics and weighting from the statistical uncertainties from each dataset. sextupole (∼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 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 φ, the simplest possible solution is wherek is a constant. The electric field is given by and expanding in x/R 0 : where k = −2k/R 2 0 . 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 V as part of the ESQ calibration.

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 where x e = ηδ = η ∆p p and the simple solution above to Laplace's equation is assumed. The average electric field along the trajectory is where L is the length of the trajectory. Using where k = n vsB R0 . 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. A8 into Eq. 13 gives the correction to ω a due to fractional momentum offset δ and betatron amplitude x β0 as The next step requires averaging C e (δ, x β0 ) over the entire momentum and CBO distribution: where δ = x e /η. If x e /η = 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 In fact, momentum offset and betatron amplitude are strongly anticorrelated, which will be discussed below. If x e = 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 For the equilibrium radial distributions typical of Run-1 (see Fig. 6b), Conservatively estimating that x 2 β0 < (45mm) 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%.

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 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. A5. The ratio of the coefficients of the quadratic and linear terms is r hyp = − 1 2R0 = −0.0703m −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 f it = b2 b1r0 = −2.71281 × 10 3 /(1.01609 × 10 6 × 0.045m) = −0.0593m −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.

Appendix B: E-field and Pitch Corrections Explored Numerically
In simulation, ω 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 ω 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 Efield 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 where the experiment actually measures the scalar quantity ω m a and not its vector approximation. However, Eq. B1 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 β × E electric field term. Define Φ a = T 0 ω a dt, there T is the total integrated time that the muon is precessing. In the absence of E, Φ a can be computed by integrating B ⊥ along the muon trajectory. The contribution from the electric field to φ a is given by where ∆p = p−p 0 . In the tracking code, it is straightforward to compute the sums in Eq. B2. The electric field correction for the trajectory, assuming that B is parallel to β × E, is In a similar manner, the expressions for the pitch correction C p developed in Sec. V can be rearranged to read To test the equivalence of the shift in precession frequency based on Eqs. B3 and B4, both were computed simultaneously along a single trajectory. The muon is propagated for 1000 2π/ω a in order to determine the spin tune by integrating the Thomas-BMT equation with sufficient accuracy, and β × E is summed along that same trajectory. The trajectories chosen to test Eq. B3 have momentum offsets but no vertical motion. To test Eq. B4, 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.
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,  FIG. 25. Top: contribution to ωa due to the electric field as computed by spin tracking (integration of the Thomas-BMT equation) and integration along trajectories (Eq. B3) 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 ωa due to the vertical oscillation computed with spin tracking and integration along the trajectory (Eq. B4) for muons with magic momentum and zero radial betatron amplitude. Note that the formula that assumes linearity overestimates the contribution at large amplitudes.
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 ω a . The as-built ESQ plate alignment is known with a precision better than ±1 mm, and the high-voltage set points to ±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 σ(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 cor- rection 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 σ(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, φ). 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 ω a frequency, which is known to sufficient precision from the E821 experiment.
The simulated time and energy distributions of the decay positron match well with the experimental measurement as shown in Fig. 26. This 2D intensity mapwith the decaying exponential removed-shows positron energy versus time in fill from 30 to 80 µs. The clear modulations are at the frequency ω 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 ω a analysis methods with different positron weightings (see Ref. [11]). In practice, the asymmetry-weighted method provides the most precise determination of ω 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.

Acceptance maps
The average acceptanceε ij for a muon decaying in transverse position bin ij is given bȳ where the first sum is over positron truth energy bins k, the second sum is over positron deposited energy bins k above a threshold k th , f k is the fraction of events with energy k, ε ijk is the geometric acceptance for muon decays in bin ij and positron energy k, R ijkk is the calorimeter response function, and w k is the analysis weighting factor [11]. R ijkk represents the probability that a positron in position bin ij and truth energy bin k ends up in the k detected energy bin. This response function is illustrated in Fig. 27.  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.

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/τ (1+A cos (ω a t + ϕ)), where ϕ is the fitted phase. The fitted phase for all decay positrons is denoted as ϕ 0 , and the fitted phase for each decay xy bin and energy bin is denoted as ϕ ijk . The value of each ijk bin in the phase map is then given by ϕ ijk − ϕ 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 phaseacceptance effect than those behind the tracker stations.
The parameters of the form A N,x,i,j and ϕ N,x,i,j correspond to the effect of the ith moment of the radial (x) beam distribution at the jth multiple of the relevant fundamental frequency (in this case, ω CBO ) on the rate normalization N . Analogous parameters model the modulation of the average asymmetry A and phase ϕ, as well as the effect of moments of the vertical (y) beam distribution.  FIG. 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.