A Precision Measurement of the Muon Decay Parameters $\rho$ and $\delta$

The TWIST collaboration has performed new measurements of two of the parameters that describe muon decay: $\rho$, which governs the shape of the overall momentum spectrum, and $\delta$, which governs the momentum dependence of the parity-violating decay asymmetry. This analysis gives the results $\rho=0.75014\pm 0.00017(\text{stat})\pm 0.00044(\text{syst})\pm 0.00011(\eta)$, where the last uncertainty arises from the correlation between $\rho$ and the decay parameter $\eta$, and $\delta = 0.75067\pm 0.00030(\text{stat})\pm 0.00067(\text{syst})$. These are consistent with the value of 3/4 given for both parameters in the Standard Model of particle physics, and are a factor of two more precise than the measurements previously published by TWIST. A new global analysis of all available muon decay data incorporating these results is presented. Improved lower and upper limits on the decay parameter $P_\mu^\pi\xi$ of $0.99524<P_\mu^\pi\xi \leq \xi<1.00091$ at 90% confidence are determined, where $P_\mu^\pi$ is the polarization of the muon when it is created during pion decay, and $\xi$ governs the muon decay asymmetry. These results set new model-independent constraints on the possible weak interactions of right-handed particles. Specific implications for left-right symmetric models are discussed.


I. INTRODUCTION
In the Standard Model (SM) of particle physics, the charged-current weak interaction violates parity maximally-only left-handed particles (or right-handed antiparticles) are affected. The TWIST experiment is a high-precision search for contributions from non-SM forms of the charged-current weak interaction, including parity-conserving currents.
The decay of the positive muon into a positron and two neutrinos, µ + → e + ν eνµ , is a purely leptonic process, making it an excellent system for high-precision studies of the weak interaction. It proceeds through the charged weak current-mediated by the W boson-and can be described to a good approximation as a four-fermion point interaction. The matrix element for the most general Lorentz-invariant, local, four-fermion description of g γ ǫµ ψ eǫ Γ γ ψ νe ψ νµ Γ γ ψ µµ , (1) where the g γ ǫµ specify the scalar, vector, and tensor couplings between µ-handed muons and ǫ-handed positrons [1], and satisfy certain normalizations and constraints. In the Standard Model, g V LL = 1 and all other coupling constants are zero. The probability Q ǫµ (ǫ, µ = L, R) for the decay of a µ-handed muon into an ǫ-handed positron is given by where δ ǫµ = 1 for ǫ = µ and δ ǫµ = 0 for ǫ = µ. The probability: (3) sets a model independent limit on any muon right-handed couplings [1,2].
The differential muon decay spectrum [3,4,5] can be described following the notation of Fetscher and Gerber [6] as where G F is the Fermi coupling constant, θ s is the angle between the muon spin and the positron momentum, E max ≈ 52.8 MeV is the kinematic maximum positron energy, x = E e /E max is the positron's reduced energy, x 0 = m e /E max is the minimum possible value of x, corresponding to a positron at rest, and P µ is the degree of muon polarization at the time of decay. P µ can be used to determine P π µ , the degree of muon polarization at its creation from pion decay, when the amount of depolarization undergone by the muon is known. The two components of the decay spectrum in Eqn. (4) are the isotropic component: (5) and the anisotropic component: R.C. represents the electromagnetic radiative corrections, which have been calculated to O(α 2 ). Corrections due to the strong interaction in loops give a fractional contribution on the order of 4 × 10 −7 [7], which is more than two orders of magnitude smaller than the ultimate precision goals of TWIST. The quantities ρ, δ, η, and ξ, often called the Michel parameters, are bilinear combinations of the weak coupling constants and describe the shape of the decay spectrum. These decay parameters can be used in combination with other muon decay measurements, such as inverse muon decay (e − ν µ → µ − ν e ) and the polarization of the decay positron, to determine limits on the weak coupling constants. Left-Right Symmetric (LRS) models [8] comprise an interesting class of extensions to the SM. These models include a right-handed weak coupling, which is suppressed by the mass of the associated gauge boson. LRS models contain four charged gauge bosons (W ± 1 , W ± 2 ) with masses m 1 and m 2 , and two additional massive neutral gauge bosons. The mass eigenstates W 1 and W 2 are related to the weak eigenstates W L and W R through a mixing angle ζ. g L and g R are the coupling strengths of the LRS weak interaction to left-and right-handed particles. LRS models affect the muon decay spectrum, including a modification to the decay parameter ρ: therefore requiring ρ ≤ 0.75. To a good approximation, the value of δ is unaffected by LRS models. Many other proposed SM extensions also lead to modifications of the Michel parameters. For example, supersymmetric models can lead to a non-zero value for g S RR [9]. The values of ξ and δ can be combined to provide a model-independent limit on right-handed couplings in muon decay: Previous measurements by TWIST had set new limits, of ρ = 0.75080 ± 0.00032(stat.) ± 0.00097(syst.) ± 0.00023(η) [10] and δ = 0.74964 ± 0.00066(stat.) ± 0.00112(syst.) [11]. This paper presents new results from a refined analysis of newer data, providing details about the TWIST experiment and analysis (Sects. II and IV) including improvements over the previous studies (Sect. V), discusses the validation of the Monte Carlo simulations (Sect. VI), describes the current state of the systematic uncertainties in the experiment (Sect. VII), and presents new measurements of ρ and δ (Sect. VIII). These new results are then incorporated into a global analysis of all muon decay data (Sect. IX).

