A Revised Experimental Upper Limit on the Electric Dipole Moment of the Neutron

We present for the first time a detailed and comprehensive analysis of the experimental results that set the current world sensitivity limit on the magnitude of the electric dipole moment (EDM) of the neutron. We have extended and enhanced our earlier analysis to include recent developments in the understanding of the effects of gravity in depolarizing ultracold neutrons (UCN); an improved calculation of the spectrum of the neutrons; and conservative estimates of other possible systematic errors, which are also shown to be consistent with more recent measurements undertaken with the apparatus. We obtain a net result of $d_\mathrm{n} = -0.21 \pm 1.82 \times10^{-26}$ $e$cm, which may be interpreted as a slightly revised upper limit on the magnitude of the EDM of $3.0 \times10^{-26}$ $e$cm (90% CL) or $ 3.6 \times10^{-26}$ $e$cm (95% CL). This paper is dedicated by the remaining authors to the memory of Prof. J. Michael Pendlebury.

We present for the first time a detailed and comprehensive analysis of the experimental results that set the current world sensitivity limit on the magnitude of the electric dipole moment (EDM) of the neutron. We have extended and enhanced our earlier analysis to include recent developments in the understanding of the effects of gravity in depolarizing ultracold neutrons; an improved calculation of the spectrum of the neutrons; and conservative estimates of other possible systematic errors, which are also shown to be consistent with more recent measurements undertaken with the apparatus. We obtain a net result of d n ¼ −0.21 AE 1.82 × 10 −26 e cm, which may be interpreted as a slightly revised upper limit on the magnitude of the EDM of 3.0 × 10 −26 e cm (90% C.L.) or 3.6 × 10 −26 e cm (95% C.L.).

I. INTRODUCTION
Measurement of the electric dipole moment (EDM) of the neutron (or of other fundamental particles) provides an extremely sensitive approach to investigating potential new sources of CP violation and physics beyond the Standard Model.
The experimental technique underlying the most sensitive measurement to date [1,2] of the neutron EDM has been discussed extensively in an earlier publication [3]. The data were collected at ILL, Grenoble, between 1998 and 2002. In this article, we focus solely on the analysis, which we have carried out anew. We begin, in Sec. II, with a description of the determination of the neutron resonant frequency via the Ramsey method of separated oscillatory fields in conjunction with the cycle-by-cycle corrections from the mercury comagnetometer.
A number of cuts were applied to the data. These are discussed in Sec. III.
For each data-taking run, of typically 1-2 days' duration, a value of the EDM was determined from the slope of a linear fit of the (field-drift corrected) neutron frequency versus applied electric field E. However, in the presence of a magnetic-field gradient, the mercury (and, to a far lesser extent, neutron) frequency acquires a component that is linear in the applied electric field [4], and thus mimics the signal for an EDM. We have no direct measurement of the applied B-field gradient, but since the neutrons have a lower center of mass than the (thermal) mercury atoms, any change in the gradient results in a change in the neutron-tomercury frequency ratio. A linear dependence of the measured apparent EDM as a function of this ratio therefore emerges, as shown in Fig. 2 of [1] and reproduced here in Fig. 1, with a positive (negative) slope for a downwards (upwards) direction of the applied magnetic fieldB. We discuss this further in Sec. IV below.
The offsets and profiles of these lines can be affected by a number of factors, including nonuniform and/or horizontal magnetic-field gradients, the gravitationally enhanced spectrum-dependent depolarization of the ultracold neutrons (UCNs), and even the rotation of the Earth [5]. These are discussed in Sec. V.
In order to measure and compensate for such effects, measurements were made in an auxiliary bottle that used a lid of adjustable height, so that the neutron and mercury frequencies could be determined as a function of height for a given B-field configuration. These measurements are described in Sec. VI.
In Sec. VII, we describe how all of these measurements are brought together in a global fit, in order to determine the necessary corrections to the initial EDM estimate obtained from the fit to the two crossing lines.
In Sec. VIII we discuss further possible contributions to systematic errors.
We close, in Sec. IX, with a final summary and conclusion.
For convenience in what follows, and following the convention of [1], we define where ν, γ are the respective frequencies and gyromagnetic ratios of the two species (neutrons and Hg), and we also introduce In this context, R 0 is usually specified in parts per million (ppm). Throughout this analysis we have used values of γ Hg ¼ 7.590118ð13Þ Hz=μT and γ n ¼ −29.1646943ð69Þ Hz=μT, based on the measurements of Cagnac [6] and of Greene et al. [7]. Both of these measurements were made relative to the shielded proton gyromagnetic ratio in pure water, for which the most recently measured value is γ 0 p ¼ 42.5763866ð10Þ Hz=μT [8]. From these results, we observe that R 0 ¼ 0 when ν n =ν Hg ¼ 3.8424560ð66Þ, which represents an accuracy of 1.73 ppm. We incorporate the recent measurement of γ n =γ Hg ¼ −3.8424574ð30Þ by Afach et al. [9] into the later stages of our analysis.

II. NEUTRON FREQUENCY DETERMINATION
The neutron frequency fitting procedure is described in [10] and [3]. For any measurement cycle (equivalently referred to as a "batch" in [3]), the first-order estimate ν n;Hg of the neutron resonant frequency is determined from the measured mercury precession frequency ν Hg by ν n;Hg ¼ γ n γ Hg ν Hg : The Ramsey line shape may then be written as where δν ¼ ν n;Hg − ν 1 is the difference between the firstorder resonant-frequency estimate ν n;Hg and the frequency ν 1 of the applied rf pulses. The linewidth Δν is the width at half height of the central fringe, and is given by where T is the free precession time and t is the duration of each of the π=2 spin-flip pulses. The polarization α is related to the fringe maximum and minimum N max ; N min as . From [1]. Measured EDM (binned data) as a function of the relative frequency shift of neutrons and Hg. The solid red line is a linear best fit.
The phase incorporates the difference between the true resonant frequency ν 0 and the first-order estimate ν n;Hg derived from the mercury. This difference arises for several reasons. First, there is the inherent uncertainty in the ratio γ n =γ Hg as discussed above. Second, and as also discussed above, ϕ alters because the neutrons were of such low energy that they populated preferentially the lower portion of the storage volume, whereas the mercury sampled the entire volume uniformly; there was a difference in height of the centers of gravity of the two systems of a few millimeters.
In the presence of a vertical gradient ∂B z =∂z in the magnetic field, the average field sampled by the mercury would have been slightly different from that sampled by the neutrons. The normalized frequency ratio is therefore related to this height difference Δh by where the þ sign corresponds to B 0 downwards. Third, as will be discussed below, gradients in horizontal components of the magnetic field can increase the neutron precession frequency while leaving the mercury unaffected. The rotation of the Earth also has a part to play, as do minor effects such as light-induced frequency shifts in the mercury.
These factors, which remained essentially constant throughout each run, caused the departure of R a [Eq. (1)] from unity. There are, in addition, two factors that could provide extremely small cycle-to-cycle variations in ϕ. These are the small mismatch in temporal overlap of the mercury and neutron frequency measurements, as discussed in Sec. III E, and a genuine neutron EDM or some effect that mimics it.
The high voltage (HV) was applied in a simple pattern of alternating sign, with an additional four cycles at E ¼ 0 between each 16 cycles of HVat a given polarity. The initial polarity for each run was chosen randomly. In principle, one could treat each E value within a run as a distinct data set, with its own independent Ramsey-curve fit. Our approach instead has been to fit a single Ramsey curve to each entire run, with the frequency shifts (which are very small perturbations, well within the noise) then being analyzed on a cycle-by-cycle basis. This has the benefit of allowing analysis of runs containing insufficient data within each E group for a reliable Ramsey fit.
While the polarity of the E field could be switched automatically under computer control, changing the direction of the magnetic field was a much more onerous procedure. This was normally undertaken every few weeks.
For each run, typically consisting of several hundred measurement cycles, the neutron counts were fitted to Eq. (4), as shown in Fig. 2, to provide values ofN s , α s and an average value of the phase ϕ s for each of the two spin states s. (For maximum sensitivity, measurements were taken repeatedly at four working points, close to the halfheight of the central Ramsey valley.) For the first iteration of each of these two fits, the uncertainty σ y allocated to each data point was simply the square root of the number of neutrons counted in that particular spin state. The uncertainty in the fitted mercury frequency, which resulted in an uncertainty σ x in the x coordinate of each data point, was incorporated by calculating an "indirect" error σ I from the slope of each first-iteration fitted curve, and adding it in quadrature to the original error bar to obtain a new uncertainty for use in a second iteration of the fit. This correction was, however, generally negligible, as the mercury usually measured the magnetic field with much higher precision than did the neutrons. An overall average (weighted as the fit uncertainties)φ of the two phases ϕ s was calculated, giving a common phase offset for the entire run and for both spin states. Equation (4) was then inverted to yield an individual phase shift δϕ s;i for each data point i and for each spin state: A phase shift δϕ i was calculated for each data point by averaging the δϕ s;i over the two spin states. This averaging helped to remove nonstatistical fluctuations in the total neutron flux, which could, exceptionally, vary by a percent or two during a run. Any shift δν 0i in the neutron resonant frequency from the value ν 0 predicted by the fitted Ramsey curve should be proportional to this phase shift δϕ i : In the presence of an EDM, this frequency shift should be directly proportional to the strength of the electric field. A straightforward linear least-squares fit of the frequency shifts δν 0i as a function of the applied electric fieldẼ yielded a value for the measured EDM d meas , with its associated uncertainty, for each run.
We have recently carried out detailed modeling of our storage-chamber geometry using the Opera [11] finiteelement analysis package, and we find that the average electric field within the volume is 1.1% lower than the nominal value determined from the applied voltage divided by the separation between the electrodes.
For completeness, we note that the statistical uncertainty to be expected from an EDM measurement based upon the Ramsey technique is [3] Care was taken to verify that the analysis delivered the correct sign of EDM: in particular, withẼ andB fields parallel, an increase in neutron precession frequency characterizes a negative EDM.