A. Experimental Setup
Highly polarized muons resulted from the decay of pions stopping at the surface of a carbon production target bombarded by 500 MeV protons at TRIUMF. The M13 beamline [12] selected positive muons with a momentum of 29.6 MeV/c, with a momentum bite of ∆p/p ≈ 1% FWHM, and delivered these in vacuum to the TWIST spectrometer [13] at a typical rate of 2.5×10 3 per second. The beam had a contamination of positrons, at the same momentum, with a typical rate of 22 kHz, as well as a small fraction of pions.
The muon beam was characterized and tuned using a low pressure (8 kPa dimethylether (DME) gas) removable beam monitoring chamber located upstream of the TWIST spectrometer [14]. The beam monitor was inserted for measurement of the beam properties, and removed during data-taking.
The TWIST spectrometer [13] was designed to measure a broad range of the muon decay spectrum, allowing the simultaneous determination of the decay parameters. The detector consisted of 44 drift chambers (DCs) and 12 multiwire proportional chambers (PCs) in a planar geometry, symmetric about a muon stopping target foil at the centre. Figure 1 shows the upstream half of the detector and the four PCs surrounding the stopping target. The chambers were placed in a highly uniform 2 T solenoidal magnetic field. The z axis was defined to be the detector axis. To reduce scattering and to allow the muons to reach the central target foil, the detector was designed to be very thin; there was approximately 140 mg/cm 2 of material from the vacuum of the M13 beamline through to the centre of the stopping target.
The muon stopping target was a 71 ± 1 µm foil of 99.999% pure aluminum, which also served as a cathode foil for the adjacent two PCs (see below). The muon stopping distribution, as determined from the last DC or PC plane in which the muon left a signal, was used in a feedback loop to control the fractions of He and CO 2 in a gas degrader at the end of the beamline, in order to maintain the average stopping position of the selected muons at the centre of the target. Muons were required to be recorded by the PC immediately before the stopping target (PC 6) and not by the PC immediately after the target (PC 7); simulations showed that 97.0 ± 1.5% of selected muons stopped in the target, with the rest stopping in the CF 4 /isobutane, the cathode foils, or the wires in the vicinity of the target. The experiment watched for muons decaying at rest into positrons; decays of muons in flight were identified and discarded.
Each PC and DC consisted of a wire plane and two cathode foil planes, all oriented perpendicular to the detector axis. Chambers were grouped together in modules as described below, and the cathode foils interior to each module were shared between adjacent chambers. Each wire plane was in either U or V orientation, and these were at right angles to each other and 45 • to the vertical.
Four PCs were located at each end of the detector, and two on either side of the stopping target. Each PC included 160 sense wires at 2 mm pitch. A fast drift gas (CF 4 /isobutane) was used in these chambers. The PCs mainly provided particle timing information; in addition the widths of the timing signals from the PCs were used to discriminate muons from positrons.
Each half of the detector included 22 DCs, with a slow drift gas (DME), which had a small Lorentz angle; these provided precise measurements of the e + position as it crossed the chamber. Each DC plane [15] included 80 sense wires at 4 mm pitch. The DC planes were very slightly asymmetric in the direction along the detector axis, in that the wires were located 150 µm off of the centre of the cell. 28 of the DCs were paired together in 14 modules, with U and V wire planes in each pair. Additional "dense" stacks of eight planes each were located near either end of the detector.
The space between the chambers was filled with a 3% N 2 and 97% He gas, minimizing the material thickness of the detector. The chamber gas and the helium mixture were maintained at atmospheric pressure. The differential pressure between the chambers and the surrounding volumes was controlled to stabilize the positions of the chamber foils.
Chamber alignments perpendicular to the detector axis were measured using 120 MeV/c pion tracks with the magnetic field off. Alignment measurements were taken at the beginning and end of the data-taking period. These data were also used to measure wire time offsets introduced by the electronics and cabling. The plane positions and rotations about the beam direction were determined to an accuracy of 10 µm and 0.4 mrad. Relative wire positions were known to 3 µm. The alignment of the chambers to the magnetic field was measured using positron tracks in the 2 T field, and had an uncertainty of 0.03 mrad.
The chambers were positioned within the 1 m bore of a liquid-helium-cooled superconducting solenoid. The magnet was placed inside a cube-shaped steel yoke, approximately 3 m per side, designed to increase the uniformity of the field within the detector region. The shape of the z component of the magnetic field was mapped using an array of Hall probes mounted on a radial arm, aligned along the solenoid axis and calibrated using an NMR probe. The probes' positions at each field sample were known to ±1 mm, and the field was mapped to a precision of ±1 × 10 −4 T. The standard operating setting for TWIST is 2 T at the centre of the solenoid; at this setting the field is uniform to within 8×10 −3 T within the tracking region. A finite-element simulation of the magnetic field was performed using Opera3d [16], providing the full three-dimensional magnetic field throughout the detector volume and well outside the measured re-gion. Within the tracking region the simulated field map agreed with the measurements to within ±2 × 10 −4 T. Event acquisition was triggered by a thin scintillator upstream of the spectrometer. Only 1.8% of events were triggered by a beam positron, in spite of the much higher rate of beam positrons compared to muons; most beam positrons were discriminated against by the trigger threshold. Remaining beam positrons were identified with high efficiency in the data since they left signals in the full length of the detector.
Preamplifier chips were mounted directly on the chambers and drove custom postamplifier and discriminator modules placed 2 m away. The chamber behavior was exceptionally stable, with less than one high voltage trip per month. All output channels were functional.
Data acquisition from scintillators and tracking chambers [17] was performed using LeCroy Model 1877 TDCs. The trigger and TDC read-out recorded signals in 0.5 ns time bins from 6 µs before to 10 µs after a muon passed through the trigger scintillator. The start and stop of each pulse was recorded so that the pulse width could be determined; these widths can be related to the amount of energy deposited in the cell. In this configuration up to eight wire signals could be recorded for each wire in any triggered event. A fixed blanking time of 80 µs was imposed after each accepted trigger, to allow each TDC to finish conversion before the next event was recorded.
The data considered for this analysis were more than 1.5×10 9 muon decays, taken during 2004. The data sets, summarized in Table I, were taken under a variety of conditions (low polarization from beam steering, rate, muon stopping position, etc.). Provided the simulation reproduced these conditions correctly, the decay parameters extracted from the data should be independent of the run conditions. The decay parameters were extracted from the data by comparing the data to a simulated spectrum, as discussed in Sect. IV D. Monte Carlo (MC) simulations were run to match the conditions of each of the seven main data sets, in addition to simulations run for studying systematic uncertainties. Each simulation included 2-3 times as many muon decays as the corresponding data set. The simulations were implemented using Geant 3.21. The output was in exactly the same format as the files produced by the data acquisition system, and were processed in the same way as real data. Space-time relationships (STR) determined with garfield [18] were used to model drift chamber response.
The theoretical decay spectrum included full radiative corrections at O(α) [19], as well as O(α 2 L 2 ) and O(α 2 L) [20,21], where L = log(m 2 µ /m 2 e ) ≈ 10.66. O(α 2 L 0 ) terms have also been calculated [22]; the effect of neglecting these last terms has been evaluated and shown to be negligible (see Sect. VII F). All radiative corrections are calculated within the Standard Model. The values of the muon decay parameters assumed by the simulation in its theoretical decay spectrum were kept hidden until the end of the study; see Sect. IV D below for details.
The simulation included all known depolarization effects, including interactions with the magnetic field as the muon enters the solenoid and depolarization in the stopping target.

IV. DATA ANALYSIS
The full analysis procedure was applied in the same way to both data and simulation. To the level that the simulation accurately represented the data, this canceled spectrum distortions due to detector response, positron energy loss and scattering, reconstruction biases, and other effects, which would otherwise lead to systematic errors in the measurement.
Due to the large amount of simulation and analysis required, the Western Canada Research Grid (WestGrid) was used. TWIST used approximately 10,000 CPU-days for this simulation and analysis.

A. Event Reconstruction
To reconstruct an event, the hits-signal times on individual wires-were first grouped based on timing information from the PCs. Tracks within these groups were then identified using the distribution of DC hits in space and time, as well as the hit widths. DC hits associated with decay positrons were then used to reconstruct the energy and angle of each positron.
Pattern recognition was performed on the decay positrons using the spatial hit distributions, to determine an initial estimate of the positron track. The track fit parameters were the position and momentum three-vectors at the DC closest to the target; the main parameters of interest were the positron's total momentum p, and the angle θ between the positron momentum vector and the detector axis. Multiple overlapping tracks could be distinguished at this stage. The initial track estimate was then refined in a χ 2 fit using the hits' drift times, in combination with maps of the STRs of the drift chamber cell as provided by the garfield chamber simulation software [18]. The positron track was assumed to deviate from a helix by continuous energy loss through the gas volumes and discrete energy loss through each foil encountered, using mean energy loss formulas [6]. Tracks were allowed discrete deflections at each DC module and in the dense stacks. The deflection angles were fit parameters, with associated penalties to the fit χ 2 , based on the method described by Lutz [23].
Event and track selection cuts were then applied and the decay spectrum was assembled. A muon hit was required in PC6, immediately upstream of the stopping target, to ensure the muon reached the target. A hit in PC7, immediately downstream, vetoed the event. Events with multiple muons were rejected. The muon was required to decay between 1 µs and 9 µs after the trigger; the delay ensured that the ionization from the incident muon was collected before a hit from the positron was recorded. Events where a beam positron arrived within 1 µs of either the muon or decay were rejected as well. Additional cuts included the muon flight time through the M13 beam line, used to reject muons from pion decays in flight, and a requirement that the muon stopped within 2.5 cm of the detector axis, which ensured that all decay positron tracks within the fiducial region (Sect. IV C) were fully contained within the detector. Note that these cuts depend on observations of the muon prior to decay or of the time of decay, and are hence unbiased in terms of the properties of the positron. Figure 2 shows a muon decay spectrum reconstructed from the "Centred Stops" data set (3×10 8 muon decays; see Table I).

B. Energy Calibration
Differences in the details of the muon stopping position within the target foil, differences in the STRs, the accuracy of the magnetic field maps, and knowledge of the wire positions resulted in differences in the energy scales of the simulation and the data. The kinematic endpoint region (near 52.3 < p < 53.4 MeV/c) was used to determine the relative energy calibration between the data and the simulation after reconstruction was complete. The relative positions of the endpoint were compared between data and simulation in bins of angle (see Fig. 3 for an example), and the differences were parameterized separately for upstream and downstream parts of the spectrum, according to where ∆p e is the difference in the momentum of the spectrum endpoint between data and simulation. Equation 9 describes an energy loss proportional to the amount of material encountered by a particle; for a planar detector the amount of material is proportional to 1/ cos θ. This form was found to describe ∆p e well, also. Uncertainties in the calibration parameters B i and A i were statistical. This relative calibration technique replaced the calibration using an analytic approximation to the end point shape used in the previous analyses [10,11,24].