III. APPLIED CUTS
In this section we list all cuts applied to the data. No systematic dependence of the EDM signal upon any of these cuts was observed. After all of these following cuts were applied, 545 runs containing 175,217 measurement cycles and 2.5 × 10 9 neutrons remained.

A. Manual cuts
After inspection, 10 runs were rejected for various reasons relating to the failure of hardware components such as amplifiers, valves, the HV supply and so on. Some individual cycles within each run were also cut for similar reasons: the majority of these were due to problems with the delivery of neutrons. A further 26 runs were rejected because they had 19 or fewer measurement cycles, which is not sufficient for the HV polarity reversal required for an EDM measurement.
These manual cuts removed about 4% of the available data, after which the above-mentioned 545 runs remained.

B. Partial runs
During part of the data taking, an adjacent experiment was using a superconducting magnet. When turned on or off, R a changed by typically 2-3 ppm. The 11 measurement runs during which this occurred were therefore split into partial runs for separate analysis.
The χ 2 =ν distribution from the cycle-by-cycle online Hg frequency fitting is shown in Fig. 3, which is reproduced from Fig. 10 of [3]. The red "expected" distribution curve is based upon the model of a perfectly constant frequency, whereas in fact we do expect some slow variation over time within a cycle. That being the case, a large χ 2 =ν is not necessarily indicative of a problem with the frequency measurement, and we therefore retain a generous portion of this distribution beyond the vicinity of the peak. A cut was made at χ 2 =ν ¼ 3, by which time the tail is fairly flat.
This cut removed nearly 7% of the remaining data. We note in passing the discontinuity at χ 2 =ν ¼ 4. This arises because, beyond this point, the online fitting procedure attempted to correct for potential hardware errors (e.g. gain saturation or a missed reading of the analog-todigital converter), as discussed in [3].

D. Mercury frequency uncertainty
The fitted Hg frequency sometimes had a large uncertainty, particularly if the depolarization time was short. The distribution of these uncertainties is shown in Fig. 13 of [3], and again here on a semilog plot in Fig. 4. A typical value is 1-2 μHz, which (when scaled with γ n =γ Hg ) corresponds to approximately a factor of 5 better than the typical inherent neutron frequency uncertainty from counting statistics. A cut was made at 25 μHz, at which point the relationship is approximately inverted, with the mercury uncertainty entirely dominating the frequency-ratio measurement.

Data Expected
Beyond this, the data are so imprecise that they make essentially no contribution at all. This cut removed 4% of the remaining data.

E. Magnetic-field jumps
The distribution of Hg frequency jumps, i.e. the difference in Hg frequency between a given cycle and the previous cycle, is shown in Fig. 14 of [3]. There are broad tails due to occasional sudden changes in field, for example due to the movement of an overhead crane or to a mechanical disturbance to the mu-metal shields.
The mercury and the neutron frequency measurements do not have perfect temporal overlap. As discussed in [3], the Hg frequency was determined by fitting a 15 s averaging period at either end of the cycle in order to determine the phase, and hence the integrated phase difference accumulated over the time between them. The corresponding phases for the neutrons, on the other hand, are determined by the 2 s spin-flip pulses at either end. The field jumps can occur at any point during the 220 s measurement cycle. If they occur outside the Ramsey sequence, they are of little concern. If they occur between the 15 s windows, they can be regarded as appropriately compensated. That leaves a potential 30 s period during which there is a risk of incomplete compensation, i.e. 1=7 of such jumps can be expected to affect the measurement to a greater or lesser extent. On average, these potentially risky field jumps would occur halfway through the 15 s window, i.e. about 1=20 of the way through the Ramsey measurement. The mercury frequency-jump distribution was truncated at AE60 μHz, corresponding to a change in neutron frequency of 230 μHz over this 1=20 of the Ramsey period. When averaged over the entire Ramsey period, the neutrons would therefore see a frequency shift of up to 11 μHz, or 0.4 ppm, that would not be compensated properly by the mercury. Including the aforementioned 1=7 probability, this corresponds (even at the extremes of the frequency-jump distribution) to a potential error in R 0 of 0.06 ppm on the rare cycles within which such jumps occur, to be compared with a typical statistical uncertainty on the neutron frequency of about 0.7 ppm.
The frequency-jump cut removed 3% of the remaining data.

F. First-cycle cut
The first cycle of any run is different from any of the others, as the neutron trap and guides are initially empty; for other cycles there is likely to be some remnant population from the previous cycle. In consequence the first cycle often has an anomalously low total neutron count.
This cut removed 0.2% of the remaining data.

G. Ramsey residuals outlier cut
The data points of Fig. 5, reproduced from Fig. 5 of [3], show the distribution of stretch values r i of the fits to the Ramsey curve: where ν i is the calculated frequency of the ith cycle, σ i is its uncertainty and ν R i is the expected frequency for that cycle as defined by the Hg frequency, the applied rf and the Ramsey curve function. Ideally, and in the absence of any EDM-like signals, this distribution would be expected to be a Gaussian of unit width. The continuous line is a Gaussian of width 1.06. The true distribution departs from this Gaussian at about 4σ. The few points lying outside this range tend to be associated with runs that have other known problems, for example with intermittent failure of the neutron delivery system. A cut was therefore made at AE4σ. Because of the symmetric way in which the data are taken, cutting the tails from this distribution cannot of itself induce a false EDM signal.  This cut removed a further 0.05% of the remaining data.

H. High voltage
A single measurement cycle with a leakage current in excess of 60 μA was removed from the data set. All remaining measurement cycles had leakage currents below 10 μA, with the great majority being of the order of a few nanoamperes; the distribution for both polarities is shown in Fig. 6, which reproduces Fig. 20 of [3]. No further cuts relating to the HV were made: the mercury magnetometer is relied upon to compensate for any residual effects from this source. Under the assumption that leakage currents would run along preferred established paths-rather than averaging out over different paths-then if such leakage currents were to generate EDM signals, measurement cycles with high leakage currents would be expected to show a greaterthan-average departure from the Ramsey curve, and thus to have high Ramsey-residual stretch values [Eq. (14)]. In order to quantify this, the data were binned by leakage current, and the distributions of stretch values within each bin were fitted to Gaussians. The widths of these Gaussians are plotted against leakage current in Fig. 7. No consistent trend is visible. Furthermore, after processing, no dependence of the measured EDM on leakage current was observed, as discussed in Sec. VIII E.

I. Frequency ratio
Measurements were undertaken with a range of different applied B-field gradients ∂B z =∂z by preadjusting currents in field-trimming coils. (Initially, this arose through trial and error as the system was optimized; we then settled upon a more consistent configuration for each B 0 direction, with a few runs later on used to explore the effects of deliberately large gradients.) As discussed in Sec. IV, these field gradients induced false EDM signals. Any nonlinearity in this effect would appear as a systematic departure of the points in Fig. 1 from their fitted lines. We define the stretch values where d meas i AE σ i is the measured EDM for run i, and d f i is the corresponding value from the fitted line. The distribution of these stretch values is shown in Fig. 8. There being no evidence of any correlation, and none expected, it was decided not to impose any restriction upon the range of R a .