C. Fiducial Region
Limitations of the reconstruction and physical aspects of the detector required the imposition of fiducial cuts, which were applied after the data and simulation were calibrated and assembled into spectra. Improvements in reconstruction allowed us to use a larger fiducial region for this analysis than was used for previous TWIST measurements [10,11,24]. The fiducial region adopted for this analysis required 0.50 < |cos θ| < 0.92, p < 51.5 MeV/c, 10.0 < p t < 39.7 MeV/c, and |p z | > 13.7 MeV/c. This region is shown in Fig. 4. The 51.5 MeV/c momentum cut removed the region of the spectrum used by the energy calibration; this also removed the only sharp feature of the spectrum, greatly reducing the sensitivity of the decay parameters to the resolution. The minimum |cos θ| cut eliminated high angle events, which underwent large amounts of multiple scattering and were difficult to reconstruct. The maximum |cos θ| cut was limited by increased reconstruction bias and a difference in reconstruction efficiency between data and simulation at small angles. The minimum p t cut eliminated tracks with radii too small to provide sufficient transverse separation in the hits for reliable fits. The maximum p t cut, in combination with the requirement that muon stopped within 2.5 cm of the detector axis, ensured that all accepted decay positrons were fully contained in the detector. The minimum |p z | cut removed tracks with wavelength comparable to a detector periodicity of 12.4 cm. The robustness of the analysis against variations in these cuts was checked using variations of 0.02 in cos θ or 0.5 MeV/c for the momentum cuts; neither the decay parameter measurements nor the systematic uncertainties were changed by more than 1 × 10 −5 .
The fiducial cuts were applied to the spectrum histogram, and were therefore influenced by the histogram binning. A histogram bin was included in the fit if the (p, cos θ) value of its centre passed the above fiducial cuts. The widths of the spectrum histogram bins were 0.5 MeV/c in momentum, and 0.02 in cos θ.

D. Spectrum Fitting
The decay parameters were determined by comparing the shape of the two-dimensional data spectrum to that of a simulated spectrum with matching conditions. The muon decay spectrum is linear in the parameters (ρ, η, P µ ξ, P µ ξδ), so the data spectrum can be described in terms of the MC spectrum as where S D = d 2 Γ D /dx d(cos θ) is the decay spectrum from data, and S M = d 2 Γ M /dx d(cos θ) is the spectrum simulated using Eqn. (4) to generate the muon decays. Note that the simulated spectrum is generated using hidden values of the decay parameters. ∂S/∂ρ etc. are the derivatives of Eqn. (4) with respect to the decay parameters. Radiative corrections were left out of these derivatives, under the assumption that the dependence of the radiative corrections on the decay parameters is negligible; the exception was ∂S/∂P µ ξ, which included the anisotropic radiative corrections to facilitate the consistent treatment of P µ ξ as a product. "Derivative spectra" were fully simulated and analyzed spectra in the same way as S M . To generate these spectra, the magnitudes of derivatives of Eqn. (4) with respect to the decay parameters (|∂S/∂ρ| etc.) were treated as probability distributions for the purposes of spectrum generation; a sign was associated with each event according to the sign of the derivative at that point, and this sign was used as a weight when the reconstructed event was included in the final derivative spectrum.
The coefficients ∆ρ, ∆(P µ ξ), and ∆(P µ ξδ) are the fit parameters, and represent the differences in the decay parameters between data and simulation. The muon decay spectrum in our fiducial does not accurately determine η because of the x 0 coefficient in Eqn. (5), so the value η = −0.0036 ± 0.0069 from the recent global analysis [25] was assumed.
Only differences between the real decay parameters in data and those assumed by the simulation were measured directly in this procedure. The values of the parameters assumed by the simulation (ρ h , δ h , ξ h ) were chosen randomly within 0.01 of the SM values while being constrained to physically meaningful values (e.g. ξδ/ρ ≤ 1), and these values were kept hidden throughout the analysis and the evaluation of systematic uncertainties. δ, and hence ∆δ, was extracted from ∆(P µ ξδ) as (ξ h δ h + ∆(P µ ξδ))/(ξ h + ∆(P µ ξ)). When measuring systematic uncertainties, it was sufficient to assume Standard Model values for the extraction of ∆δ. In this way the measurement outcome remained unknown until the values of the hidden parameters were revealed.
The ability of the fitted Monte Carlo spectra to describe the data can be viewed in terms two distributions: the angle-integrated momentum distribution of the decay positrons, which is governed by ρ (with η), and the momentum dependence of the decay asymmetry, governed by δ. The momentum spectrum of accepted events from a standard data set is compared with that from the corresponding fitted simulation in Fig. 5. The decay asymmetry as a function of momentum for accepted events from the same data set is compared with the fitted simulation in Fig. 6. Here asymmetry is defined as (B − F )/(B + F ), where B and F are the number of counts in the backward (cos θ < 0) and forward directions in a given momentum bin. The normalized residuals demonstrate the quality of the fit and the lack of bias. This study used the same data analyzed by Jamieson et al. for the direct measurement of the decay parameter P π µ ξ, published in 2005 [24]. However, the P π µ ξ analysis did not include the study of the ρ or δ systematic uncertainties and corrections. Furthermore, the present mea-surement used more advanced analysis and simulation techniques than were available for the previous study, as described below. A new, independent set of hidden decay parameters was generated for use in the simulation, and this new set was used exclusively for this analysis. Therefore, as with the previous measurements from TWIST, the measurement outcome could not be known until the hidden decay parameters were revealed.

V. SUBSTANTIAL CHANGES SINCE PREVIOUS STUDIES
There were a number of improvements to the experimental apparatus and data-taking techniques between 2002, when data were taken for the first TWIST measurements of ρ and δ [10,11], and 2004, when the data were taken for both the first TWIST measurement of P π µ ξ [24] and the present measurement of ρ and δ. The geometry and material of the muon stopping target used in 2004 was much better known than the target used previously. Moreover, additional monitoring and feedback controls were implemented for the taking of data in 2004. These improved the stability of the dipole magnets in the beamline, the flatness of the chamber foils, average muon stop position, and other aspects of the experiment, significantly increasing the quality of the data. The new muon beam monitoring chamber allowed significantly better characterization of the muon beam for input into the simulation, as well.
Since the previous measurements, the simulation and analysis were modified to account for the asymmetric construction of the drift chambers, where the wires were not centred between the foils (see Sect. II A). The track reconstruction used by the present analysis accounted for energy loss by the positrons, reducing reconstruction bias by several keV/c; other reconstruction improvements further reduced the reconstruction biases, particularly at low angles. Energy calibration used by the present analysis was performed as a relative calibration between simulation and data, resulting in a more accurate calibration.
As this analysis was nearing completion, a technique was developed to extract STRs directly from the data. This occurred too late to be incorporated into this entire analysis, but it was used to estimate the systematic uncertainty due to mismatch in the chamber response between data and MC (see Sect. VII B).

VI. SIMULATION VALIDATION
Since the data are fitted with simulated spectra to measure the muon decay parameters, the simulation must correctly reproduce physical and detector effects in an unbiased way; of particular importance are energy loss, multiple scattering, and the probabilities for the production of delta rays and bremsstrahlung. Direct comparisons between simulation and data were performed with specialized data not used for decay parameter measurements. In this mode, labeled "upstream stops", the momentum of the muon beam was reduced so that the muons stopped upstream of the DCs. Downstream-going decays then passed through the entire detector, and were reconstructed twice, first using only the upstream half of the detector, then using only the downstream half. Energy loss, scattering, helix fitter biases, and reconstruction resolution all resulted in differences in the properties of the two tracks. Distributions of ∆p = p DS − p US , ∆θ = θ DS − θ US , and other differences were used to examine reconstruction and physical effects, independent of the shape of the decay spectrum, and to compare the simulation directly to the data. About 8 × 10 7 upstream stop events were considered for this analysis, with 1.5×10 6 events accepted after cuts. About ten times that amount of simulated decays were generated.
As mentioned earlier, to first order energy loss is proportional to the amount of material a particle encounters; in TWIST's planar geometry this was proportional to 1/ cos θ. Multiplying ∆p by cos θ removes this firstorder angle dependence. Figure 7 shows the distributions of (∆p)(cos θ), for both data and simulation; the simulation has been normalized to data by integrated counts. The most probable value was −28.4 ± 0.1 keV/c for data and −29.65 ± 0.04 keV/c for MC; the FWHM width was 155.9 ± 0.1 keV/c for data and 141.64 ± 0.04 keV/c for MC. The slight difference in most probable energy loss was compensated by the energy calibration. The difference in the width of the distribution leads to a systematic error in the measurement of the decay parameters of 1.3×10 −4 , as discussed in Sect. VII E, and a correction was applied (Table III). Figure 8 shows the same distri- butions but for events where (∆p)(cos θ) < −1 MeV/c, demonstrating very good agreement between data and simulation over several orders of magnitude. This was used to determine part of the systematic uncertainty due to the simulation of positron interactions (Sect. VII D), which was found to be less than 0.7×10 −4 for both ρ and  Figure 9 shows the distributions of ∆θ for data and simulation; the simulation was normalized to the data as before. The most probable value was −0.97 ± 0.02 mrad for data and −0.581 ± 0.007 mrad for MC; the FWHM width was 29.75 ± 0.02 mrad for data and 29.159 ± 0.007 mrad for MC. The difference in the most probable ∆θ values between data and MC is very small compared to the approximately 10 mrad angular resolution and so did not affect the decay spectrum. Unlike the energy loss, the slight difference in the width of the ∆θ distributions had a negligible effect on the measurement of the decay parameters. Figure 10 shows the same distributions but including events with large ∆θ, again demonstrating very good agreement between data and simulation over many orders of magnitude. The upstream stops technique was also used to compare reconstruction efficiencies in data and simulation, as a function of momentum and angle. A simplified definition of reconstruction efficiency was used to avoid introducing artifacts at the edge of the fiducial region and other complications. Let r d (p, cos θ) be the number of downstream tracks in a given (p, cos θ) bin for which an upstream track was also reconstructed, and let T d (p, cos θ) represent the total number of downstream tracks in that bin. Then the upstream efficiency was defined as ǫ u (p, cos θ) = r d /T d . Downstream efficiency was similarly defined. Figures 11 and 12 show the (Data-MC) differences in reconstruction efficiency, as functions of momentum and cos θ respectively; the differences of upstream efficiencies are shown, and the downstream efficiencies were similar.
Each bin represents an average across the fidu- cial region. The region around the upstream stops beam momentum of ∼ 25 MeV/c has been excluded to avoid an artifact in the efficiency measurement caused by the beam positron phase space. The simulation and data show good agreement in reconstruction efficiency over the entire fiducial region, and the slight differences are independent of momentum and angle. Overall, we find that all but (5.7 ± 0.2) × 10 −4 [(5.12 ± 0.05) × 10 −4 ] of the upstream stops events in the data [simulation] that contain a reconstructed downstream track within the fiducial region also contain a reconstructed track in the upstream half of the detector; similarly, all but (8.2 ± 0.2) × 10 −4 [(7.87 ± 0.06) × 10 −4 ] of the events that contain a recon- structed upstream track within the fiducial region also contain a reconstructed track in the downstream half of the detector. The differences in upstream and downstream reconstruction efficiency give rise to a systematic error of less than 0.6×10 −4 , as discussed in Sect. VII B.