A. Introduction
A false EDM signal d n;Hg;f arises when the trapped particles experience a gradient ∂B z =∂z in the presence of the electric field E [4]. The effect is easiest to understand in the case of an azimuthally symmetric field with a vertical   gradient ∂B z =∂z, i.e., a slightly trumpet-shaped field, as shown in Fig. 9.
Because∇ ·B ¼ 0, field lines either enter or emerge from the side walls of the cell, giving a radial field of strength that is proportional to the radius r: Consider now a particle moving at speed v that crosses the storage cell close to its diameter, as shown in Fig. 10. As it travels through the electric field, the particle experiences in its own rest frame an additional magnetic field [12] B above and beyond the laboratory magnetic fieldB 0 . At the start of its trajectory, just left of center at the bottom, it is subject to the radial field B r as well as to the sidewaysB v component, yielding a diagonal resultant. As it traverses the trap, the B r component shrinks and then reverses direction, causing a smooth rotation of the net additional effective B field. Eventually, the particle reaches the far end of the trap, is reflected from the wall and begins its trajectory back. However, theB v component then faces in the opposite direction: therefore, after a discrete jump at the point of reflection, the additional net effective field component continues to rotate in the same direction. The particle thus sees a rotating field in the x-y plane, which, through the Ramsey-Bloch-Siegert [13,14] mechanism, "pulls" its resonant frequency away from the central value. When averaged over both directions of a given trajectory, the net frequency shift is proportional to E, and it therefore mimics an EDM signal. The size of this false-EDM effect depends upon the relative magnitudes of the orbital-trajectory frequency of the particle and its Larmor precession frequency. The neutrons are in the adiabatic limit, where their rapidly precessing spins can follow the variations in field. In that case, the false EDM is given by [4] δν where v xy is the average transverse particle velocity. The origin of this shift is a geometric (or Berry's) phase [15], and therefore-as pointed out explicitly in [16]-is independent of the coupling to the magnetic field; the gyromagnetic ratio is absent from this equation. The mercury lies in the nonadiabatic regime, and its frequency shift in such a field gradient is [4] δν ¼ where r B is the radius of the storage bottle. This has been independently measured [16]. A more general expression, valid for arbitrary magnetic fields in the nonadiabatic limit, was derived more recently [17]: where the brackets refer to the volume average over the (arbitrarily shaped) trap. FIG. 9. From [4]. Showing the shape of the B 0 field lines, when there is a positive gradient ∂B z =∂z, shown in relation to an outline of the trap used to store 199 Hg atoms and ultracold neutrons. Additional fields having lines that both enter and leave through the sidewalls, like the one on the right-hand side, do not affect the false EDM signals that are generated.
FIG. 10. From [4]. Showing the B xy fields (in black) seen by a particle going back and forth close to the y-axis. Going towards positive y, the B xy field rotates steadily anticlockwise by about 70°as drawn. The first reflection of the particle towards negative y causes an instantaneous anticlockwise rotation by about 110°as drawn. The same two rotations occur on the path to, and at, the second reflection. The size of the rotations depends on the size of B 0r =B v .
Since the mercury is used to compensate for shifts in the magnetic field, any EDM-like component contributing to the mercury frequency will affect the measurement of the neutron EDM. We note that [17] contains a sign error in its Eq. (20) and Fig. 2, arising from this transferral of the mercury false EDM to the neutron signal.
In this experiment, the contribution to d meas of the false-EDM effect in the mercury is about 50 times larger than the geometric-phase induced false-EDM effect of the UCNs.

B. False-EDM analysis: First iteration
Combining Eq. (8) with Eq. (19), the mercury's false-EDM contribution to d meas is seen to be where r B is the trap radius and the þ sign again corresponds to B 0 downwards. The notation d n;Hg;f is drawn from Eq. (87) of [4]. It follows that we can write where d 0 n is the true d n plus all other systematic effects discussed below, and R 0 0 is the value of R 0 where ∂B z =∂z ¼ 0. Equation (22) defines two straight lines, one with a positive slope for B 0 down and one with a negative slope for B 0 up. Naïvely, one would expect that the crossing point ðR × ; d × Þ would occur at ∂B z =∂z ¼ 0, and would therefore provide an estimator of d 0 n free of d n;Hg;f . In reality, as we discuss below, various effects can induce shifts in this crossing point, and we therefore must apply appropriate corrections.
As discussed above, Fig. 1 shows the data (binned for clarity) for d meas as a function of R 0 for each direction of B 0 . The solid red lines represent a least-squares fit to all 545 of the (unbinned) run results, using as free parameters the two intercepts c 1 , c 2 and a common absolute slope k. This was done by minimizing where the terms in the sum that correspond to B 0 down (up) used the þ (−) sign and the intercept c 1 (c 2 ). The fit yielded a crossing point at This fitted value of the slope was particularly influenced by one run (number 1900) that was taken at a large applied ∂B z =∂z (giving it significant leverage) and which departed from the fitted line by 2.7σ: excluding that run would have changed the value of the fitted slope to 1.63 AE 0.27 × 10 −26 e cm=ppm and increased the crossing-point EDM value by 0.22 × 10 −26 e cm, but there appeared to be no a priori reason to do so. The run was therefore included both in the initial 2006 analysis [1] and in the analysis that follows here. According to Eq. (21), the slope k is expected to be inversely proportional to Δh. We shall discuss the estimation of Δh in Sec. IV C. The relationship between R and d n;Hg;f is also affected substantially by the phenomenon of gravitationally enhanced depolarization [18], and this will be discussed in Sec. IV D below. In addition, the slope k can be altered by a few percent (although still remaining highly symmetric under B 0 reversal) by various other mechanisms including the UCNs' own false-EDM signal (a 2% effect); a slight reduction in mean free path due to cavities and grooves in the electrodes as well as to the presence of 10 −3 torr of He gas to prevent sparks [4,19] (1%); and a possible bias in the volume-averaged frequency measurement due to asymmetric surface relaxation of the Hg, e.g. if there is a preferential depolarization on one electrode or the other (up to 5%).

C. Estimating Δh
As seen in Eq. (21), the slope k of the lines in Fig. 1 is determined, to first order, by the height difference Δh between the centers of mass of the mercury and neutron ensembles. Initial estimates [1] This was based on a misinterpretation of the frequency response in a trap with a variable-height lid, arising from a lack of understanding at the time of the way in which gravitationally enhanced depolarization affects the measured neutron-to-mercury frequency ratio. A detailed Monte Carlo simulation has now been carried out, taking into account the known properties of the neutron source, guide tubes, elbow, polarizer foil, and storage chamber; it reveals that the spectrum was considerably softer, i.e. of lower mean energy, than had originally been assumed.
The simulated spectrum is also essentially consistent with a largely analytic calculation [3], which concluded that at the start of the storage period the UCN density ρ ϵ ðzÞ at the level z ¼ 0 of the lower electrode could be roughly approximated by a softened Maxwell spectrum, where V f ¼ 88 neV is the Fermi potential of the quartz insulator, above which energy ρ drops to zero. Henceforward, when we refer to the UCN energy, it is to be taken to mean the kinetic energy of the UCN at the level of the bottom electrode; or, equivalently, the total (kinetic plus potential) energy at any point in the bottle, where the potential energy is again referenced to the level of the bottom electrode. Equation (24) estimates relative densities at the bottom electrode, rather than total numbers in the bottle: the latter would have a slightly different distribution, since higher-energy neutrons extend further in height. In order to account for this, we note that for a monoenergetic group having energies between ϵ and ϵ þ dϵ, where the energy ϵ is given in terms of the maximum height that the neutrons could reach if traveling vertically upwards, the available phase space dictates that the density variation with height ρ ϵ ðzÞ has the form [20] and one must integrate this function over the bottle height (or, if less, the attainable height ϵ) to establish the relative numbers of UCNs of each energy trapped within the bottle. Although our best estimate of the spectrum is that given by the simulation-which will underlie the analysis that follows-we have carried out a complete analysis for both spectra, in order to see whether there was any significant sensitivity of the final result upon the initial spectrum. These two initial spectra are shown as the uppermost line and set of points in Fig. 11.
The spectrum softens during storage, since higherenergy UCNs not only have a higher interaction rate with the walls; they also have a higher loss probability per bounce, as given by [21] where f ¼ W=V is the ratio of the imaginary (W) to the real (V) part of the material's Fermi potential. In fact, although V is in general well known, W is difficult to determine; and in any case, losses are likely to be dominated by hydrogen that has diffused into the containing surfaces. This hydrogen can-because of its extremely high incoherent scattering cross section-substantially influence loss rates without significantly altering the surface potential V.
The value of f used in the simulation was therefore adjusted until, at f ¼ 4 × 10 −4 , the loss rate (including β decay) more or less matched the observed loss rate of neutrons in the bottle, as shown in Fig. 12. This is broadly similar to the value f ¼ 3 × 10 −4 obtained from an independent data-simulation comparison undertaken with the apparatus circa 2006 [22]. It was assumed in these calculations (as in [22]) that f would be similar for both the insulator and the electrodes. The softened-Maxwell spectrum was propagated in a similar manner (although for this f ¼ 3 × 10 −4 gave a better fit to the storage-time data, as shown also in Fig. 12). The resulting two spectra after 135 s of storage, by which time the Ramsey sequence is complete, are shown along with the initial spectra in Fig. 11. The resulting final simulated spectrum is slightly firmer than the final softened-Maxwell spectrum. Throughout this analysis the storage trap was modeled as a simple cylinder, with an additional cavity at the bottom center, 4.0 cm deep and 3.4 cm in radius, within which is set the UCN entrance door to the trap.
Given that the calculated Δh after storage is 3.7 mm, to first approximation the anticipated slope k of the lines of Fig. 1 should decrease from 1.57 × 10 −26 e cm=ppm (in [1]) to 1.18 × 10 −26 e cm=ppm, before taking gravitationally enhanced depolarization into account.

D. Gravitationally enhanced depolarization
Under the influence of gravity, neutrons of different energies effectively sample different regions of the storage trap. The (sometimes surprising) consequences of this effect have only recently been studied in detail [18,23], and have now been validated by comparison with data [24,25]. This stratification results in a dephasing of the ensemble in addition to that arising from the intrinsic depolarization that naturally occurs within each energy bin. Since lower-energy neutrons preferentially populate the bottom of the bottle, this dephasing is asymmetric, and results in a (strongly spectrum dependent) frequency shift. Calculations in [18] were based upon a Maxwellian velocity distribution with a sharp 93 cm height-equivalent energy cutoff. The simulated spectrum shown in Fig. 11 is much softer than this, showing a clear peak at 20-30 cm height-equivalent energy. Under these conditions the gravitational-depolarization effect distorts the shape of the crossing lines, significantly enhancing the slope for measurements beyond a few ppm from the origin, as shown explicitly in Fig. 13. The solid red line is the best linear fit to the data, mirroring the lines shown in Fig. 1, whereas the dashed green line is the expectation from Eq. (21). The dotted blue line includes the effect of gravitationally enhanced depolarization; we designate its functional form as The data of Fig. 1 were fitted to thus allowing the crossing-point coordinates R × , d × to vary, as well as allowing ξ to be multiplied by a factor k 0 as a consistency check of the input spectrum, since the slope of the function is completely determined by the spectrum via Δh. As well as crossing-point values of R × ¼ 2.62AE 0.96 ppm, d × ¼ ð−0.55 AE 1.51Þ × 10 −26 e cm, the fit yielded k 0 ¼ 1.15 AE 0.16; and, when it was repeated without inclusion of the aforementioned possibly anomalous datum of Run 1900, we found k 0 ¼ 0.93 AE 0.19. These numbers are consistent with unity, giving confidence in the simulated spectrum (for comparison, the softened-Maxwell spectrum yielded k 0 ¼ 1.37 AE 0. 19), and so, going forward, we set k 0 ¼ 1. We therefore fitted the data instead to This yielded crossing-point values of R × ¼ 2.11 AE 0.89 ppm and d × ¼ ð−0.59 AE 1.53Þ × 10 −26 e cm. χ 2 =ν had a value of 650=542 ¼ 1.19, unchanged from the earlier fit. We note that the change from three to two free parameters was accommodated almost entirely by shifting R × , with virtually no change in d × ; we deduce that the crossing-point value of the measured EDM d × is relatively insensitive to such changes in the slope, and therefore also to the detailed shape of the input spectrum. By eye it may appear that the straightforward linear fit to the data is so close to the predicted curve that it should make little difference to any analysis. However, the central linear region of the predicted curve is now shallower by a factor 1.18=1.88 ¼ 0.63 compared to the slope of the original linear fit, and any systematic shifts in R will therefore have a correspondingly reduced effectiveness in producing systematic shifts in d × .
A number of assumptions were made while calculating the size of the predicted false-EDM effect-for example, the gravitational-depolarization calculation here assumes only a uniform vertical B-field gradient; and no allowance has been made for the effect of grooves in the electrodes (which are necessary to locate the insulating trap walls). Effects of intrinsic depolarization were explicitly included in this model, although perhaps not perfectly since (as discussed in [18]) it has a significant dependence upon the specularity of reflections within the trap. However, its contribution to the frequency-shift effect is negligible [24]; primarily, it simply smears the phases. Additionally, it is worth noting that in principle the effective UCN spectrum could be slightly influenced, leading to different slopes k ↑ , k ↓ for the two B 0 field directions, if the spin relaxation were substantially different for the two configurations. This extra relaxation, if it existed, would be expected to be UCN-velocity dependent, and so would change the effective UCN energy spectrum. However, the observed polarization product α (i.e. the visibility of the Ramsey fringes) is identical at the 2% level for the two B-field polarities. This suggests [24] that the average UCN energy, and therefore Δh (which is approximately inversely proportional to the energy) should be the same to better than a few percent between the two polarities. It is difficult to conceive of a mechanism that would change the slope more significantly than this upon field reversal.
We have one further small correction to make at this point. Since the Ramsey-resonance measurements were carried out with an oscillating rather than a rotating transverse rf field, the measured frequencies are subject to a Ramsey-Bloch-Siegert (RBS) shift [13,14,26] given by FIG. 13 (color online). Plot showing the anticipated false EDM d n;Hg;f as a function of R 0 , for B 0 upwards. The dashed green line is the expectation from Eq. (21), and the blue dotted line shows the revised expectation when gravitational depolarization is taken into account. The red line is the linear best fit to the data, as shown in Fig. 1.
where ω 0 is the resonant frequency (182 rad/s), t is the period (2 s) for which each of the rf pulses is applied, and T is the (130 s) period of free precession between the rf pulses. This amounts to a mere 0.09 ppm.
Overall, therefore, we regard the match between the measured and expected curves to be in reasonable agreement. Going forward, we use the (RBS-corrected) value R × ¼ 2.02 AE 0.89 ppm, along with d × ¼ ð−0.59AE 1.53Þ × 10 −26 e cm, as input to our subsequent analysis.

E. Statistical sensitivity
For the majority of the data-taking runs, the upper electrode was held at voltages of circa AE80 kV, giving an applied electric field of E ¼ 7 kV=cm. (Voltages as high as 130 kV were sustainable towards the end, whereas in early runs it was not possible to go much above 60 kV.) The average polarization was α ¼ 0.58, and the Ramsey coherence time was T ¼ 130 s. Removing from consideration the cycles measured at E ¼ 0 kV=cm, Eq. (13) yields an anticipated sensitivity of 1.34 × 10 −26 e cm. Our achieved sensitivity is a few percent larger than this, for a number of reasons: for example, field changes meant that not all cycles were taken at the intended points on the Ramsey curve; some cycles had frequency errors enlarged by relatively rapid depolarization of the mercury; and any changes in field gradient resulting in changes in R 0 during the run (which at the ∼0.2 ppm level or below would not be detectable) could also add noise at the percent level.

V. MODIFICATIONS TO THE FALSE EDM SIGNAL
As indicated above, there are some processes that can displace the crossing point, and thus interfere with this technique of removal of false-EDM effects-essentially any process that changes R a and/or d n;Hg;f without conforming to the ratio between the two given by Eq. (21), and where, in addition, the changes differ with the direction of B 0 . These processes may broadly be divided into two categories: those that shift frequencies, and hence move the lines of Fig. 1 horizontally, and those that move the lines vertically by altering the gradient-induced false EDM signal relative to the expectation from h∂B z =∂zi. In this section we summarize these effects, before describing in more detail in Sec. VI the additional diagnostic measurements that were carried out first to characterize and then to compensate for them. Other systematic effects, which contribute to the overall uncertainty but that do not introduce a bias to the data, are considered in Sec. VIII.
A. Uncertainty in γ n =γ Hg An offset in γ n =γ Hg shifts the two lines of Fig. 1 sideways in the same direction by the same amount, leaving d meas unaffected. This therefore has no direct effect on the final result, although it can influence slightly the overall fit to the data.

B. Horizontal quadrupole fields
Fields with finite ∂B x =∂y and/or ∂B y =∂x but with (∂B x =∂x þ ∂B y =∂yÞ ¼ 0 ¼ −∂B z =∂z cause R a to increase quadratically [4] without contributing to d n;Hg;f . This arises because the neutron spins follow the total field direction adiabatically, with a correspondingly increased precession frequency, whereas the (nonadiabatic) Hg atoms average out the transverse components and are generally sensitive only to the B z component. We consider here a quadrupole field aligned with z, with B x ¼ qy, B y ¼ qx. This gives rise to a shift in R a of where, as before, r is the radius of the trap. As discussed in [3], the region in the vicinity of the EDM measurement volume was scanned with a fluxgate magnetometer. This revealed the presence of such horizontal field components, quadrupolar in form and a few nT in magnitude throughout the 12 cm height of the measurement bottle. The magnitude of these quadrupole fields also depends slightly upon the orientation (upwards or downwards) of the B 0 holding field. This is to be expected, since the B 0 return flux passes through the innermost mu-metal shield: indeed, about 40% of the B 0 field arises not from the B 0 coil directly but from the magnetization of the shields. The geometry of the 1.15 m diameter mu-metal shield around the B 0 coil has irregularities on the scale of a few millimeters, so it seems reasonable to expect that aberrations in the B 0 field near the outer edges of the trap should be of the order of a few nT. Unlike the remanent fields arising from any permanent magnetization of the shields themselves, any B 0xy fields will reverse when B 0 is reversed. Contributions from these two different sources may therefore add in one field direction and subtract in the other, leading to a differential shift of all R a values, and thus of the two lines, thereby changing the EDM value d × at the crossing point of Fig. 1.
Further evidence for the existence of such quadrupole fields came some years after the EDM measurements described in this paper, when the shields were disassembled, shipped to the Paul Scherrer Institut (PSI) and rebuilt. The results of a detailed fluxgate scan taken in 2010 within the rebuilt shields is shown in Fig. 14, and here a quadrupole-type field pattern can be seen clearly. A detailed harmonic analysis of the field shape yielded an average q value of −1.29 ð3.24Þ nT=m for B 0 up (down), which would imply corresponding quadrupole shifts of δR Q ¼ 0.02 ð0.14Þ ppm for the two respective B 0 directions. We do not include these figures in our analysis at this point, but we will return to the issue briefly during our conclusions.

C. Rotation of the Earth
The rotation of the Earth shifts all of the frequency ratio measurements R a to lower values by 1.33 ppm when the B 0 field is upwards, and to higher values by 1.33 ppm when the B 0 field is downwards [2,5]. This therefore acts in a manner very similar to that of the horizontal-quadrupole shifts.

D. Localized losses
B-field averaging in the trap is affected by localized loss of UCN and Hg particles, and by polarization loss in the presence of the 10 −3 fractional B 0 inhomogeneities, which may change with B 0 direction. However, we estimate that the resulting R a shifts are < 0.1 ppm and < 0.01 ppm for the UCN and Hg respectively. Furthermore, they will be indistinguishable from the quadrupole shifts. We do not consider them further.

E. Magnetic dipoles
The field of a permanent magnetic dipole (PMD) close to the trap results in a nonuniform ∂B z =∂z, and therefore induces a frequency-ratio shift δR dip , moving the lines of Fig. 1 horizontally by equal and opposite amounts for the two directions of B 0 . As demonstrated by Eq. (20), it also adds to d n;Hg;f a false EDM d dip [27], thus shifting the lines vertically by this amount.
Our fluxgate magnetometer surveys of the trap cannot rule out PMD fields of less than 1 nT at 2 cm from the inner surface. Large areas of the trap are SiO 2 or Al, backed by large voids, and do not come under suspicion; but the mercury and UCN doors involve a heterogeneous collection of small parts close to the trap.
Based upon our simulations [4,27] (since confirmed to be perfectly consistent with the analytic prediction of [17]), we included in our 2006 analysis a d dip uncertainty of AE6.0 × 10 −27 e cm to allow for an undetected 1 nT PMD at the mercury door.
In the case of the UCN door we have better diagnostics. As discussed above, this door sits at the bottom of a small cavity, 4.0 cm deep and 6.8 cm in diameter, at the center of the lower electrode. As detailed in Sec. VI, measurements were made of the neutron and mercury response in a storage trap with a ceiling of variable height, and the results provided clear evidence for the existence of both dipole and quadrupole fields. This conclusion was reinforced by measurements undertaken some years later when, as part of the process of characterizing the system during its removal to PSI, scans were made of the components. A crank-andsliding-plate assembly that converted rotation of the driving shaft into the linear motion of the neutron door was revealed to be slightly magnetic, with summary records showing that a sensor several centimeters away registered a field of 2.5 nT. A disk that constituted the lower part of the door itself also yielded a total field of about 100 pT at a distance of ∼13 cm from the sensor. Various other parts, designated "screws," "piston," "Al frame" and so on, were tabulated as having associated fields of a few pT to a few tens of pT. Detailed records of orientations and positions were not kept, since the primary goal at that time was merely to determine whether or not the components were magnetic. It is therefore not possible at this point to specify the dipole position and strength precisely from the information now available, so as in the 2006 analysis we rely upon the contemporary measurements and regard the later scans merely as qualitative supporting evidence, but it seems impossible to avoid the conclusion that there were fields of ∼ nT strength and approximately dipolar in nature in the vicinity of the UCN door. The picture is also consistent with measurements undertaken (after the EDM measurements, but prior to moving the apparatus FIG. 14 (color online). Results of a detailed scan of the magnetic field inside the shields, taken some years after the EDM measurements described in this paper. By this time the shields had been disassembled, moved to PSI and reassembled. The quadrupole-field pattern is clear, and is obviously associated with the shields themselves. The scales on the x and y axes represent the horizontal coordinates, in millimeters. The magnitude of the field is represented both by the lengths of the arrows and, equivalently, by their colors, as shown in the scale (in nT) at the bottom.
to PSI) with a fluxgate magnetometer located in the vertical neutron guide immediately underneath the neutron door. This registered a slight change (circa 1 nT) in the ambient magnetic field when the door mechanism was operated.
The dipole-field analysis that we have carried out assumes a vertically orientated dipole centered on the axis. There is no particular reason to assume that any PMD would have this particular orientation or precise radial position. However, we note that the associated frequency shifts arise via h∂B z =∂zi, which would be much lower in the case of a horizontal orientation (not least due to cancellation of the respective field components from opposite ends of the dipole). Furthermore, the analysis of [17] (and in particular Fig. 2 therein) shows both that a horizontal alignment (for a dipole within a few centimeters of the axis) yields a much smaller EDM contribution than does a vertical alignment, and also that the EDM contribution (of the vertically aligned dipole, again within a few centimeters of the axis) has very little sensitivity to the radial position. We therefore conclude that our model accounts in a reasonable manner for the most important contributions.
The post hoc magnetic-field scans at PSI revealed no magnetic contamination in the vicinity of the mercury door, although a rescan in 2011 showed a small contaminant outside the storage volume, at a radius of 31 cm, within the threaded hole connecting the ground corona ring to the lower electrode. If present throughout our measurements, this could have contributed a false EDM of up to AE4.0 × 10 −28 ecm (according to the Pignol-Roccia analysis). In order to be conservative, though, we retain the full AE6.0 × 10 −27 e cm contribution to the systematic uncertainty, in order to accommodate both the possible undetected presence of a (further) PMD and the fact that we have not modeled the actual door dipole perfectly.

F. Mercury light shift
The presence of the mercury reading light can, via the RBS mechanism [13,14], shift the resonant frequency of the mercury atoms [28,29]. Such shifts are produced by any small component, parallel to B 0 , of the 204 Hg probe light beam passing through the precessing 199 Hg atoms. This component, and the consequent R a shift, reverse sign on reversal of B 0 . False EDM signals can thenarise in two distinct ways: (i) Direct: If any changes in intensity are correlated with the electric field direction, the resultant frequency shift would mimic an EDM. (ii) Indirect: Even if the intensity of the light is completely constant, the light shifts move the data points of the two lines in Fig. 1 in opposite (horizontal) directions, and thus contribute a vertical offset to the crossing point. A slight dependence of R a on the incident light intensity was indeed found during the 2006 analysis, the magnitude ∼0.2 ppm being in agreement with theory. (Similarly small shifts were found in a recent measurement of the ratio of the neutron to 199 Hg gyromagnetic ratios, using the upgraded nEDM apparatus at PSI with similar Hg discharge bulbs but with a new storage chamber [9].) A correction to d meas was made on a run-by-run basis, leading to an overall correction of ð3.5 AE 0.8Þ × 10 −27 e cm. For this analysis we have simply applied exactly the same set of (almost imperceptible) shifts to the data, prior to the crossing-lines fitting procedure. Since the dependence should thereby have been removed, we expect no remaining net bias from this source. Details of the light-shift analysis are discussed extensively in [3], and will not be addressed further here.

VI. AUXILIARY MEASUREMENTS
Separate ν n , ν Hg and R a measurements (without E fields) were made in an auxiliary trap with a roof that could be raised or lowered to change the height H. Critically, the trap was built on the same lower electrode and door mechanism that were used for EDM data taking. Full details of the apparatus and of the measurements made are available in [30].
The trap had a smaller radius than the main data-taking bottle, 18.5 cm rather than 23.5 cm, and because the floor of the trap consisted of an aluminum plate that rested on the lower electrode, the door-cavity depth became 6.0 rather than 4.0 cm. These differences would have reduced the quadrupole shifts by a factor Q ∼ 0.59. (Note that, although the quadrupole shifts are reduced in this way, the shifts due to Earth's rotation are not.) The roof and floor were coated with deuterated polystyrene (DPS) rather than the diamondlike carbon used for the main data-taking trap, but the DPS tended to flake away, leaving some bare aluminum exposed: the resulting UCN spectrum would therefore have had a slightly lower end-point cutoff than in the data-taking trap, although after some tens of seconds of storage this would have made little difference. Measurements were taken at storage times of several tens of seconds' duration rather than the full 135 s, because with the reduced-height bottle the loss rates were much higher. Spectrum calculations were carried out exactly as described above for the data-taking bottle, assuming the same initial spectrum and per-collision loss rates. Appropriate Ramsey-Bloch-Siegert frequency-shift corrections [13,14], typically amounting to 0.1-0.2 ppm, were applied to each measurement.
Assuming B z ðzÞ ¼ b 0 þ b 1 z þ b 2 z 2 , one can show that h∂B z =∂zi V ¼ 0 for a trap height H when ν Hg is the same for roof settings at H=2 and H. This situation was approximated by adjusting the trim coils until the same ν Hg was obtained for trap heights of 109 and 59 mm. With this field established, R 0 was measured for trap heights H of 34, 59, 84 and 109 mm for B 0 down. The procedure was repeated for B 0 up, although in this case measurements were only made for trap heights of 59 and 109 mm. We denote the values of R 0 in this "gradient-compensated" environment as R 0 gc ðHÞ. Naturally, it was not possible in reality to trim the gradients perfectly, but studying the Hg response enabled appropriate corrections to be made. For this analysis we have restricted ourselves to the 11 measurement runs for which such gradient corrections to R 0 amount to less than 0.5 ppm.
The resulting values of R 0 gc ðHÞ are listed in Table I. Their averages for each height and magnetic-field configuration are shown in Fig. 15. The original intention of this series of measurements had been to determine a value for the ratio of gyromagnetic ratios γ n =γ Hg . However, the strong variation with height, and in particular the separation of B 0 up and down values as the height was reduced, led to the conclusion that there was a dipole field of strength ∼1 nT penetrating into the door cavity. Such a field could in principle have been produced by a pointlike contaminant more or less at the limit of detectability by a fluxgate magnetometer, such as the aforementioned door-actuator crank. It could also in principle have arisen from a more dispersed (and therefore undetectable) source such as, for example, a slight residual magnetization within the BeCu door plate or from a thermoelectric current as the door mechanism operates in vacuo.

VII. GLOBAL FIT
In this section we detail the manner in which each contributing component was modeled and incorporated into the fit, before moving on to discuss the outcome of the fit and the consequent corrections that must be applied to the crossing-point EDM measurement d x .

A. Dipole field
A simple model incorporating a dipole of moment p situated a distance z dip below the surface of the door was found to fit the form of the auxiliary-trap data extremely well. The average z component hB z i of such a field over a circle of radius r a distance z above the dipole is By dividing the trap volume numerically into 1 mm thick discs, the additional contribution to B z from the dipole field and from an additional uniform applied gradient could be calculated as a function of the height z. Assuming a spectral distribution of UCN as given in Sec. IV C (in contrast to the 2006 analysis, which assumed a Maxwellian velocity distribution up to the quartz Fermi-potential cutoff), the relative height distribution of the neutrons in each energy bin was calculated. The contributions to the mercury and UCN frequencies, and hence the expected ratio R a , could then be calculated as a function of the height H of the trap roof. Each iteration of the fit used its selected values for the strength and position of the dipole, for R γ , and for the quadrupole shifts. An appropriate uniform gradient was then applied so as to yield the same calculated mercury response at H ¼ 59 mm and 109 mm, thus matching the experimental configuration of the auxiliary bottle. The R a values for each bottle height H were then calculated and compared with the data before proceeding to the subsequent iteration of the fit. We designate the shifts in R a due to the gradient-compensated dipole field as δR dip ðHÞ.

B. Quadrupole fields
The quadrupole-field shifts δR Q↑ , δR Q↓ (where the subscript arrows indicate the respective B 0 directions) are assumed to be independent of height (in reasonable agreement with the fluxgate measurements), and to be  15 (color online). Values of R 0 (in ppm) for a "gradientcompensated" field (see text) for each B 0 -field polarity, as a function of the height of the storage volume. Where more than one datum was taken for a particular configuration, the point shown represents the average. The curved dotted lines represent an approximation to the fitted function: they are based upon the spectrum calculated for 75 s of storage in a bottle of height 84 mm, whereas the fit uses spectra specific to each datum. proportional to a scale factor Q that depends in principle only upon the (volume-averaged) square of the trap radius, thus yielding Q ¼ 0.59 for the auxiliary trap relative to the main data-taking trap as discussed in Sec. VI above. As noted above, these quadrupole-type fields almost certainly arise from the mu-metal shields. It is possible in principle that the dependence upon radius may not be as expected, particularly towards the outside of the data-taking bottle, which extends to about 40% of the (58 cm) radius of the inner mu-metal shield. As a test, we allowed the Q values for the auxiliary trap to be reduced by a further factor f Q , which was introduced as a free parameter in the fit. Ultimately this made no difference, as the fitted value of f Q ¼ 0.95 was close enough to unity that the final EDM limit was unaltered. We therefore omitted this from the final analysis; i.e. we set f Q ¼ 1.
The quadrupole-field contributions δR Q to the R 0 values within the fit were therefore Q · δR Q↑;↓ for the auxiliary bottle and δR Q↑;↓ for the main data-taking bottle.

C. Earth's rotation
As noted above, the rotation of the Earth shifts all of the frequency ratio measurements R a to lower values by δR E ¼ 1.33 ppm when the B 0 field is upwards, and to higher values by the same amount when the B 0 field is downwards. For either the auxiliary or the data-taking trap in isolation this effect would be indistinguishable from the differential quadrupole shifts, and in the 2006 analysis it was not taken into account separately. However, when extrapolating from one trap to the other a correction has to be made since the Earth's rotation shift is independent of the trap radius and therefore of Q. The shifts δR E are therefore built explicitly into the fit in this analysis.

D. Literature γ ratios
With no applied gradient and with any dipole and quadrupole contributions appropriately accounted for, the true value R γ of R 0 should be recovered. The average of the two independent literature values of jγ n =γ Hg j, namely 3.8424574(30) from [9] and 3.8424560(66) derived from [6] and [7], is 3.8424572 (27). This was used as an input estimator of R γ . It is equivalent to R 0 ¼ 0.31 AE 0.71 ppm.

E. Crossing point
By symmetry, the average quadrupole shift δR Q ¼ ðδR Q↑ þ δR Q↓ Þ=2 away from the true value of R γ should correspond to the crossing point R × of the lines of Fig. 1: It is convenient at this point to introduce the parameter which represents the quadrupole-field splitting.
We note that, as shields were opened, closed and demagnetized many times throughout the four years of data taking, the profile of the field could in principle have changed during this time-indeed, the same trim-coil settings could sometimes yield variations in R a values of AE3 ppm. However, since the overall χ 2 ν of 1.2 for the fit to the lines of Fig. 1 was not excessively large, such changes were clearly dominated by alterations in the vertical gradient, which would result in movement along the lines rather than displacement of data points from the lines. This result is not surprising-the largest holes in the magnetic shields were on top and underneath, and this would have provided the primary access route for ingress of external fields.
We account for these random variations by increasing the uncertainty of the datum for the crossing point R × by a factor ffiffiffiffi ffi χ 2 ν p ¼ 1.09 in the global fit.

F. Fitting procedure and outcome
The rather simple underlying model of a dipole and a quadrupole field yields a fairly complex functional shape that can only be evaluated numerically. The (five) parameters that were allowed to vary were the dipole strength and z position, R γ , the quadrupole shift δR ↓ , and δR QQ . The (13) input data points were (a) the literature value of R γ , (b) the crossing point R × from the earlier fit to the crossing lines, and (c) the 11 auxiliary-bottle measurements R 0 gc of R 0 in the gradient-compensated fields at various bottle heights H and storage times t s : The fit used the Levenberg-Marquardt algorithm. The resulting best-fit function is represented, together with the data, in Fig. 15. For clarity, the data in that plot are reduced to six points by averaging over each height and B 0 up/down configuration, and the plotted function was generated with the spectrum appropriate to a 75 s storage time in an 8.4 cm high bottle.
The fit yielded output values of R γ ¼ −0.22 AE 0.64 ppm, quadrupole-field shifts of δR Q↓ ¼ −0.5 AE 1.1 ppm and δR Q↑ ¼ 2.7 AE 1.1 ppm, and a downwards-pointing dipole of moment m ¼ 4.2 × 10 −7 Am 2 situated 1.3 cm below the surface of the door-quite consistent with its being a part of the door mechanism. Such a dipole would contribute a field 1.1 nT in magnitude 2 cm above the surface of the door and 0.3 nT at the level of the electrode surface, and as such would be at or below the limit of detectability. In fact, the fit was fairly insensitive to the dipole strength and position, with different starting values of the parameters yielding a range of dipole strengths from 2 to 8 × 10 −7 Am 2 correlated with a range of positions from 0.8 to 2.2 cm below the door. R γ , R × and δR QQ , on the other hand, were robust.
We note in passing that the value of R × ¼ 0.86 AE 1.08 ppm that emerges is significantly shifted (approximately 1σ) from the input value of 2.02 AE 0.97 ppm. However, since d × is essentially independent of such small changes in R × , as discussed above, we do not pursue this line of inquiry further-what we are interested in is shifts relative to R × .
The overall χ 2 of the global fit was 8.9 for the 8 degrees of freedom. The quality of the fit is therefore rather good, particularly bearing in mind the complexity of the fitting function, and suggests that the model is both sound and complete-there is no evidence of contamination from fields of a different form.

G. Field-shift corrections
In the absence of the dipole field, we would expect the R 0 values at which the vertical magnetic-field gradient within the data-taking bottle is zero to be determined by the respective quadrupole-field offsets for each of the two B 0 directions, together with the shifts due to Earth's rotation: where the upper (lower) sign corresponds to B 0 down (up). It would therefore be appropriate under those circumstances to determine the EDM d 0 at either of these points rather than at the crossing point (R × , d × ). By symmetry, the values R 0 either side of the crossing point, i.e. at R 0 ¼ 1.14, 0.58 ppm respectively. They therefore share a common EDM value d 0 . Using our fitted function from Sec. IV D, we find that this shift in R 0 away from R × requires us to add δd ¼ ð−0.33 AE 0.14Þ × 10 −26 e cm to our earlier crossing-point value d × ¼ ð−0.59 AE 1.67Þ × 10 −26 e cm (where we have expanded the uncertainty on d × by ffiffiffiffi ffi χ 2 ν p as we did for R × ), yielding d 0 ¼ ð−0.92 AE 1.68Þ × 10 −26 e cm as the EDM value corrected for quadrupole shifts and Earth's rotation.
Further corrections must now be made to accommodate the dipole field. The formulation of Pignol and Roccia [17] predicts that the dipole characterized by the parameters of our fit would produce, in our storage bottle with its door cavity, a false EDM of d dip ¼ 0.44 × 10 −26 e cm. The crossing lines of Fig. 1 are therefore higher by this amount than they would otherwise be, and we must compensate by subtracting d dip from all of the measured data, and thus, ultimately, from d × .
In addition, this dipole field applies a volume-averaged magnetic-field gradient ∂B z =∂z of 0.26 nT/m, which shifts the normalized frequency ratio R 0 by δR dg ¼ AE0.98 ppm (where the subscript "dg" indicates "dipole gradient," and the upper sign once again corresponds to B 0 down). This shift of the lines in opposite directions moves the crossing point downwards by 1.15 × 10 −26 e cm. The net dipolerelated correction that must be applied is therefore now ð0.71 AE 0.07Þ × 10 −26 e cm. As mentioned above, the dipole strength and position arising from the fit are strongly correlated, so the uncertainty on this total dipole-related shift was calculated by studying the behavior of χ 2 as a function of the shift for a variety of dipole strengths and positions.
We also at this point add in quadrature to the uncertainty the contribution of 0.6 × 10 −26 e cm from a possible mercury-door PMD, as discussed in Sec. V E.
Taking all of these dipole-related components into account, the net EDM thus far-compensated for the false-EDM effects arising from field gradients, dipole and quadrupole fields, and Earth's rotation-is d fec ¼ ð−0.21 AE 1.79Þ × 10 −26 e cm, where the subscript "fec" stands for "field-effect compensated."

H. Consistency check: Polarization data
Because UCNs of different energies have different values of Δh, the surviving UCN polarization within a trap decreases in proportion to ð∂B z =∂zÞ 2 . One would therefore expect that α would be maximized close to the R 0 value (for each B 0 direction) at which h∂B z =∂zi is zero, as given by Eq. (36) with an additional correction for the dipole shift δR dg . In fact, our calculations show that, due to the dipole field, ð∂B z =∂zÞ 2 is minimized at 0.920 ppm below R 0 0↑ for B 0 up, and 0.920 ppm above R 0 0↓ for B 0 down. We therefore expect the polarization α to be maximized at R 0 ¼ 0.2 AE 1.1 ppm, 1.5 AE 1.1 ppm for B 0 up, down respectively.
During runs prior to January 2000 the gradient was varied frequently and over a fairly wide range. We have plotted in Fig. 16 the polarization α as a function of R 0 for this subset of runs, which represent 14% of the total statistical weight of the EDM data analyzed here. Although other factors (e.g. loss of polarization during filling) affect the data to a greater or lesser extent, the envelope of these data points ought to peak in the region where the magneticfield gradient is minimized.
ln the 2006 analysis, bands of points along the tops of these distributions were fitted to Gaussians to find the positions of the peaks. The fit was carried out multiple times, each time repeatedly removing points that lay more than pσ below the curve, where p was varied between 4.0 and 1.0, such that the final remaining points would represent the envelope. In recognition of other factors contributing to depolarization, and in order not to allow the few data points with particularly small uncertainties on α to skew the fits unduly, an additional uncertainty of 1% was added in quadrature to each error bar. Points with error bars are those used in the fits represented by the solid lines.
Although these results were used in the global fit in the 2006 analysis, we feel upon reflection that the large scatter of data points underneath each peak casts doubt upon their validity for this purpose. In addition, the fitted peak position has some dependence both upon the number of data points included and upon the functional form of the fitted curve-which, as we now understand [24], can resemble a sharp peak rather than a Gaussian. We therefore present them here only as a consistency check, and note that, although the R 0 values at which the fitted curves peak appear to be arguably just a little higher than expected, they are close enough to be regarded as being in broad agreement.

VIII. OTHER SYSTEMATIC ERRORS
We now consider systematic errors that do not involve the field-gradient induced false-EDM effect, and that do not bias the results. For those that are discussed in detail in [3], we summarize the nature of the effect and state their contributions to the overall uncertainty. We begin with firstorder v × E effects.

A. v × E effects
The Lorentz transformation of electric and magnetic fields to a moving reference frame is [12] where γ ¼ 1= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 − v 2 =c 2 p . As noted above [Eq. (17)], in the context of the false-EDM effect, a particle moving fairly slowly (γ ≈ 1) through an electric fieldẼ therefore experiences an additional magnetic field above and beyond the laboratory magnetic fieldB 0 . Such fields are clearly linked to the electric field, and the necessity of controlling the consequent spurious effects provides a strong constraint on the design of the experiments. In particular, to generate a systematic error in the EDM requires thatṽ,Ẽ andB 0 have components that are mutually perpendicular. It was therefore necessary in this experiment to keepẼ andB 0 closely aligned and the velocityṽ averaged over the storage time as small as possible. Only the motion of the center of mass of the UCN gas and any net rotation of the gas about the center of mass can contribute to this first-orderṽ ×Ẽ effect. The effect has been clearly seen in atomic beam experiments [31][32][33] and has been a cause for concern in neutron beam experiments.
If the neutrons have an average center-of-mass velocity of η perpendicular toẼ, andẼ has a component ϵE perpendicular toB 0 , theṽ ×Ẽ component of the magnetic field will give a systematic error in the EDM of jdj ¼ μ n c 2 ηϵ ð40Þ The only source of indefinitely sustained motion that can give rise to a finite average velocity would be a net upwards movement of the center of mass due to warming of the neutron gas as a result of inelastic collisions with the walls. The effects of such collisions, in which the neutron energy changes by only a small amount so that it remains in the UCN range, have been looked for by Richardson [30]. For a greased surface, an upper limit of 50 feV for the energy change per wall collision was established, from which it can be calculated that the maximum movement will be no more than 1 mm during the 130 s Ramsey measurement period. If the volume-averaged angles between E, B and v are each as high as 2°¼ 0.035 radians, the induced false EDM will be 1 × 10 −30 e cm. (Note that such warming is distinct from the changing Δh arising from the preferential loss of faster UCNs: the latter do not ultimately contribute to the EDM measurement, and so cannot generate a systematic error.) If there is a difference in a transverse coordinate x between where the neutrons enter the trap and their average center of mass during the storage time, the net average velocity is η ¼ x=T s . The effect of any such motion will be greatly reduced by the fact that it is only motion during the storage time between the two oscillating fields that affects the precession frequency. Since the trap is filled for 20 s, which in our apparatus is roughly 1.3 filling-time constants, and this is followed by a 5 s settling period, the stored neutrons have spent an average of one filling-time constant in the trap before the first oscillating field is applied. Collisions with the wall are expected to destroy any nonrandom motion within 100 bounces, i.e., about 10 s, so that any such motion that results from the filling process will mostly have died away before the first oscillating field is applied. The position of the center of mass during storage is expected to be very close to the axis of the 67 mm diameter guide tube. In the unlikely event that the neutrons fill preferentially at x ¼ 1 cm towards one side of the guide tube, and ignoring any damping of the motion in the trap during filling, we would be left with a net motion during the 130 s Ramsey measurement of up to ∼5 mm. This effect would then contribute 2ϵ × 10 −28 e cm ¼ 6 × 10 −30 e cm. This dominates the earlier result associated with warming, and therefore becomes our total uncertainty for first-order translational v × E; to be conservative we round the value up to 1 × 10 −29 e cm.
Higher-order v × E effects, in particular from the net increase in jB 0 j as the v × E component is added to it in quadrature [34], contribute < 10 −30 e cm.
In a similar manner to the translational effect, any net rotational flow of the UCNs in conjunction with a radial component of the E field may lead to an induced EDM signal. However, any such flow of UCNs is expected to be attenuated by wall collisions before the first Ramsey pulse is applied. This systematic effect requires the electric field to have a radial component ϵE 0 , such as may exist near the outside of the trap if the insulator of the trap becomes charged. If a net fraction f rot of the neutrons has such a flow with a velocity of 1 m/s, and if the flow persists for 10 s after the storage period has begun, and hence 5 s into the Ramsey measurement, the neutrons will have moved a net ∼5 m through this field within the Ramsey measurement. The effect on the EDM is then 5f rot ϵ × 10 −24 e cm. It is reasonable to assume that only a very small fraction such as f rot ≤ 0.001 could persist in orbits through this anomalous field region, and the radial field fraction ϵ is unlikely to be larger than 10% of the primary field, giving a net error from this source of below 5 × 10 −28 e cm.
There will also be additional cancellations to the v × E effect from the reversal of B 0 , but in order to be conservative we do not consider these.

B. Uncompensated magnetic-field fluctuations
There may in principle be residual effects from B field fluctuations, such as hysteresis in the mu-metal shield following disturbances in the stabilized B 0 coil current supply caused by pickup from the high-voltage changes. As discussed in [3], we would expect this to manifest itself as an approximately dipolar field B d originating at the HV feedthrough. Because of the effective height difference Δh between the centers of the neutron and mercury systems, the experiment retains a slight vulnerability to such a systematic. The shielding factor is expected to be about where r ≈ 55 cm is the distance from the source of the field to the center of the trap. In our 2006 analysis, and in [3], we used Δh ¼ 2.8 mm, and thereby concluded that the protection factor was about 66. In light of our revised calculations of the spectrum, which yield Δh ¼ 3.7 mm, we now reduce this to a factor of 50. Figure 18 of [3] showed the apparent EDM signals of the neutron and mercury channels individually, plotted against one another. Figure 17 here shows the same data, but this time binned for clarity.
The neutrons alone yielded a net uncompensated EDM signal of ð17 AE 4Þ × 10 −26 e cm. We therefore expect the mercury-magnetometer compensation to shield us from this systematic effect to a level of ð17 × 10 −26 Þ=50 ¼ 3.4 × 10 −27 e cm.

C. Mercury atom EDM
The measured EDM values d meas are obtained from a linear fit to the ratio ν n =ν Hg versus E. In principle therefore, d meas contains a contribution from the intrinsic EDM d Hg of the 199 Hg atom. The true d Hg has been shown FIG. 17 (color online). After [3], Fig. 18. Apparent neutron EDM signals (due to uncompensated random magnetic-field fluctuations) as a function of the corresponding apparent mercury EDM signals. Data are binned for clarity. to be ð0.49 AE 1.29 stat AE 0.76 syst Þ × 10 −29 e cm [35], so the systematic error thereby introduced into d meas is a negligibly small ð−2 AE 6Þ × 10 −29 e cm.

D. Electric forces
Another possible source of systematic error arises from electrostatic forces, which may move the electrodes slightly. Discussion of this issue in [3] concluded that the associated systematic uncertainty is 0.4 × 10 −27 e cm.

E. Leakage currents
If the leakage current that flows through (or along the surface of) the insulator between the electrodes has an azimuthal component, a component of the magnetic field due to the current would be parallel (or antiparallel) toB 0 and would produce a frequency shift that changes sign when the polarity of the electric field is reversed, giving rise to a systematic error in the EDM. As discussed in [3], the false signal that would result is likely to be no larger than 0.1 × 10 −27 e cm. As is clear from Fig. 18, which is drawn from [3], no dependence of frequency shift upon the leakage current is apparent in the data.

F. Sparks
High-voltage breakdown within the apparatus can produce localized high current densities, and can in principle lead to permanent changes in the residual magnetization of the magnetic shield. If these sparks occur preferentially for one direction of the electric field, systematic effects could be induced in the precession frequency, and hence in the EDM signal. A similar effect could be produced by hysteresis in the innermost shield following a disturbance to the B 0 supply. However, the mercury would naturally compensate for any such effect, just as with any other shifts in the magnetic field.
As sparks invariably disrupt the mercury frequency measurement, cycles that contain them are excluded from the analysis, so beyond the residual effects just discussed the sparks themselves cannot contribute to any artificial EDM signals. We therefore do not associate any additional systematic uncertainty to this effect.

G. HV AC ripple
A "ripple" on the high voltage would generate an oscillating displacement current in the storage chamber and thereby an oscillating B field. Through the Ramsey-Bloch-Siegert mechanism [13], this could in principle lead to changes in precession frequency. A detailed discussion of this effect in [3] concludes that such a mechanism would in this case not contribute to the overall uncertainty at a level higher than 1 × 10 −29 e cm.

H. Effect of nonuniform mercury depolarization
The detrimental effect of the high voltage upon the mercury depolarization time could result in a false signal if (a) the average depolarization time were different for the two HV polarities, and (b) the mercury frequency had some small dependence upon the depolarization time. As discussed in [3], evidence for such an effect was sought in the data, and it was ruled out at the level of 1.2 × 10 −29 e cm.

I. Artifacts of the measurement process
In this class of systematics there is usually no real change in precession frequency-only an apparent change due to a malfunction of the measurement process. An example might be bursts of false counts in the neutron detection channel derived from a high-voltage spark. No evidence of any such effects has been seen, and in any case any artificial signals from such effects would tend to cancel upon reversal of B 0 . We do not associate any further uncertainty with this mechanism.

J. Accuracy of Hg frequency measurements
A number of mechanisms can affect the frequency measurement of the Hg magnetometer. These are discussed extensively in [3], and are already implicitly dealt with by our analysis. No further uncertainty is assigned to any of these effects.

K. Stability of results
As a final check for possible unforeseen systematic errors, we consider the stability of the measurements. On a run-by-run basis, one can consider the stretch values [Eq. (15)] from the geometric-phase line fit of Fig. 1 as a function of time, as shown in Fig. 19: this removes the natural variation in the measured EDM due to the changing R a values. If there were an offset in one series of data due, for example, to a magnetic anomaly that was later removed, one would see points typically above the line for a series of

L. Systematic error summary
Table II lists all of the systematic errors that we have discussed, and gives their uncertainties. Its last line specifies the net shift in EDM value, and the total systematic uncertainty. For the avoidance of doubt, we note that the shifts listed represent the offsets generated by the biases in question. The compensating corrections that have been applied are therefore equal in magnitude and opposite in sign to the stated values.

IX. RESULTS AND CONCLUSION
We have reanalyzed the data that were used in 2006 to calculate an upper limit on the magnitude of the neutron electric dipole moment [1]. As in the 2006 analysis, earlier data from the interim result of [36] have been excluded. 1 We have here taken into account recent developments, including an understanding of the process of gravitational depolarization and more detailed calculations of the spectrum of the stored UCNs, and we have also taken into consideration measurements made in the meantime on the apparatus. The picture that emerges appears to be perfectly self-consistent.
Once the crossing-point EDM value d × ¼ ð−0.59 AE 1.53Þ × 10 −26 e cm is corrected for the systematic biases listed in Table II, as shown explicitly in Table III for each stage of the analysis, we obtain (by adding in quadrature the statistical and systematic uncertainties) a final value of d n ¼ ð−0.21 AE 1.82Þ × 10 −26 e cm. The 90% and 95% confidence limit (C.L.) ranges are therefore −3.2 < d n < 2.8 × 10 −26 e cm and −3.8 < d n < 3.4 × 10 −26 e cm respectively. Taking limits that are symmetric about zero yields jd n j < 3.0 × 10 −26 e cm (90% C.L.) and jd n j < 3.6 × 10 −26 e cm (95% C.L.).
This analysis was also repeated with a number of variations: for example, with the spectrum derived from the softened-Maxwell approximation, which yielded (with the same uncertainty) central values of 0.86 and 0.79 × 10 −26 e cm with the slope multiplier k 0 kept fixed and allowed to vary, respectively; and without the apparently anomalous Run 1900 (resulting in central values of −0.05 × 10 −26 e cm with the simulated spectrum, and þ0.43 × 10 −26 e cm with the softened-Maxwell spectrum). Furthermore, we have also looked carefully at the implications of the fields as defined by the maps produced by the 2010 scans at PSI. We do not consider it appropriate to apply the results directly to our data, not only because the field scans were carried out after a lapse of several years and following dismantling, shipping and reassembly of the shields, but also because the accuracy of measurements of small transverse field components, typically a fraction of a percent of the main vertical field strength, is limited by the mechanical precision of the mapper system itself. Nonetheless, we considered it useful to carry out the analysis with the quadrupole shifts δR Q fixed at the values defined by these more recent field measurements. The   1 We note in passing that reanalysis of the earlier data has revealed a sign error in the central value of Eq. (4) in [36], as a consequence of which its reported upper limit (90% C.L.) should have been 7.4 rather than 6.3 × 10 −26 e cm. However, this is of historical interest only, since it obviously takes no account of the systematic effects detailed in the current work resulting value of d n was 0.76 × 10 −26 e cm (again with the same uncertainty), although as might be expected χ 2 =ν for the global fit deteriorated somewhat to 1.48. The range of variation of all of these results is of the order of the uncertainty on our systematic errors, and all lie comfortably within our confidence-limit range.
Ultimately, and in conclusion, after a thorough reevaluation of the previously known systematic effects and the inclusion of various other newly established effects, the overall corrections to the 2006 result have been found-both individually and collectively-to be relatively small compared with the statistical uncertainty. As a direct consequence, our new limit represents a departure of only a few percent from the earlier limit of jd n j < 2.9 × 10 −26 e cm (90% C.L.) [1] despite the very real differences between the two analyses.
Recently, the apparatus employed in this experiment has enjoyed something of a renaissance. It has been very substantially upgraded, including in particular with the provision of an array of high-precision cesium magnetometers [37], for a further nEDM measurement [38] at the new UCN source at PSI [39]. Furthermore, an innovative new spin-echo based technique [25] allows both the measurement of the UCN spectrum within the apparatus and an accurate determination of magnetic-field gradients. It is currently running with a greater sensitivity than ever before, and, with potential systematic errors now being very tightly constrained, new results will be limited only by statistics for some time to come.