A. General Procedure
A systematic uncertainty is an uncertainty in the decay parameters as the result of an unknown or unaccountedfor value or variation of some experimental parameter. A summary of the systematic uncertainties for this measurement is given in Table II by category; these are described in more detail below. A systematic error repre- sents a known experimental bias rather than an uncertainty; where these errors were found to be significant, a correction was applied. Very small errors were simply included in the total systematic uncertainties. Table III lists the corrections applied to the present measurement. Several systematic uncertainties important to the measurement of P π µ ξ were not studied for this measurement. The leading systematic uncertainty in the previous measurement of P π µ ξ by TWIST [24] arose from uncertainty in our knowledge of the muon beam emittance, and its impact on the muon depolarization as the beam entered the magnetic field. The depolarization of the muons due to the muon beam emittance was found to be 5-7×10 −3 depending on the beam used, and the associated systematic uncertainty was ∆P π µ ξ = 0.0034. The muons depolarized further while at rest in the stopping target at a rate of (1.6 ± 0.3) × 10 −3 µs −1 , leading to an associated systematic uncertainty of ∆P π µ ξ = 0.0012. We have no additional data to improve our knowledge of the muon beam used at the time. Furthermore, whereas this issue had a very large impact on the determination of P π µ ξ, its impact on the measurement of ρ and δ is negligible.
Systematic uncertainties and errors were studied using the spectrum fitting technique described by Eqn. (10). A new simulated decay spectrum was produced by exaggerating an experimental parameter at the simulation or analysis stage, and this was fit against a simulation generated using the standard parameters. The fit results ∆ρ and ∆δ, scaled by the amount of exaggeration, provided a direct measurement of the systematic uncertainty associated with the experimental parameter.
All systematic errors and uncertainties were determined prior to revealing the hidden decay parameters assumed by the simulation.

B. Chamber Response
To determine the effect of inaccuracies in the garfield-generated STRs, analyses using STRs derived directly from the data and the simulation were fit against the corresponding standard analyses. Data-derived and simulation-derived STRs for TWIST have only recently been developed. They were derived for data and simulation from the means of the time residuals (the difference between reconstructed hit time and measured hit time) as a function of hit position within the drift cell, during reconstruction using garfield-generated STRs. These residuals were used to modify the garfield-generated STRs, and the process was iterated until it converged. The data-derived STRs automatically account for variations in temperature, foil positions, etc., within the de-tector, as well as compensating for some residual biases in the helix reconstruction. They lead to reduced χ 2 values from the helix fitter. They also provide improved momentum resolution and a better match between data and simulation compared to that described in Sect. VI. Thus, we conclude that they are more correct for use in the analysis than the STRs generated by garfield.
The tests determined by how much the measured decay parameters were affected by the use of garfieldgenerated STRs in the standard analysis. The difference between the effect in simulation and the effect in data demonstrated the amount by which the measured decay parameters in a standard data-MC fit would shift; a correction was applied to the final measurement to account for this (Table III). The uncertainty on this difference, at 3×10 −4 for ρ and 5×10 −4 for δ, represented the systematic uncertainty due to inaccuracies in the STRs. This dominated the chamber response uncertainties.
Other aspects of the chamber response studied included the uncertainties in the DC foil positions, the asymmetries in the decay reconstruction efficiency, and electronics-related wire time offsets. None of these represented uncertainties greater than 0.6 × 10 −4 for either ρ or δ, and most were much smaller.

C. Energy Scale
To translate uncertainties in the energy calibration parameters into uncertainties in the decay parameters, new spectra were created by applying energy calibration, with each calibration parameter exaggerated, to a standard simulation. These exaggerated spectra were fit against a spectrum created using the standard calibration to measure the effects on the decay parameters, which were found to be 3×10 −4 for ρ and 4×10 −4 for δ. The systematic uncertainty due to the energy calibration was the dominant uncertainty related to energy scale.
Other aspects of the energy scale studied included the assumed behavior of the energy calibration with momentum, and errors in the shape of the magnetic field map used for analysis. The energy calibration was assumed to be independent of the positron momentum; using a calibration proportional to the momentum changes the decay parameters by less than 0.8 × 10 −4 for both ρ and δ. The simulated magnetic field used for analysis and simulation was compared against the measured magnetic field map; the differences were found to influence the decay parameters at the level of 0.7×10 −4 for both ρ and δ.

D. Positron Interactions
The production of delta rays and bremsstrahlung were the most important discrete positron interactions for TWIST. The simulation of these was validated using the upstream muon stops described in Sect. VI above. The rate of delta ray production in simulation was compared to that in data using the ratio R δ = N 12 /N 11 , where N 12 is the number of events with one upstream track and two downstream tracks, and N 11 the number of events with one upstream track and one downstream track. We find R δ = (1.432 ± 0.003) × 10 −2 in the upstream stops data. R δ was also measured in several additional upstream stops simulations with the delta ray production cross-section multiplied by various factors, as shown in Fig. 13. A linear fit to R δ as a function of the cross- section multiplier was used to translate the R δ value from data into a difference in delta ray cross-sections between data and simulation, conservatively assuming that all events of this topology are due to delta rays. This is not true in practice, since R δ does not go to zero in the simulation when the delta ray production cross-section in GEANT is set to zero (see Fig. 13); however, the conservative assumption was sufficient for estimating the systematic uncertainty. The effective delta ray production cross-section determined from this method was found to be approximately 18% lower in the simulation than was seen in data.
The rate of bremsstrahlung production (R B ) in simulation was compared to that in data by counting the number of through-going positrons whose change in momentum was p DS − p US < ( ∆p − 1 MeV/c), normalized to the total number of reconstructed events; here ∆p represents the most probable value of ∆p. We find R B = (1.42±0.01)×10 −2 . As with the study of delta ray production, R B was also measured in several additional upstream stops simulations with the bremsstrahlung production cross-section increased by various factors, as shown in Fig. 14. As before, a linear fit was used to translate the R B value from data into a difference in bremsstrahlung cross-sections between data and simulation; the simulation was found to agree with the data to within 2%.
The simulations produced with the delta ray or bremsstrahlung cross-sections increased by a factor of three were fit against standard simulation. The systematic uncertainty due to errors in the simulation of delta ray production was found to be 1.5×10 −4 for ρ and 0.9×10 −4 for δ; that due to errors in the simulation of bremsstrahlung production was less than 0.7×10 −4 for either parameter.
Materials outside the sensitive region of the detector could scatter particles back inside, resulting in extraneous hits that the track reconstruction must sort out. The number of backscattered positrons, normalized by the number of muons, was studied in the data and the simulation. A new simulation was produced with an additional plate of aluminum downstream of the detector as a backscattering source, and this was fit against the standard simulation. The effect on the decay parameters was less than 0.4×10 −4 .

E. Resolution
Reconstruction resolution in momentum and angle smeared the decay spectrum. The RMS angle resolution was about 5-15 mrad, and the simulation agreed with the data to within 3 mrad. The RMS momentum resolution was 0.040-0.120 MeV/c, and was smaller in the simulation by 0.007 MeV/c. Both types of reconstruction resolution varied with energy and angle of the decay positron. Since the spectrum was inherently smooth, however, these resolutions did not significantly distort the decay parameters directly. The largest effect was at the endpoint, where the momentum resolution changes the shape of the edge; the 7 keV/c difference in resolution between simulation and data had a small effect on the energy calibration, which resulted in a distortion in the decay parameters. Resolution was measured in the data and MC as a function of momentum and angle using decays from upstream muon stops, by measuring the widths of the energy loss and scattering distributions. New spectra were produced by smearing the reconstructed momentum or angle of each event in the simulation before including it in the spectrum and calculating the energy calibration. The effect of the angle resolution on the decay parameters was found to be less than 0.4×10 −4 for both ρ and δ. The effect of the momentum resolution on the decay parameters was 1.2×10 −4 for ρ and 1.3×10 −4 for δ, and a correction was applied (Table III).

F. Other Sources of Error
A number of other possible sources of systematic error were studied and found to be small.
The effect of the uncertainty in the chamber alignment on the decay parameters was tested by analyzing a simulation using distorted chamber alignments, and fitting against the standard simulation. These were found to affect the decay parameters by less than 0.03×10 −4 . Errors in the transverse and longitudinal length scales of the detector translated directly into errors on the transverse and longitudinal components of the momentum; these were tested by distorting these momentum components for reconstructed decays before including them in a new spectrum. The effect was less than 0.3×10 −4 .
New simulations were produced with exaggerated muon or positron beam rates, to test the sensitivity of having an incorrect pile-up rate in the simulation. This showed the systematic uncertainty due to the simulation of beam intensity to be less than 0.2×10 −4 for both beam components.
As described above, the highest order of radiative corrections used for this analysis was O(α 2 L). The O(α 2 ) radiative corrections [22] represented the theoretical error in the O(α 2 L) corrections, and this determined the theoretical uncertainty for this analysis. A new simulation was produced with the O(α 2 L) radiative corrections exaggerated, and this was fit to the standard simulation. The results show that the theoretical uncertainty in the measured decay parameters is less than 0.3×10 −4 .
During decay parameter fits, η was normally fixed to the world average value of -0.0036 [25]. The standard fit between data and MC was repeated with η raised or lowered by one standard deviation (δη = ±0.0069), giving the correlations ∂ρ/∂η = 0.0162, ∂δ/∂η = 0.0015, and ∂P µ ξ/∂η = 0.0155. This led to an uncertainty in ρ due to the uncertainty in η of 0.00011. Future improvements in the knowledge of η can be used to reduce this systematic uncertainty directly.
Since this measurement was finalized, development of the simulation and analysis software has continued for the analysis of additional data. During this process, an error in the particle identification was found, which caused good events to be misclassified and discarded, affecting about 0.6% of events in the fiducial region. This distortion was primarily linear in |cos θ|, which is a shape not reflected in any of the derivative spectra; furthermore, the distortion was the same in the analysis of both data and simulation, leaving the derived decay parameters unaffected. This software error has been fixed for future analyses, but has been listed here for completeness.

VIII. RESULTS
The values of ρ and δ measured from fitting each data set with its corresponding simulation are listed in Table IV. After taking weighted averages and applying the corrections in Table III, we find where the last uncertainty is due to the uncertainty in η, and δ = 0.75067 ± 0.00030(stat) ± 0.00067(syst).
Both results are consistent with the Standard Model values of 3/4. These represent a factor of two improvement over the previous TWIST measurements [10,11]. The typical correlation coefficient between ρ and δ from the decay parameter fits is +0.15 for the fiducial range adopted here. Correlations also exist between the systematic uncertainties in ρ and δ. The contributions to σ 2 ρδ associated with the chamber response and positron interaction systematics are also positive and somewhat larger than the statistical contribution from the decay parameter fits. In contrast, the energy calibration systematics for ρ and δ are strongly anti-correlated, and provide the dominant contribution to σ 2 ρδ in the present measurement. Overall, we find that the correlation coefficient between the total uncertainties in ρ and δ is −0.16.
The decay parameter fits also determined a value for P µ ξ, which can be converted to a value for P π µ ξ by correcting for known depolarization effects. After applying corrections corresponding to the ρ and δ corrections in Table III, the fits give P π µ ξ = 1.0025 ± 0.0004(stat). As explained in Sect. VII A, the systematic uncertainties important to the measure of P π µ ξ were not revisited during this analysis, and this value should not be considered a revised measurement. When differences in the two analyses are considered, the value of P π µ ξ found here is consistent with the published TWIST measurement of P π µ ξ = 1.0003 ± 0.0006(stat) ± 0.0038(syst) [24].

IX. GLOBAL ANALYSIS OF MUON DECAY
In 2005 Gagliardi et al. performed a global analysis of all available muon decay data [25], including earlier TWIST measurements of ρ and δ [10,11]. That analysis has been repeated, incorporating the TWIST measurement of P π µ ξ [24] and the new ρ and δ measurements presented here; the correlation factor of −0.16 between the TWIST ρ and δ measurements has also been included in the calculation. All other input values are the same as in the analysis of Gagliardi et al.
In brief, the global analysis used a Monte Carlo method similar to that of Burkard et al. [26] to map out the joint probability distributions for 9 independent variables, Q RR , Q LR , Q RL , B LR , B RL , α/A, β/A, α ′ /A, and β ′ /A. Each of these parameters is a bilinear combination of the weak coupling constants g γ ǫµ . The decay parameters could then be written in terms of these independent variables. Table V shows the results of this global analysis, as well as the results of the analysis of Gagliardi et al.
The 90% confidence limits are given for the independent variables listed above, and global best-fit values of the decay parameters ρ, δ, and η are given. The present analysis represents significant improvements in the limits on Q LR and B LR , and tightens several of the other limits. It is interesting to note that the global analysis significantly reduces the uncertainty in the value of ρ, from a total of 0.00063 to 0.00035. The values of the Q ǫµ from this global analysis can be used in Eqn. (2) to place limits on the magnitudes of the weak coupling constants g γ ǫµ ; the exceptions are g V LL and g S LL , which are determined more sensitively from inverse muon decay, e − ν µ → µ − ν e . The limits determined with this method are listed in Table VI, along with the values from the previous global analysis [25]. The present analysis represents a reduction of approximately 16% in the limits for |g S LR | and |g T LR |, and a 30% reduction for |g V LR |.

X. CONCLUSION
This new measurement of the muon decay spectrum is a factor of two more precise than previous measurements [10,11]. It is consistent with Standard Model predictions, placing more stringent limits on "new physics" in the weak interaction. New indirect limits on the value of P π µ ξ can be obtained using these values of ρ and δ, in combination with the the measurement of P π µ ξδ/ρ > 0.99682 at 90% confidence by Jodidio [27,28] and using the constraint ξδ/ρ ≤ 1 (required to obtain positive definite decay probabilities from Eqn. (4)). Accounting for the correlation coefficient of −0.16 between the ρ and δ uncertainties, we find 0.99524 < P π µ ξ ≤ ξ < 1.00091 at 90% confidence. This is a significant improvement over the previous indirect limit of 0.9960 < P π µ ξ ≤ ξ < 1.0040 [11].
The quantity Q µ R = Q RR + Q LR represents the total probability for a right-handed muon to decay into any type of electron, a process forbidden under the Standard Model weak interaction. The new measurements of ρ and δ lead to the new limits on Q RR and Q LR shown in Table V, and hence to a new 90% confidence limit upper bound on the combined probability Q µ R < 0.0024, a slight improvement over the limit of Q µ R < 0.003 from the previous global analysis [25].
Under left-right symmetric models, ρ > 0.75 is forbidden, so the general measurement of ρ can be converted into a 90% confidence limit lower bound within LRS mod-els: ρ > 0.7493 (compared with ρ > 0.7487 from the ρ measurement previously published by TWIST). The relation ρ ≃ 3 4 (1 − 2ζ 2 g ) then gives a 90% confidence limit of |ζ g | < 0.022, a significant improvement over the limit of |ζ g | < 0.030 for the previously published TWIST value of ρ.
The final phase of TWIST is in progress. Additional data are in hand, taken in 2006 and 2007, and new analysis is underway. Further improvements to the simulation and analysis are being implemented, including the use of measured drift time tables and improved track reconstruction algorithms, and the resulting better agreement in the momentum resolution between simulation and data. These improvements are expected to lead to additional factor-of-two reductions in the uncertainties on ρ and δ, providing another incremental improvement to searches for new physics.