Methods, Analysis, and the Treatment of Systematic Errors for the Electron Electric Dipole Moment Search in Thorium Monoxide

We recently set a new limit on the electric dipole moment of the electron (eEDM) (J. Baron et al., ACME collaboration, Science 343 (2014), 269-272), which represented an order-of-magnitude improvement on the previous limit and placed more stringent constraints on many CP-violating extensions to the Standard Model. In this paper we discuss the measurement in detail. The experimental method and associated apparatus are described, together with the techniques used to isolate the eEDM signal. In particular, we detail the way experimental switches were used to suppress effects that can mimic the signal of interest. The methods used to search for systematic errors, and models explaining observed systematic errors, are also described. We briefly discuss possible improvements to the experiment.


Introduction
Symmetries play a vital role in physics and experimental tests of symmetries have revealed insights into physical theory. Perhaps the most famous early example is the experiment of Michelson and Morley [1], now understood as an early demonstration of Lorentz invariance. Similarly, observed violations of parity (P ) symmetry [2] and charge-parity (CP ) symmetry [3] have informed and motivated understanding of the weak and strong forces [4,5]. The recent discovery of the Higgs boson [6] is a confirmation of a predicted spontaneously broken gauge symmetry [7], and the LHC continues to probe physics at high energies, looking for evidence of physics beyond the Standard Model (SM). On a complementary front, precision measurements of charge-parity-time (CP T ) invariance and Lorentz invariance using low-energy techniques continue to test these fundamental symmetries [8][9][10][11][12][13][14][15][16].
Precision measurements in atomic and molecular systems are well suited to testing fundamental physics, and searches for EDMs of fundamental particles have been at the forefront of such tests [17,18]. Measurements of the EDMs of the electron, neutron [19] and atomic species such as mercury [20], are complementary tests of beyond-SM physics and of fundamental symmetries [21]. As discussed in section 2.1, an EDM of a fundamental particle can only exist if time-reversal (T ) symmetry is broken, which is equivalent to CP violation for CP T -invariant models [22]. For many theories, intrinsic CP violation is predicted to manifest as eEDMs at an experimentally accessible level [21,23,24]. Consequently, discovering an eEDM, or further constraining its value, can inform our understanding of particle physics at high energy and help to shed light on outstanding issues such as the baryon asymmetry problem [25,26]. The current best limit on the eEDM was reported by ACME in 2014 [27]: † d e ≤ 9.3 × 10 −29 e · cm (90% conf. level). (1) Many extensions to the SM predict eEDM values many orders of magnitude higher than the SM prediction of < 10 −38 e·cm [17,21,30], meaning measurement of an eEDM at current experimental sensitivity would be a signature of new physics. Supersymmetry is an example of an extension to the SM that predicts a large, potentially measurable eEDM. The current eEDM limit constrains the parameter space associated with supersymmetry such that it is often considered unnatural [31,32].
In most models, the eEDM arises as a radiative correction (Feynman loop diagram) due to CP -violating interactions with new particles. An example of such an interaction within generic supersymmetric theory is shown in figure 1. The CP violation is associated with the presence of non-trivial complex phases in the ¡ Figure 1: Example of a supersymmetric 1-loop contribution to the eEDM. The symbols e L and e R represent the left and right helicities of the electron,ẽ a selectron,γ a photino and γ a photon. This generic diagram illustrates how a CP -violating phase (represented by the + symbol) can be produced in a straightforward manner by SM extensions. Note that a detailed discussion of associated high-energy theory is beyond the scope of this paper. theory. For a given CP -violating phase φ CP , one can make a generic estimate of the mass scale Λ of new physics being probed, according to the following formula for an n-loop process [24]: where e is the electron charge, m e the electron mass, α is the fine structure constant and φ CP is a CPviolating phase. Assuming that sin(φ CP ) ∼ 1 [21], we find that our most recent result interrogates energy scales for one-loop processes of around 10 TeV. Similar analysis shows that our result was sensitive to twoloop effects at around the 1 TeV mass scale. While any such estimates are inherently model-dependent, we see that using an apparatus that fits in a room we have been able to probe fundamental physics at energy scales usually associated with the largest particle accelerators.

Theory
The eEDM, d e , is a vector quantity that is aligned along (or against) the axis of the electron's spin, s [17]. By convention, we write d e = 2d e s, such that a measurement of any Cartesian component of d e yields a value of ±d e . (Here and throughout, we set = 1.) For an electron moving non-relativistically, the eEDM interacts with an electric field E via the Hamiltonian Under time reversal, T , s reverses direction, but E is unchanged. Similarly, under space inversion, P , s is unchanged, but E reverses direction. Hence H EDM is odd under P and T .
To measure the eEDM, one looks for an energy shift due to the interaction in equation 3. Since 1964, every improvement in experimental sensitivity to d e has been obtained by measuring this shift for electrons bound in a neutral atom or molecule [33][34][35][36][37][38][39][40][41][42]. This might seem surprising at first glance, since Schiff's theorem states that there can be no net electric field acting on a non-relativistic point particle bound in a neutral system [43]. However, in 1958 Salpeter showed that, when relativistic effects are taken into account, a neutral species can experience an energy shift due to an eEDM when an external electric field E ext is applied [44]. In 1965 Sandards showed, strikingly, that the size of the resulting energy shifts can be much larger than d e E ext [45].
More detailed explanations of this relativistic enhancement can be found elsewhere, e.g. [17,46,47], but we summarise the basic principle here. Taking into account the relativistic length contraction of the eEDM for a moving electron, its interaction with a total electric field E (the sum of an external, applied field and an intra-atomic/molecular field) takes the form where β = v/c is the dimensionless velocity and γ is the Lorentz factor [46]. The first term in this expression is the non-relativistic EDM interaction, whose expectation value vanishes by Schiff's theorem. The second, relativistic term can result in a nonzero net interaction when the electron's velocity and the electric field are non-uniform in space (as, for example, when the electron travels near a charged nucleus in an atom or molecule), and when the atom or molecule is polarised by an external electric field. This interaction can be expressed in terms of an 'effective electric field', E eff , defined in analogy to equation 3 such that Detailed calculations show that this 'effective electric field' within an atom or molecule can be significantly larger in magnitude than the applied external field. The size of E eff is maximal for systems where a valence electron has significant wavefunction amplitude near a highly-charged nucleus. In such species with a nucleus of atomic number Z, E eff scales approximately as [17] E eff ∝ P E Z 3 R(Z), where P E ∈ [0, 1) is the degree of electric polarisation of the state and R is a relativistic factor that is roughly constant for Z α −1 , but grows quickly as Z approaches α −1 [17,[47][48][49]. For fully polarised systems with Z ≈ 90 (as with our molecule of choice, ThO), the effective electric field can reach values as large as E eff ≈ 100 GV/cm. In practice, the maximum polarisation attainable with atoms, even in the highest laboratory static electric fields (∼ 100 kV/cm), is P E ∼ 10 −3 . Nevertheless, this can lead to values of E eff nearly 1,000 times larger than the applied laboratory field (e.g. E eff ∼ 70 MV/cm in Tl atoms [41]). Using this kind of enhancement, the limit on d e was reduced by six orders of magnitude by the first atom-based eEDM measurement [33]. Further improvement is afforded by working with polar molecules, which are much more polarisable than atoms due to having much more closely spaced levels of opposite parity (associated with their rotational motion). In practice, polarisation P E ∼ 1 is achievable with molecules [47][48][49]. In ThO, the effective electric field is E eff ≈ 78 GV/cm [28,29].
To measure the eEDM, the electron spin is prepared in a state oriented perpendicular to E eff , i.e. in a superposition of states parallel and antiparallel to E eff . After an interaction time τ , the eEDM energy shift in equation 3 produces a relative phase accumulation 2φ EDM = −2d e E eff τ between these states; this is equivalent to a precession of the spin orientation about E eff by an angle 2φ EDM .
For a shot-noise-limited measurement, the uncertainty in the eEDM, δd e , is given by where N is the number of measurements. The large values of E eff accessible in many molecules have motivated several recent eEDM searches [42,50,51]. This and other advantages associated with the molecule ThO are discussed in the following section.

ThO Molecule
ThO has a number of properties that make it well-suited to an eEDM measurement, both by enhancing statistical sensitivity and by suppressing systematic errors. We performed our measurements in the H electronic state of ThO, which has two valence electrons in a (σδ) 3 ∆ 1 state. Such states were first proposed for use in an eEDM measurement by Meyer et al. in 2006 [52]. The σ orbital valence electron wavefunction has a large amplitude near the heavy Th nucleus, facilitating the large E eff required for a large eEDM sensitivity, as described in section 2.1. The H state of ThO has one of the largest calculated values of E eff ≈ 77.6 GV/cm [28,29]. We note that the value of E eff in our experiment with ThO is more than 5 times larger than that attained in experiments using YbF, which set the previous eEDM limit [53][54][55], and over 1,000 times larger than that in experiments using Tl atoms [41]. All 3 ∆ 1 states have very small magnetic moments [56] since the δ 3/2 orbital valence electron serves to nearly cancel the magnetic moment of the σ 1/2 orbital. The actual magnetic moment of H deviates from zero primarily because of mixing with other states [57]. We express ThO molecule states using the basis |Y, J, M, Ω , where Y is the electronic state, J is the total angular momentum, M is the projection of J onto the laboratoryẑ-axis, and Ω is the projection of the electronic angular momentum onto the internuclear axis, n, which points from the lighter nucleus to the heavier nucleus. We used the |H, J = 1, Ω = 1 rotational manifold for our measurement, for which the magnetic moment is µ 1 = g 1 µ B M , where g 1 = −0.00440(5) is the associated g-factor [57,58] and µ B is the Bohr magneton. This small g-factor, generic to all molecules with this structure, ensures that the H state is particularly insensitive to spurious magnetic fields.
States with nonzero Ω have closely spaced pairs of opposite-parity levels with identical values of J called 'Ω-doublets', which are split by energy ∆ Ω due to the Coriolis effect in the rotating molecule [59][60][61]. The application of an electric field E mixes the M = 0 opposite-parity levels via the Stark interaction, − D · E, where D is the electric dipole operator, and from here on E is the applied (laboratory) field. In the limit | D · E| ∆ Ω , the molecule is fully polarised, the internuclear axis is nearly aligned or anti-aligned with the applied electric field, and the alignment orientation is described by quantum numberÑ ≡Ê ·n = ±1. This structure is shown for the H state of ThO in figure 2.
The use of molecules with Ω-doublet structure for an eEDM measurement, first explored in [48,63] in the context of experiments using PbO, is of great importance to our experiment. The |H, J = 1 manifold has an Ω-doublet splitting ∆ Ω,1 ≈ 2π × 360 kHz † [64] and an electric dipole moment D 1 ≡ | H, J = 1, M = ±1, Ω| D ·ẑ|H, J = 1, M = ±1, Ω | ≈ 2π × 1 MHz/(V/cm) [65]; this permits full (> 99 %) polarisation of the state in small applied electric fields, E 10 V/cm, allowing us to take full advantage of 10 V/cm, the associated Stark shift is linear and given by −Ñ D 1 |E|, where D 1 ≈ 2π × 1 MHz/(V/cm) (black arrow/lines) is the expectation value of the molecular electric dipole moment in these states [58]. Additionally, a magnetic field B causes a Zeeman shift ≈ − M g 1 µ B B z , with g 1 µ B ≈ −2π × 6 kHz/G (red arrows/lines) [57,62]. A nonzero eEDM would result in an additional energy shift ≈ −MÑ d e E eff (blue arrows/lines) whereẼ = −1 (+1) when the applied E field is (is not) reversed. The orientation of E eff (green arrows), the spin of the electron in the σ orbital (black arrow next to molecule), the external electric field E, and the external magnetic field B are shown relative to the laboratoryẑ direction which is oriented upwards on the page. Diagram not to scale.
the huge E eff in ThO. The Ω-doublet structure is also useful in rejecting systematic errors since it allows for spectroscopic reversal of E eff ∝ −n by addressing differentÑ states without reversing the applied electric field [66]. This is discussed in greater detail in section 5. 4.
The H state in ThO is metastable with a lifetime ≈1.8 ms [67], limiting our measurement time to τ ≈ 1 ms. We note that this is comparable to previous beam-based eEDM measurements where the atomic/molecular states used had significantly longer lifetimes [42,67,68].
As with many other species, ThO proved nicely compatible with a new approach to creating molecular beams, the hydrodynamically enhanced cryogenic buffer gas beam [69][70][71]. This method provides a cold, high-flux and low-divergence beam [72] yielding a large number of molecules in the few lowest-lying quantum states. The molecule beam's forward velocity (≈180 m/s) was also lower than a typical supersonic beam, which helped minimise the apparatus length for a given coherence time. For more details on the beam source, see section 3.2.2.

Basic Measurement Scheme
We performed a spin precession measurement, resembling previous beam-based eEDM experiments [40][41][42], on 232 Th 16 O molecules in a pulsed molecular beam generated by a cryogenic buffer gas beam source. Figure 4 shows a simplified schematic of the measurement. The molecules fly at velocity v ≈ 200 m/s into a magnetically shielded region with nominally uniform and parallel electric E and magnetic B fields. Molecule population is transfered from |X 1 Σ + , J = 1, M = ±1 in the electronic ground state to the metastable |H, J = 1, M = ±1, Ω =ÑẼM ≡ |±,Ñ state manifold (in the |±,Ñ nomenclature we use ± to refer to M = ±1) by optical pumping through the short-lived |A 3 Π 0 + , J = 0, M = 0 state with a 943 nm laser. This results in an even distribution of population in an incoherent mixture of the four |±,Ñ states in H. † Figure 3 shows the electronic states of ThO relevant to the eEDM measurement.  [67,73,74]. Solid arrows indicate transitions we address with lasers, wavy arrows indicate spontaneous decays of interest. For more details on how these transitions were used, see the main text.
In the absence of any experimental imperfections, we describe our system in terms of coordinate axes +ẑ along + E (for a specified sign of applied field that we denote as positive, pointing approximately east to west in the lab) and +x along the direction of the molecular beam (which travels approximately south to north) such that +ŷ is approximately aligned with gravity (cf. figure 4). Note that when we reverse the direction of the electric field, by construction the laboratory coordinate system does not change and the orientation of the electric field can be described byẼ ≡ sgn(ẑ · E) = ±1. Analogously, we reverse the direction of the magnetic field between twoB ≡ sgn(ẑ · B) = ±1 states. Since the directions of the fields are encoded byẼ andB, we define the magnitudes of the fields simply as B z ≡ |B z | and E ≡ | E|.
A superposition of the M = ±1 sublevels is prepared by optically pumping on the transition at 1090 nm between states |±,Ñ and |C 1 Π 1 , J = 1, M = 0 (|Ω = +1 −P|Ω = −1 )/ √ 2 ≡ |C,P , whereP = ±1 is the excited state parity ‡, with laser light linearly polarised in the xy plane. The resulting state corresponds to having the total angular momentum of the molecule aligned in the xy plane. Because the σ electron's spin is aligned with J, by the Wigner-Eckart theorem this is equivalent to aligning the spin [75], and we use this shorthand from here on. The state preparation laser frequency is tuned to spectroscopically select the molecule alignmentÑ , while the nearly degenerate M = ±1 states remain unresolved. The excited state C, which decays at a rate γ C ≈ 2π × 0.3 MHz, decays primarily (≈75 % [65]) to the ground state so that one superposition of the two |±,Ñ states is optically pumped out of H and the remaining orthogonal superposition, which is 'dark' to the preparation laser beam, is the prepared state. The linear polarisation † A glossary of symbols used throughout this paper is provided in section Appendix B. ‡ In this paper we follow the convention given in [60].
of the state preparation laser beam,ˆ prep , sets the relative coupling of each of the two |±,Ñ states to |C,P and determines the spin alignment angle of the remaining state in the laboratory frame. The bright superposition |B(ˆ prep ) is pumped away, and the orthogonal dark superposition |D(ˆ prep ) remains.
For the moment, we consider the specific caseP = +1 andˆ prep =x, (the general case will be discussed in section 3.1.2). In this case, the prepared state has the electron spin aligned along theŷ axis. As the molecules traverse the spin precession region of length L = 22 cm (which takes a time τ ≈ 1 ms), the electric and magnetic fields exert torques on the electric and magnetic dipole moments, causing the spin to precess in the xy plane by angle 2φ; this corresponds to the state where φ is given approximately by the sum of the Zeeman and eEDM contributions to the spin precession angles, φ = −(Bg 1 µ B B z +ÑẼd e E eff )τ.
The sign of the eEDM term,ÑẼ, arises from the relative orientation between E eff and the electron spin as illustrated in figure 2.
At the end of the spin precession region, we measure φ by optically pumping on the same H → C transition with the linearly polarised state readout laser beam. The polarisation alternates rapidly between two orthogonal linear polarisationsX andŶ , such that each molecule is subject to excitation by both polarisations as it flies through the detection region, and we record the modulated fluorescence signals F X and F Y from the decay of C to the ground state at 690 nm. This procedure amounts to a projective measurement of the spin ontoX andŶ , which are defined such thatX is at an angle θ with respect tox in the xy plane. To determine φ we compute the asymmetry, We set B z and θ such that φ − θ ≈ (π/4)(2n + 1) for integer n, so that the asymmetry is linearly proportional to small changes in φ and maximally sensitive to the eEDM. A simplified schematic of the experimental procedure just described is shown in figure 4.
By repeating the measurement of φ after having reversed any one of the signsÑ ,Ẽ orB, we may isolate the eEDM phase from the Zeeman phase. In practice, we repeat the phase measurement under all 2 3 (Ñ ,Ẽ,B) experiment states to reduce the sensitivity of the eEDM measurement to other spurious phases, and we extract the phase φ N E = −d e E eff τ = φ EDM . Here, we have introduced the notation φ u , discussed in detail in the next section, which we use throughout this document to refer to the component of φ that is odd under the set of switches listed in the superscript u, and implicitly even under those which are not listed (see section 3.1.2 and equation 23 for a rigorous definition). A component which is even under all switches is considered to be 'non-reversing' and is given an 'nr' superscript.

Measurement Scheme in Detail
To fully describe the method by which we extracted d e from the data in section 4, and to describe the systematic error models in section 5, we must introduce some additional formalism to describe the spin precession measurement to generalize the simple case described in the previous section.
We work in the regime in which the Stark shift in H is approximately linear, E Stark ≈ −Ñ D 1 E, which holds when the Stark interaction energy is large compared to the Ω-doublet energy splitting ∆ Ω,1 but small compared to the rotational energy scale, described by the H-state rotational constant B H ≈ 2π× 9.8 GHz, i.e. ∆ Ω,1 D 1 E B H . In this regime, the molecular alignment is approximately related to Ω byÑ =ẼM Ω; this relation is assumed throughout this document. This is a good approximation, but it is notable that due to the Stark interaction at first order in perturbation theory, each |M,Ñ state is a superposition of all four |H, J, M, Ω states with J = 1, 2 and Ω = ±1. This effect is discussed further in sections 5.2.6 and 5.6.2.

State
J alignment:Ñ 2φ 2φ/τ Figure 4: Simplified schematic of the measurement scheme; numbers next to energy levels label J. 1.
Molecules in the |X, J = 1 state are optically pumped via the A state into |H, J = 1 by a retroflected (and offset in x) laser beam (blue arrows into/out of page), polarised alongx andŷ (blue arrows). 2.
Molecules from one of theÑ states are then prepared in a superposition of M sublevels (M = −1, 0, +1 from left to right) by a linearly polarised laser beam (red) addressing the H → C transition. This aligns the molecule's angular momentum, J, which in turn aligns the spin of the eEDM-sensitive σ electron, which is on average aligned with J. 3. The angular momentum (and hence electron spin) then precesses due to the electric and magnetic fields present (into the page) by an angle φ. This precession is dominated by the magnetic interaction but also includes a term linear in d e (see equation 10). 4. The spin state is projected onto orthogonal superpositions of the M sublevels by laser beams polarised alongX,Ŷ (red arrows). The resulting fluorescence is determined by the population in each superposition state and hence the precession angle φ.
Let us consider the preparation of a spin-aligned state again. Starting from an incoherent mixture of the four |±,Ñ states, we perform optical pumping on the electric dipole transition between |±,Ñ and |C,P , for a specificÑ , with laser light of polarisationˆ prep that is nominally linear in the xy plane. This step depletes the bright superposition state (see e.g. [76]) whereˆ ±1 = ∓ (x ± iŷ) / √ 2 are unit vectors for circular polarisation. The corresponding dark state (with which the laser does not interact) is the orthogonal superposition This dark state serves as the initial state, |ψ(0),Ñ = |D(ˆ prep ,Ñ ,P = +1) , for the spin-precession experiment, where we fixed the state preparation laser frequency to address the excited state with paritỹ P = +1. The state preparation laser polarisation can be parameterised aŝ where Θ prep ≈ π/4 defines the ellipticity Stokes parameter (S 3 /I) prep = cos 2Θ prep ≈ 0, and θ prep defines the linear polarisation angle with respect tox in the xy plane. From here on, we refer to the ellipticity Stokes parameter as S ≡ S 3 /I. There is a one-to-one correspondence between the dark state superposition and the projection of the laser polarisationˆ prep onto the xy plane. If the laser polarisation does not lie entirely in the xy plane, equations 12 and 13 are still appropriate, but require normalization. Note that if the laser is linearly polarised, switching the excited state parityP has the same effect on the dark state as rotating the laser polarisation angle by π/2.
Following the initial state preparation, the molecules traverse the spin-precession region with their forward velocity nominally alongx. In this region there are nominally uniform and parallel electric ( E) and magnetic ( B) fields, which produce energy shifts given by where D 1 is the electric dipole moment of |H, J = 1 . Here η = 0.79(1) nm/V accounts for the E-dependent magnetic moment difference between the two sets ofÑ levels in |H, J = 1 [57], as described in section 4.2.2.
The energy shift terms that depend on the sign of M contribute to the spin precession angle φ, which is given by: This phase is dominated by the magnetic (Zeeman) interaction. The Stark shift, proportional to |M |, does not contribute. The state then evolves to: (recall |ψ(0),Ñ = |D(ˆ prep ,Ñ ,P = +1) per equation 13) and molecules enter a detection region where the state is read out by optically pumping again between the |H, J = 1 and |C, J = 1 manifolds. This optical pumping is performed alternately by two laser beams with nominally orthogonal linear polarisationsˆ X and Y . † These beams excite the projection of |ψ(τ ),Ñ onto the bright states (with the sameÑ that was addressed in the state preparation optical pumping step, but with an independent choice ofP) with probability P X,Y respectively. In the ideal case in which all laser polarisations are exactly linear, this probability is given by where θ X,Y are the linear polarisation angles of the state readout beams, with respect tox. The result is a signal that varies sinusoidally with the precession angle φ. To measure these probabilities, we observe the associated modulated fluorescence signals, F X,Y = f N 0 P X,Y , where N 0 is the number of molecules in the addressedÑ level at the state readout region, and f is the fraction of total fluorescence photons that are detected.
To distinguish between molecule number fluctuations and phase variations, we normalize with respect to the former by rapidly switching the state readout laser between the two orthogonal polarisations,ˆ X,Y , every 5 µs. This is significantly quicker than fluctuations in the molecule number and is sufficiently quick that every molecule is interrogated by both polarisations (see section 4 or [62] for more details). We then form an asymmetry A, which is immune to molecule number fluctuations, given by where we have assumed that the readout polarisations are exactly orthogonal, given by θ X = θ read and θ Y = θ read + π/2, and where we have defined θ ≡ θ read − θ prep . ‡ In this equation and from now on unless otherwise noted,P refers to the excited state parity that is addressed by the state readout laser, not to be confused with the excited state parity addressed by the state preparation laser, which is kept fixed. The value of B z and the state preparation and readout laser beam polarisations are chosen so that |φ − θ| ≈ π/4. This corresponds to the linear part of the asymmetry fringe in equation (20), where A is most sensitive to, and linearly proportional to, small changes in φ (cf. figure 23). A variety of effects including imperfect optical pumping, decay from C back to H, elliptical laser polarisation and forward velocity dispersion, reduce the measurement sensitivity by a 'contrast' factor with |C| ≤ 1. We measure this parameter by dithering θ = θ nr + ∆θθ (where θ nr is the average or 'nonreversing' polarisation angle) between states ofθ = ±1, with amplitude ∆θ = 0.05 rad. We found that typically |C| ≈ 0.94. We then extract the measured phase, Φ = A/(2C)+qπ/4, by normalising the asymmetry measurements according to the measured contrast -see section 4 for more details on the data analysis methods used to evaluate this quantity. In the ideal case, the measured phase matches closely with the precession phase, Φ ≈ φ. However, a variety effects that are investigated closely in section 5 lead to slight deviations between these two quantities, which can contribute to systematic errors in the measurement. Unless explicitly specified, C is assumed to be an unsigned quantity from here on. In particular, when averaging over multiple states of the experiment, |C| is used.
To isolate the eEDM term from other components of the energy shift in equation (15), the experiment is repeated under different conditions that are characterised by parameters whose sign is switched regularly during the experiment. The spin precession measurement is repeated for all 2 4 experiment states defined by the four primary binary switch parameters:Ñ , the molecular orientation relative to the applied electric field (changed every 0.5 s);Ẽ, the direction of the applied electric field in the laboratory (2 s);θ, the sign of the readout polarisation dither (10 s); andB, the direction of the applied magnetic field in the laboratory (40 s). For each (Ñ ,Ẽ,B) state, the asymmetry A(Ñ ,Ẽ,B), contrast C(Ñ ,Ẽ,B), and measured phase Φ(Ñ ,Ẽ,B) are determined as described earlier. The data taken under all 2 4 = 16 experimental states derived from these four binary switches constitutes a 'block' of data.
We can write the phase Φ(Ñ ,Ẽ,B) in terms of components with particular parity with respect to the experimental switches: We refer to these components as 'switch-parity channels'. A channel is said to be odd with respect to some subset of switches (labelled as superscripts) if it changes sign when any of those switches is performed. Thus it will also change sign if an odd number of those switches is performed. It is implicitly even under all other switches. We use this general notation throughout this document to refer to correlations of various measured quantities and experimental parameters with experiment switches. To generalize, if we have k binary experiment switches (S 1 ,S 2 , . . . ,S k ) such thatS i = ±1, and we perform a measurement of the parameter X(S 1 ,S 2 , . . . ,S k ) for a complete set of the 2 k switch states, then the component of X that is odd under the product of switches S aSb . . . is given by The switch parity behavior of a given component is expressed in the superscript which lists the experimental switches with respect to which the component is odd. We order the switch labels in the superscripts such that the fastest switches are listed first and the slowest switches are listed last. Some components give particularly important physical quantities. Most notably, the eEDM precession phase is extracted from thẽ NẼ-correlated component of the measured phase: that is, in the ideal case Φ N E = −d e E eff τ . Additionally, the Zeeman precession phase is nominally given by Φ B = −µ B g 1 B z τ . Recall we label 'non-reversing' components with an 'nr' superscript. In a few cases, we drop the superscript parity because it is redundant. For example, we drop the superscript on the dominant components of the applied electric and magnetic fields, E ≡ E E and B z ≡ B B z . Many other experimental parameters are also varied between blocks of data to suppress and monitor systematic errors (figure 5). These 'superblock' switches include: excited-state parity addressed by the state readout laser beams,P (chosen randomly after every block, with equal numbers ofP = ±1); simultaneous change of the power supply polarity and interchange of leads connecting the electric field plates to their voltage supply,L (4 blocks); a rotation of the state readout polarisation basis by θ read → θ read + π/2 to interchange the roles of the X and Y beams,R (8 blocks); and a global polarisation rotation of both state preparation and readout lasers by θ read → θ read + π/2 and θ prep → θ prep + π/2,G (16 blocks).
Additionally, the magnitude of the magnetic field, B z , was switched on the timescale of 64-128 blocks (∼1 hour), and the magnitude of the applied electric field, E, and the laser propagation direction,k ·ẑ, were changed on timescales of ∼1 day and ∼1 week, respectively.
On these longer timescales, we also alternated between taking eEDM data under Normal conditions, for which all experiment parameters were set to their nominally ideal values, and taking data with Intentional Parameter Variations (IPVs), during which some experimental parameter was set to deviate from ideal so that we could monitor the size of the known systematic errors described in section 5.2.6. We took IPV data in which we varied (a) the non-reversing electric field E nr and (b) theÑẼ-correlated Rabi frequency, Ω N E r , to measure the sensitivity of the eEDM measurement to these parameters and we varied (c) the state preparation laser detuning ∆ prep to monitor the size of the residual E nr . These systematic errors are discussed in more detail in sections 3.2.5 and 5.2.6.
The details of the data analysis required to extract the eEDM-correlated phase Φ N E are described in section 4. A lower bound on the statistical uncertainty δΦ N E of the eEDM-correlated phase is given by photoelectron shot noise to be δΦ N E = 1/(2|C| √ N ) for N detected photoelectrons [17,77]. In the case where shot noise is the sole contribution, we can express the statistical uncertainty δd e in our measurement of the eEDM as whereṄ ≈ fṄ 0 is the measurement rate (equivalent to the photoelectron detection rate) and T is the integration time (recall f is the fraction of fluorescence photons detected and N 0 is the number of molecules in the addressedÑ level). Further discussion of the achieved statistical uncertainty is presented in section 4.

Overview
In this section we provide an overview of our experimental procedure and the important components of our apparatus. The reader should consult subsequent subsections for further details. A schematic of the experimental apparatus is shown in figure 6. ThO molecules were produced via pulsed laser ablation of a ThO 2 ceramic target. This took place in a cryogenic neon buffer gas cell, held at a temperature of ≈16 K, at a repetition rate of 50 Hz. The resulting molecular beam was collimated and had a forward velocity v ≈ 200 m/s. In the state readout region the molecular pulses had a temporal (spatial) length of around 2 ms (40 cm). The buffer gas beam source is described in detail in section 3.2.2.
After leaving the buffer gas source, the molecules had a velocity distribution and rotational level populations consistent with a Maxwell-Boltzmann distribution at a temperature of ≈4 K. This was lower than the cell temperature due to expansion cooling, which enhanced the number of usable ThO molecules in the relevant rotational state. Further rotational cooling was provided via optical pumping and microwave mixing (see section 3.2.3). The molecules then passed through adjustable horizontal and vertical collimators consisting of a double layer of razor blades affixed to linear translation vacuum feedthroughs. Under normal running conditions, these collimators were withdrawn so that they did not affect the profile of the molecule beam in the spin-precession region; however, they were used to modify the spatial profile of the molecule beam during systematic checks to investigate the effect of molecule beam position and pointing. Just before the field plates, 126 cm from the beam source, the molecules passed through a 1 cm square collimating aperture, which determined the beam profile in the spin-precession region and prevented particles in the beam from being deposited on the field plates.
As described in section 3.1, a spin precession measurement was performed where the precession angle provided a measure of the interaction energy of an eEDM with the effective electric field, E eff , in the molecule. A pair of transparent, ITO-coated glass plates provided an electric field that polarised and aligned the molecules. Laser beams passed through these plates to perform state preparation and readout. Around the vacuum chamber were coils that provided a uniform magnetic field in the +ẑ direction, and five layers of ect of a theory errorbar.
ect of a theory errorbar.  magnetic shielding which shielded against environmental magnetic fields. The electric and magnetic fields are discussed in detail in sections 3.2.5 and 3.2.6. The fluorescence induced by the state readout laser beam was collected by a set of eight lenses and transferred out of the spin-precession region using fiber bundles and light pipes (see section 3.2.7), where it was detected by photo-multiplier tubes †.

Buffer Gas Beam Source
The basic operation of our beam source [64,69,71,72,[78][79][80][81][82][83][84][85][86][87][88][89] is depicted in figure 7. Neon buffer gas was flowed at a rate of ≈ 30 SCCM (standard cubic centimetres per minute) through a copper cell held at a T ≈ 16 K. The inside of the cell was cylindrical with a diameter of 13 mm and a length of 75 mm. Within the cell ThO was introduced at high temperature via laser ablation: overlapped beams of light with wavelengths † Hamamatsu R8900U-20.  Figure 6: A schematic of the overall ACME experimental apparatus. A beam of ThO molecules was produced by a cryogenic buffer-gas-cooled source. After exiting the source, the molecules were rotationally cooled via optical pumping and microwave mixing and then collimated before entering a magnetically shielded spin-precession region where nominally uniform magnetic and electric fields were applied. Using optical pumping, the molecules were transferred into the eEDM-sensitive H state and then a spin superposition state was prepared. The spin precessed for a distance of ≈22 cm and was then read out via laser-induced fluorescence. The fluorescence photons were collected by lenses and passed out of the chamber for detection by photomultiplier tubes. See main text for further details. 532 nm and 1064 nm emitted by a pulsed Nd:YAG laser † were focussed onto a 1.9 cm diameter ThO 2 target fabricated from pressed and sintered powder [90,91]. The laser pulses had a duration of a few ns, a pulse energy up to approximately 100 mJ and a repetition rate of 50 Hz. The resulting hot plume of ejected particles, which contained ThO along with various other ablation byproducts, was cooled by collisions with the neon buffer gas, became entrained, and then exited the cell. The cell temperature was maintained by a combination of a pulse tube refrigerator † and a resistive heater. The cell was surrounded by a 4 K copper shield that protected the cell from black-body radiation and cryopumped most of the neon emerging from the cell. This shield was also partially covered with activated charcoal that acted as a cryopump for residual helium in the neon buffer gas. We observed a background pressure of 10 −7 Torr without any mechanical pumping of the beam source when cold and with no buffer gas flow. The 4 K shield had a stainless steel conical collimator with a circular aperture of diameter 6 mm, located 25 mm from the cell aperture, by which distance the expanding beam was sufficiently diffuse that intra-beam collisions were negligible and most trajectories were ballistic. This collimator thus functioned as a differential pumping aperture without affecting the beam's cooling, acceleration or expansion [72]. The collimator had a thermal standoff relative to the 4 K shield to which it was mounted so that it could be kept at a temperature above the freezing point of neon by a resistive heater to prevent ice buildup on the collimator adversely affecting the beam dynamics. Another layer of shielding surrounded the 4 K copper shield, constructed from aluminium and held at a temperature of 60 K. Both the 4 K and 60 K radiation shields were thermally connected to the pulse tube by heat links made of flexible copper rope.
The aluminium vacuum chamber that housed the buffer gas beam source ‡ had windows on each side, providing optical access for both the ablation laser and spectroscopy lasers, the latter allowing characterisation and monitoring of beam properties. The ThO beam's forward velocity distribution was roughly Gaussian with mean v ≈ 200 m/s and standard deviation σ v ≈ 13 m/s, corresponding to a temperature of ≈5 K. The rotational temperature was T rot ≈ 4 K (rotational constant B X ≈ 0.33 cm −1 ), meaning that ≈90% of the population was contained in the levels J = 0-3. Upon exiting the cell, the beam had a FWHM angular spread of ≈ 45 • . Several stages of collimation were applied before reaching the spin-precession region. The final collimator subtended a solid angle of ≈ 6 × 10 −5 sr, meaning 1 in ∼20, 000 molecules exiting the cell reached the spin-precession region, where the precession measurement was performed (see figure 6).
ThO yields from a given ablation spot decreased significantly after ∼10 4 − −10 5 YAG pulses (∼10 mins), at which time the laser spot was moved to an un-depleted region via a motorised mirror to re-optimise the beam flux. Each target was found to provide acceptable levels of molecule flux for around 300 hours of continuous running (≈5 × 10 7 shots) before requiring replacement.

Rotational Cooling
We observed that ≈2 cm downstream of (further from) the buffer gas beam source cell aperture, J-changing collisions were 'frozen out' [72], and the distribution of rotational state populations was fairly well described by a Boltzmann distribution with temperature T rot ≈ 4 K. At this temperature the resulting fractions of molecules in the J = 0-3 levels were estimated to be 0.1, 0.3, 0.3 and 0.2 respectively.
As described in section 3.2.4, we sought to transfer as much of the initial ground state population as possible into |H, J = 1 via optical pumping. To enhance the population which was transferred, we accumulated population in a single rotational level of the ground state before state preparation. The scheme used to achieve this, which we refer to as rotational cooling, is illustrated schematically in figure 8 and discussed in detail in [92]. The first stage of the process was the optical pumping of molecules out of |X, J = 2 (|X, J = 3 ), via |C, J = 1 (|C, J = 2 ) into |X, J = 0 (|X, J = 1 ) using laser light at 690 nm. The natural linewidth of the X → C transition is ≈2π × 0.3 MHz, however the usable molecules had a ≈ ± 0.7 m/s transverse velocity spread, corresponding to a 1σ Doppler width of ≈2π × 1.5 MHz at 690 nm. Because the lasers used had linewidths of 1 MHz, to completely optically pump these molecules we relied on a combination of power broadening and extended interaction time. Optical pumping occured in a magnetically unshielded region where a background field B ≈ 500 mG was present; however, the magnetic moment of X (C) is ∼µ N (≈µ B /J(J + 1)), the nuclear magneton, which led to a Zeeman shift of ∼2π × 400 Hz ( 2π × 400 kHz) such that the M sublevels were not resolved by our lasers. The |C, J = 1 state has an Ω-doublet splitting of ∆ Ω,C,J=1 ≈ 2π × 51 MHz [93]. This splitting scales as ∆ Ω,C,J ∝ J(J + 1), meaning we could spectroscopically resolve the Ω-doublets for all |C, J . In addition, having no E-field present meant that the M sublevels of C and X remained unresolved and the energy eigenstates remained parity eigenstates. The X state is also insensitive to E-fields due to the lack of Ω-doublet substructure; opposite parity states are separated by † Cryomech PT415. ‡ Precision Cryogenic Systems Inc.  Figure 8: Schematic of the rotational cooling process. Numbers label J and M y (projection of total angular momentum along y) sublevels are unlabelled but are −1, 0, +1 from left to right. Population was first optically pumped out of the J = 2 and J = 3 levels (C-state Ω-doublet structure and M y sublevels omitted for clarity) in a nominally field-free region. Next, population was equilibrated between |J = 0 and |J = 1, M y = 0 via microwave pumping. An electric field of ≈40 V/cm alongŷ was empirically observed to lead to an increased population in |X, J = 1 . Grey dots represent population before these pumping processes. The schematic on the right represents the populations inside the spin-precession region (after pumping).
∼10 GHz and were hence unmixed. Laser beams with linear polarisation alternating betweenx andŷ were used to ensure that all population in |X, J = 2, 3 was addressed. This was achieved by directing around 10 passes of the beam, offset in x, through the vacuum chamber, passing through a quarter-wave plate twice in each pass, over a distance of around 2 cm. The laser light for rotational cooling was derived from home-built extended cavity diode lasers (ECDLs). The lasers were frequency-stabilised using a scanning transfer cavity with a computer-controlled servo [94]. Frequency-doubled light at 1064 nm from a frequency-stabilised Nd:YAG laser, locked to a molecular iodine line via modulation transfer spectroscopy [95], provided the reference for the transfer cavity.
After this first stage of rotational cooling, there was significantly greater population in the |X, J = 0 state than in any of the |X, J = 1, M sublevels. We obtained a ≈25 % increase in the J = 1 population by applying a continuous microwave field, resonant with the J = 0 → J = 1 transition; a sufficiently high microwave power combined with the inherent velocity dispersion of the molecule beam led to an equilibration of population between the coupled levels [92]. In this second stage of rotational cooling it was empirically observed that applying an electric field to lift the M y sublevel degeneracy was necessary to obtain the increased population in |X, J = 1 . A pair of copper electric field plates (spacing ≈ 4 cm) provided a field of ≈40 V/cm in theŷ (vertical) direction. We applied microwaves resonant with the Stark-shifted |J = 0 → |J = 1, M y = 0 transition at a frequency of 2π × 19.904521 GHz from an ex vacuo horn. Between the rotational cooling and spin-precession regions of the experiment (see figure 6) there was not a well-defined quantisation axis, and we observe that the populations of the |J = 1, M magnetic sublevels were equalised by the time the molecules reached the state preparation region.
Overall, we find that rotational cooling provided a factor of between 1.5 and 2.0 increase in the molecule fluorescence signal F in the state readout region. This gain factor was observed to vary slowly over time, possibly due to variations in the rotational temperature of the molecule beam, with significant changes sometimes observed when the ablation target was changed.

State Preparation and Readout
Following rotational cooling, the molecular beam passed into the spin-precession region, where the molecules experienced a nominally uniform electric field, E, which was nominally collinear with a magnetic field, B. Note that since neither of the states X 1 Σ + nor A 3 Π 0+ have Ω-doublet structure, parity remained a good quantum number for these levels for the small (∼100 V/cm) electric fields we applied.
We transferred the molecules into the H electronic state via optical pumping, as illustrated in figure 9. A 943 nm laser beam nominally propagating alongẑ excited molecules from the |X, J = 1 to |A, J = 0 . The laser beam passed through a quarter-wave plate, was retroflected and offset in x, then passed again  through the quarter-wave plate, such that the molecules were pumped by two spatially separated laser beams of orthogonal polarisations, allowing all population in both the |X, J = 1, M = ±1 levels to be excited. After excitation to A, the molecules could spontaneously decay into the |H, J = 1 manifold of states. We observed a transfer efficiency from X to H of ≈0.3 [92]. In this decay, five out of the six sublevels were populated; 1/6 of the population decayed to each of |H, M = ±1,Ñ = ±1 and 1/3 to |H,P = −1, M = 0 (see sections 2.2 and 3.1 for definitions ofÑ andP); decay to |H,P = +1, M = 0 is forbidden. Of these five populated states, only one corresponded to the desired initial state described by equation 13, and only 1/6 of the population in the H state was in this desired state. We estimated a total transfer efficiency from |X, J = 1, M = ±1 to the state in equation 13 of 30% × 1/6 = 5%.
The 943 nm laser light was derived from a commercial ECDL and then amplified by a commercial tapered amplifier †, generating ≈ 400 mW. As with the rotational cooling lasers, we verified that the power was sufficient to drive optical pumping to completion across the entire transverse velocity distribution of the molecular beam. This laser was also stabilised via the previously described (section 3.2.3) transfer cavity. The frequency of the laser light was monitored every 30-60 mins by scanning across the molecular resonances, allowing for independent fine-tuning and compensation of long-term frequency changes ( 2π × 100 kHz per half hour) due to e.g. temperature drifts in the cavity.
Around 1 cm downstream of the optical pumping laser beam that transferred population to H, we prepared the initial state of H (equation 13) by driving the transition between |H, M = ±1,Ñ and |C,P = +1 (see section 3.1 for more details) using laser light at 1090 nm. A distance L = 22 cm downstream of the preparation laser, a second 1090 nm laser beam was used to read out the molecule state via the same transition (but with the option to excite to eitherP state). This laser light was also derived from a commercial ECDL. It was then amplified using a fiber amplifier ‡, increasing the power to ≈250 mW. AOMs were then used to split and frequency shift the light to address bothÑ states in the H state, allowing spectroscopic selection of molecular alignment, and of bothP levels in the C state. Switching between these frequencies was achieved with either RF switches § or a DDS synthesizer . Given the linear Stark shifts D 1 E ≈ 2π × 146 MHz (2π × 37 MHz) in H with an applied electric field strength |E| = 141 V/cm (36 V/cm), and the excited state Ω-doublet splitting ∆ Ω,C,J=1 ≈ 50 MHz in C, these transitions were spectroscopically well-resolved. We fixed the nominal frequency of the state preparation laser to only addressP = +1, but periodically switched the state readout laser frequency to addressP = ±1 (∼1 min period). The transition frequencies of the state preparation and state readout laser beams were changed synchronously to always address the sameÑ level, with a switch betweenÑ levels every 0.5 s. The state preparation and readout laser beams were then independently amplified with a pair of fiber amplifiersP, providing ∼3-4 W of power. Immediately before † Toptica DL Pro and BoosTA. interrogating the molecules, the polarisation of the state readout laser beam was rapidly (100 kHz) switched between two orthogonal linear polarisations. The scheme for producing theÑ andP switches, and this fast polarisation switch, together with the corresponding laser transitions, is shown in figure 10. We now describe in detail how the appropriate frequency laser light was produced.   where ω N L is half the splitting between the twoÑ states; these AOMs were switched on and off to perform theÑ switch. The two frequency-shifted beams were combined and overlapped. For state preparation (lower branch of diagram), another AOM shifted the light by +ω L,1 , into resonance with the lower Ω-doublet in C (P = +1). This light was then amplifed again and passed through an AOM to vary the power (used as a systematic check). For the state readout (upper branch of diagram), a single AOM switched frequency to produce shifts +ω L,2,3 for the twoP states. A relative detuning between state preparation and readout laser beams (not shown) was also implemented with this AOM. (Shifts common to both beams were made by changing ω 0 .) The light was then amplified again and passed through an AOM to vary the power. Finally, polarisation switching was achieved with two AOMs switched on and off at 100 kHz, π out of phase with each other; light not diffracted (and frequency shifted by −ω L,PS ) by the first AOM was diffracted (and also frequency shifted by −ω L,PS ) by the second AOM. The diffracted light from each path was combined on a polarising beam splitter such that the linear polarisation of the final output beam alternated.
Based on the notation above we can now write the components of the frequencies of the state preparation and readout laser beams which do not reverse with any experimental switch as ω nr L,prep = ω L,0 + ω L,1 and ω nr L,read = ω L,0 + (ω L,2 + ω L,3 )/2 − ω L,PS , respectively. We can also write theP-correlated frequency component of the state readout laser as ω P L,read = (ω L,2 − ω L,3 )/2. We then write the detuning components as ∆ i = ω L,i − ω HC where i ∈ {prep, X, Y } indexes the laser and ω HC is the transition frequency between the line centres of the |H, J = 1 and |C, J = 1 manifolds †. We can rewrite this overall detuning in terms of various switch parity components: In the above equations we have defined detuning components of given switch parities -we shall now explain each component in turn.
is the mismatch between the Stark shift D 1 E(x i ) and the AOM frequency ω N L used to switch between resonantly addressing the twoÑ states, where x i is the x position of laser beam i. ∆ N E i = D 1 E nr (x i ) is a detuning component correlated like an eEDM signal which is due to a non-reversing component of the applied electric field. To understand this relation, consider figure 11. Recall that ∆ Ω,C,J=1 is the Ω-doublet splitting of the C state. For a E nr = 0, |E|, and hence the splitting between theÑ levels in H, depends onẼ. If the laser frequency for eachÑ is set assuming E nr = 0, a nonzero E nr leads to blue or red detuning from resonance, correlated withẼ. Because the sign of the Stark shift is correlated withÑ , the resulting detuning is also correlated withÑ . ∆ P i = ω P L + ∆ Ω,C,J=1 /2 is the mismatch between the excited state parity splitting and the AOM frequency, ω P L = (ω L,3 − ω L,2 )/2, used to switch between the two states (δ i,{X,Y } is the Kronecker delta, 1 if i = X or i = Y , zero else). We observed that ∆ N (∆ P ) was typically less than 2π × 20 kHz (2π × 50 kHz). Although we could measure ∆ N with ∼2π × 1 kHz precision, fluctuations in the Stark splitting, likely caused by thermally-induced fluctuations of the field plate spacing, limited our ability to zero out this correlated detuning.
We define ∆ nr = (∆ nr prep + (1/2)(∆ nr X + ∆ nr Y ))/2 as the average non-reversing detuning of the state preparation and readout laser beams; its value typically fluctuated by ∼2π × 0.1 MHz over several hours. Every 30-60 minutes the value of ∆ nr was scanned across the molecular resonance in the readout region using the ∆-tuning AOM (see figure 10), as an auxiliary optimisation. ∆ nr was set to the value where the fluorescence signal was maximum. This ensured that the average detuning of the state readout laser beams, (∆ nr X + ∆ nr Y )/2, was zero, however, if the state preparation and readout laser beams were not exactly parallel, there could be a difference between ∆ nr i due to the resulting difference in Doppler shifts. The effect of a detuning difference between the two state readout polarisations ∆ XY = (∆ nr X − ∆ nr Y )/2 is discussed in section 5.3. Additionally, each day we scanned the frequency of the preparation laser across the molecule resonance while monitoring the contrast of our fluorescence signal to ensure ∆ nr prep was kept below 2π×0.2 MHz (an example scan is shown in figure 24). The ways in which detuning components can contribute to systematic errors are discussed in detail in sections 5.2.3 and 5.2.6.
Other polarisation switches of the state preparation and readout laser beams (R andG) were controlled independently via half-wave plates mounted in high resolution rotation stages †. These switches and their use in the experiment are described in detail in section 4. Both beams were shaped using cylindrical lenses to be extended in y so all molecules in the beam were addressed. The Gaussian standard deviations of the beam intensities were 1.1 mm and 7.5 mm in the x and y directions, respectively [92]. The preparation laser beam was temporally modulated at 50 Hz with a chopper wheel, synchronous with the molecule beam pulses, to minimise the incident power on the field plates so as to reduce an important systematic error, described in sections 5.2.3 and 5.2.4.

Electric Field
The applied E-field was generated with a pair of 43 cm × 23 cm parallel conducting plates composed of ≈1.25 cm thick Borofloat glass, coated with a ∼200 nm layer of indium tin oxide on the inner faces ‡. The plates were transparent to the X → A optical pumping laser (943 nm), the H → C state preparation and readout lasers (1090 nm), and the C → X molecule fluorescence (690 nm). The outside faces of the electric field plates were prepared with a broadband anti-reflection coating with a specified <1% reflectivity at normal incidence from 600-1000 nm. The plates were made much larger than the precession region in order to minimise inhomogeneity of the field through which the molecules passed, and to enable large solid angle collection of fluorescence through the plates. One of the field plates was mounted in an aluminium frame fixed to the base of the vacuum chamber. The other field plate was secured a distance of 2.5 cm away in a kinematic aluminium frame. On the inward-facing surfaces, a frame of gold-plated copper clamped each field plate to the aluminium mounts and also functioned as a 'guard ring' electrode, suppressing the effect of fringing fields near the edges of the plate. The field plates were protected from impinging molecular beam particles by a 1 cm × 1 cm square collimator fixed to the entrance of the assembly.
The applied electric field was controlled by a 20-bit DAC, amplified to produce up to ±200 V §. The field plate assembly was referenced to the vacuum chamber ground. Equal and opposite voltages, ±V , were applied to each side of the assembly. The direction of the field (theẼ switch) was reversed every 1-2 s by † Newport URS50BCC. ‡ The plates were fabricated by Custom Scientific, Inc. § PA98A Power OpAmp.
reprogramming the output of the DAC channels to reverse their polarity. The configuration of the electrical connections between the amplified voltage and the field plates, denoted byL, was reversed via a pair of mercury-wetted relays every 2.6 minutes †. Data were also taken with two different values of E = 36 and 141 V/cm, varied on a ∼1 day time scale.
We measured the homogoneity of the electric field in a number of ways which we shall describe in turn now. Firstly, an indirect measure was obtained by determining the spatial variation of the field plate separation d using a 'white light' Michelson interferometer [96]. A schematic of the setup is shown in figure 12. We directed a light beam at normal incidence through the electric field plates. This resulted in d Light source Figure 12: Schematic of the apparatus used to perform an interferometric measurement of the electric field plate separation. A spectrally broad light beam is reflected perpendicularly off the field plates and passes into a conventional Michelson interferometer setup with one fixed arm (length L 2 ) and one movable arm (length L 1 ). An example of a pair of beam paths of interest is shown as solid and dashed red lines. If the two paths are slightly tilted relative to each other, a spatial interference pattern (inset) is observed on the CCD detector when the path length difference between the two beams is less than the coherence length, e.g.
multiple reflected beams, but we restrict discussion to the reflections from the conducting surfaces as these are of primary interest and were efficiently experimentally isolated from all others. The reflected beams passed into a Michelson interferometer with one arm of fixed length (L 2 ) and one with length adjustable via a micrometer (L 1 ). Constructive (destructive) interference occured whenever the lengths of two reflected beam paths differed by an integer (odd half-integer) multiple of the wavelength of the light. This condition was restricted further by the use of a broadband superluminescent diode ‡ with a short coherence length L c (nominally L c ≈ 15 µm). Thus the interference was only substantial when the two beams differed in length by L c . This occurred when L 1 = L 2 (for reflections off the same surface) or when L 1 = L 2 ± d (for reflections off surfaces spaced by d). The case where both beams reflected off the same surface was used as a reference to determine the position L 1 = L 2 . A measure of this interference was achieved by producing a spatial interference pattern (inset figure 12) through a slight tilting of the arms of the interferometer. Analysis of the spatial Fourier components of the resulting interference pattern provided a quantitative measure of the interference fringe contrast; a plot of contrast vs. arm position L 1 yielded a peak with width δL 1 ≈ L c . By performing this analysis while varying the path length L 1 , the plate separation was deduced. This entire procedure was then performed over a range of transverse (x, y) positions on the field plates. The resulting data are shown in figure 13.
This measurement clearly showed a bowing of the electric field plates; the plate separation varied approximately quadratically with the position in x. This is shown in the left-hand plot of figure 13. In thex direction we observed a maximum variation in the plate separation of around 20 µm. We saw a roughly † Note thatL constitutes a reversal of the supply voltages as well as a reversal of the leads connecting the power supply to the field plates, such thatẼ is unchanged. ‡ QPhotonics QSDM-680-2. 80 µm variation in theŷ (vertical) direction but note that the collimated molecular beam extended only over ±5 mm in y so the biggest plate spacing variation at a given x was ≈10 µm. From these measurements and a typical applied voltage of V = ±177 V, we expected E to vary by around 100 mV/cm in thex direction and 15 mV/cm in the y direction in the region sampled by the molecules. The indirect measurements of the spatial variation of the applied electric field provided by interferometric mapping of the field plate separation were later corroborated by direct measurements of E(x). Spatial variation of E could lead to the accumulation of geometric phases during the spin precession measurement [97]. There are known mechanisms by which such phases can contribute to eEDM-like systematic errors, as described in section 5.4, though simple estimates show that these effects are several orders of magnitude below the sensitivity of this measurement. However, additional E-field imperfections such as non-reversing fields, due to e.g. variations in the ITO coating, which could produce patch potentials, are known to contribute to eEDM-like systematic errors and are only revealed by more direct measurements of the electric field, which we will now describe.
We can write the electric field present in the precession region in the following manner: where, as usual,Ẽ = sgn(ẑ · E) is the direction of the field in the spin-precession region andL represents the binary state of the physical leads connecting the voltage supply to the field plates. The terms on the right-hand side are: EẼ, the intentionally applied electric field; E nr , a non-reversing electric field; E L , a nonreversing electric field component from the power supply that can be reversed by switchingL; and E ELẼL , a component of the applied field that is reversed by switchingẼ orL. We directly measured the components of E using the molecules themselves, in three different ways. The first method used Raman spectroscopy, driving a two-photon Lambda-type transition betweenÑ levels in |H, J = 1 as shown in figure 14. The Raman transfer was performed at positions between, but close to, the state preparation and readout laser beams, where there was sufficient optical access. The procedure was as follows: first, anx-polarised state preparation laser beam depleted a superposition |B(x,Ñ = +1,P = +1) (recall |B is the bright state as defined in section 3.1) by exciting it to the C state. Next, at a point downstream, two co-propagating,x-polarised Raman beams were used to repopulate this depleted superposition by driving population from the otherÑ state, via the transition |B(x,Ñ = −1,P = +1) → |C,P = +1 → |B(x,Ñ = +1,P = +1) . The frequencies of the two Raman beams were tuned with a pair of AOMs. The state readout laser then addressed the same transition as the preparation laser and excited the repopulated superposition to the C-state from which it spontaneously decayed back to X and fluoresced at 690 nm. The single-photon detuning is given by ∆ + δ/2 and the two-photon detuning is given by δ/2. D 1 |E| is the magnitude of the Stark shift due to the applied electric field. Right: Example scans for oppositeẼ states obtained by varying the two-photon detuning δ/2 and observing fluorescence, with Gaussian fits to the data.
Efficient transfer of population between the twoÑ states occurred for zero two-photon detuning (δ/2 in figure 14). This condition was indicated by a peak in fluorescence, giving a measure of the Stark shifted energy, and hence the absolute size of the applied field, |E|. This procedure was repeated for different positions of the Raman laser beams along thex direction. The non-reversing component of the electric field was found by repeating the measurement after reversing the applied voltages. An example of such a pair of scans is shown on the right of figure 14.
Using this method we measured the electric field at x positions where there was sufficient optical access, i.e. near the state preparation and readout laser beams. TheẼ-correlated two-photon detuning δ E = 2π × 13 kHz (2π × 11 kHz) allowed us to extract a value of the non-reversing electric field component, in the state preparation (readout) region. We did not observe any significant variation within the individual regions. We also observed that this non-reversing component did not vary with the size of the reversing electric field.
The second method used to measure the electric field had the greatest utility because it allowed for spatially resolved measurements along x in the spin precession region with comparable precision to the Raman method without perturbing the experimental apparatus. This was achieved via microwave spectroscopy. A schematic of the experimental setup is shown in figure 16.
The measurement procedure began with optical pumping of molecules into the H-state. The molecules travelled through the spin-precession region until it was entirely occupied by the molecule pulse. At this time, a π-pulse of microwaves at 2π × 39 GHz with nominalŷ polarisation was applied counter-propagating to the molecule beam. When on resonance, this transferred population from |B(ŷ,Ñ ,P) to |H, J = 2, M = 0,P (excitation to (from) eitherP (Ñ ) state was permitted) as shown in figure 15. State readout was performed as usual (see section 3.1) by optically pumping with alternating polarisationsx andŷ. The measured asymmetry (as defined in equation 11) served as a measure of the microwave transfer efficiency. The x position of the molecules at the time of the microwave pulse was mapped onto their arrival time in the detection region 2π × 39 GHz Figure 15: The transition driven by microwaves during a measurement of the electric field. We usedŷpolarised microwaves of frequency 2π × 39 GHz to drive a rotational transition between |H, J = 1 and |H, J = 2 . The M = 0 levels are labelled with their parity. We applied a moderate E-field such that  and, with knowledge of the longitudinal molecular beam velocity, v , could be extracted. Thus, the spatial dependence of the resonant frequency, ω MW (x), was provided by the time-dependence of the asymmetry, A(t). Due to the DC Stark shift, ω MW was linearly proportional to the electric field magnitude and |E(x)| could be directly extracted.
We observed a resonance linewidth of ≈2π × 25 kHz ≈ 2π/T which was limited by the microwave π-pulse duration of T = 40 µs. With our signal-to-noise, we were able to fit the resonance centre to a precision of ∼2π ×1 kHz, typically using ∼50 detuning values and averaging over ∼50 molecule pulses per detuning value.
Example data obtained via this method are shown in figure 17. In these data, it is evident that the resonant frequency of the microwaves varied across the molecule pulse by around 2π × 60 kHz. The position x of the molecules at the time of the microwave pulse was assumed to be linearly related to the molecule arrival time in the state readout region. The observed spatial variation of E was roughly consistent with expectations based on the measured variation of the plate spacing described above.
By switchingÑ andẼ between measurements of the E-field we were able to extract E nr from thẽ NẼ-correlated component of ω MW . These measurements, shown in figure 18, were used to evaluate the corresponding systematic error in equation 89.
We clearly saw a non-uniform E nr across the spin precession region. The spatial variation shown in figure 18 was reproducible for the period of several weeks over which these measurements of the electric field were taken. We are unsure as to the origin of the E nr but believe it may have been caused by patch potentials [98] present on the electric field plates. We observed unexplained disagreement between the two measurement methods (Raman spectroscopy vs. microwave spectroscopy), but note that both report non-reversing fields of a few mV/cm with the same sign.
The mapping between arrival time in the detection region and x position during the microwave pulse was approximate, suffering from spatial averaging due to a variety of effects. For example, velocity dispersion led to averaging of dx × σ v /v , where σ v is the longitudinal velocity spread of the molecular beam and dx is the distance between microwave interrogation and state readout. This averaging distance was largest, ≈1.6 cm, at the state preparation region. Spatial averaging also occurred across the ≈0.7 cm distance traversed during the T = 40 µs microwave pulse. Finally, there was averaging of the spatial position of the molecules due to the finite size of the state readout laser beam and the polarisation switching; molecules were optically pumped (with varying probability) throughout the ≈0.5 cm wide laser beam. In addition to spatial averaging, uncertainty in the mean longitudinal velocity also contributed an uncertainty in position. Changes of ≈10 m/s between molecule pulses were quite typical over the course of the E-field measurement, giving an estimated position uncertainty of 1 cm.
By adding the above contributions in quadrature we concluded that the range of positions from which the microwave-induced signals could have originated increased from around ≈1.3 cm at the state readout beam to ≈2.1 cm at the optical pumping beam. These ranges are shown as horizontal error bars at the extrema of position in figure 18.
We used a third method to measure E and E nr in situ throughout the eEDM dataset by performing 'intentional parameter variation' tests with large ∆ prep (denoted by 'c' in figure 5). Detuning the state preparation laser resulted in a reduction in the measured contrast |C| as shown in figure 24 (B). Setting ∆ prep ≈ ±2 MHz gives |C| ≈ 0.5, and the contrast was then approximately linearly proportional to ∆ prep with a sensitivity of about 1/γ C ≈ 1/(2π × 2 MHz). Any variation in the electric field would change the Stark shift, and thus also ∆ prep , resulting in a change in contrast. Thus, using the previously described spin precession scheme, we indirectly measured parity components of the electric field from the appropriate parity components of the contrast: We looked for variation of E or E nr every 3-4 hours. Measurements of E nr were consistent with the microwave measurements, with a constant value E nr (x prep ) = −4.8 ± 0.9 mV/cm. However, the mismatch ∆ N = D 1 E − ω N L between the Stark shift D 1 E and theÑ -correlated laser frequency shift, ω N L , was found to drift significantly on the scale of around 2π × 20 kHz/day. This drift of ∆ N was servoed by tuning ω N L after each measurement, ensuring ∆ N < 2π × 30 kHz at all times [92], see sections 5.2.6 and 5.6 for more details.

Magnetic Fields
Our experimental scheme did not require the application of a magnetic field. This was not the case with some previous eEDM experiments, where the magnetic field was used to define a quantization axis [40,41], or to cause the precession of spin to a direction associated with maximum sensitivity [42,99]. Instead we used the electric field to define a quantization axis, and we used the relative polarisations of the state preparation and readout lasers to define the basis in which we read out the electron's spin precession with maximal sensitivity.
However, we regularly applied a magnetic field B in order to perform searches for systematic errors. The phase accumulation induced by an eEDM δd e ≈ 5 × 10 −29 e · cm would have the same size as a Zeeman phase produced by a magnetic field of B ≈ 0.2 µG, which is small compared to some of the magnetic field imperfections in the experiment. However, phases associated with magnetic-field-induced precession were distinguished from eEDM-induced precession by the use of the switches at our disposal (e.g. electric field reversal). Nevertheless, it was important to investigate, quantify and minimize the effects of such magnetic fields, as they could have coupled with other experimental imperfections to give eEDM-like phases.
Under normal operating conditions we ran the experiment at three different magnetic field magnitudes, corresponding to a relative precession phase of φ B ≈ q π 4 for q = 0, 1, 2. The required z-component of the field was then B z = qB 0B , where B 0 = π 4 1 g1µBτ ≈ 20 mG. We also had the ability to apply transverse magnetic field components alongx andŷ, and all five linearly independent first-order gradients. The various coils that we used are illustrated in figure 19. Figure 19: A schematic of the magnetic field coils used. The main coils consisted of rectangular cosine coils (orange) wound on the surface of a cylindrical plastic frame together with additional end coils (red) to correct for the low aspect ratio (length/diameter) in our system; a second set of these coils, mirrored in the xy plane, is not shown. Also wrapped around this frame are a pair of circular auxiliary coils shown in yellow. The other auxiliary coils are shown in blue and green and consist of rectangular coils above and below the vacuum chamber. See the main text for descriptions of the functions of all of the coils.
The primary magnetic field, B z , was produced by two sets of rectangular coils, shown in orange in figure 19. These were wound on the surface of two hemicylindrical plastic shells, on the ±z sides of the spinprecession region. The coils were designed to maximize field uniformity and minimize distortion due to the boundary conditions imposed by the magnetic shielding. It was also possible to apply a ∂B z /∂z gradient with these coils. Two end coils (red in figure 19), located on the ±x ends of the spin-precession region, enhanced the uniformity of the B-field along x and enabled application of a ∂B z /∂x. The main coils were powered by two separate commercial power supplies †, and the end coils were powered by custom power supplies. The current flowing through these coils was continuously monitored throughout the course of the experiment by † Krohn-Hite 521/522 Table 1: A summary of the magnetic fields and magnetic field gradients that we could produce. The coil colours refer to figure 19.
Coil colour Fields produced Field gradients produced measuring with a digital multimeter the voltage dropped across precision resistors. We used three sets of auxiliary magnetic field coils in systematic error searches. A pair of circular Helmholtz coils (yellow in figure 19) were wrapped around the same frame used for the main coils and were formed from ribbon cable. They provided a magnetic field in the ±x directions and could also provide a ∂B x /∂x. Above and below the spin-precession region chamber (±y) there were four sets of rectangular coils (blue and green in figure 19). These allowed us to produce a field in the ±ŷ directions as well as all three associated first-order gradients. Note that the three first-order magnetic field gradients that we could not apply could be inferred from Maxwell's equations. A summary of the fields that we could apply is given in table 1.
Several measures were taken to minimize stray magnetic fields affecting the molecules. The simplest was to ensure no magnetized objects were placed within the spin-precession region. To ensure this, all components were fabricated from non-magnetic materials (e.g. no stainless steel). The magnetization of all objects was also checked before installation by passing them across an AC-coupled magnetometer sensitive to 0.1 mG field variations.
The ambient B-field in the laboratory was dominated by that from the Earth's core (∼500 mG approximately alongx +ŷ). To suppress this and other DC/low-frequency fields, the spin-precession region was surrounded by a set of five concentric cylindrical magnetic shields constructed from ≈1.6 mm thick mumetal †. Each layer of shielding should have provided around a factor of 10 reduction in the DC magnetic field [67]; however, residual magnetisation of the mu-metal was found to limit the field components to 20 µG for B x and B y , and 500 µG for B y . ‡ Each shielding layer was divided into two half-cylinders and two end caps. The outermost (innermost) shield was 132 cm (86 cm) long and had a diameter of 107 cm (76 cm). These shields had holes to allow lasers to pass through in the z direction, and to accommodate the molecule beam. There were also holes for the light pipes to extract molecule fluorescence, and some electric connections, in the x direction. Measurements and simulations showed that these holes had a negligible impact on the shielding efficiency. The shielding factor remained approximately constant up to an AC frequency ∼ 2π × 3 GHz for which the wavelength becomes comparable to the size of any apertures in the shields, ∼10 cm, and the magnetic field noise starts to penetrate the shields. However, our measurement was only sensitive to magnetic field noise at frequencies up to roughly the inverse of the spin precession time 1/τ ≈ 2π × 1 kHz [77]. The aluminium vacuum chamber also shielded AC magnetic noise above a frequency ∼1/πσt 2 µ ≈ 2π × 100 Hz, where σ ≈ 3.5 × 10 7 S/m is the electrical conductivity, t ≈ 1 cm is the thickness and µ is the permeability ≈µ 0 , the vacuum permeability [100,101].
The relatively large (B ∼ 10 mG) fields applied by the B z coils caused the inner magnetic shields to become slightly magnetized, inducing a non-reversing magnetic field, B nr ≈ 30 µG. In order to suppress this remanent field we performed a degaussing procedure on the magnetic shields by passing a 200 Hz sinusoidal current through sets of loosely wound ribbon cable coils which wrapped axially (in the xy plane at z = 0) between the shield layers. The maximum current amplitude was 1 A, sufficient to drive the mu-metal to saturation, and the amplitude was decreased with an exponential envelope over a period of 1 s. To fully degauss all layers of the magnetic shielding takes around 4 s. There was also a 1 s period of 'dead time' during which the main magnetic field was turned back on and allowed to settle. This degaussing procedure † Amuneal Inc. ‡ We later found that the residual By could be reduced to a level comparable to Bx and Bz by performing degaussing with a higher current. was repeated every time the applied magnetic field was changed, which occured approximately every 40 s.
Variations in the magnetic fields present were continuously measured throughout the experimental procedure. This was achieved using a set of four three-axis fluxgate magnetometers †, which were mounted in a tetrahedral configuration outside the spin-precession region vacuum chamber (but inside the magnetic shielding). We also used an additional fluxgate magnetometer which was positioned at a distance of around 1 m from the apparatus and outside of the magnetic shielding. By continuously recording the measurements provided by these magnetometers we were able to search for correlations of our data with the magnetic field present. In particular, we checked for the presence of a magnetic field correlated with the electric field, B E , which would have been characteristic of a leakage current flowing between the electric field plates -an effect known to contribute a significant systematic error in previous eEDM experiments [68,99]. The left-hand plot shows the reversing components of field whilst a nominal B z was applied. The right-hand plot shows the corresponding non-reversing components. The data are fit by polynomial curves.
Additional measurement of the magnetic fields was carried out by opening the vacuum system and passing a rotatable flux-gate magnetometer into the chamber. This allowed for measurement of the fields directly along the beam line. The freedom to rotate the magnetometers was crucial to distinguish between electronic offsets and B nr for fields 1 mG. From these measurements we were able to directly characterise most of the magnetic fields and first-order field gradients, including non-reversing components. Example data obtained from these measurements are shown in Figure 20. We saw that the applied fields were all flat to within 1 mG, and the non-reversing components, with the exception of B nr y , were less than 50 µG. Systematic uncertainty due to these fields is discussed in section 5.7.

Fluorescence Collection and Detection
As previously described, our experimental data consisted of laser-induced molecule fluorescence, emitted in all directions (with a well-defined angular distribution [62]) when the molecules were interrogated by the state readout laser beam. The apparatus for collecting this light is illustrated in figure 21. The fluorescence light passed through the transparent electric field plates, whose inner (outer) faces are ITO (anti-reflection [AR]) coated. Behind each field plate was a set of four AR-coated lens doublets, which collimated and then focussed the light. The optical axes of the doublets intersected a ray path from the centre of the fluorescing The lens assembly was mounted on rails and the entire assembly sat on a breadboard which was fastened to the vacuum chamber. The view on the right also shows the approximate position of the state readout laser beam as it passes through the apparatus. molecule region, accounting for refraction through the electric field plates. The first (second) lens of each doublet was a 75 mm † (50 mm ‡) diameter spherical lens of focal length 50 mm (35 mm). On each side (±z), each of the four lens doublets focussed light onto one of four sections of a 'quadfurcated' fiber bundle § whose input ends were 9 mm in diameter and fastened in lens tubes. The output of the fiber bundle was connected to a 19 mm diameter fused quartz light pipe with optical couplant gel in between. The light pipe passed out of the spin-precession region vacuum chamber and magnetic shields and directed the light onto a PMTP. Bandpass filters + were used to suppress backgrounds from e.g. scattered light. Detailed tests of the light collection were carried out [92] which estimated that ≈14 % of the fluorescence photons were collected. The major contributions to this efficiency were the finite solid angle subtended by the collection lenses (≈50 %), finite coupling efficiency into fiber bundles (≈60 %) and finite coupling efficiency between the fiber bundles and the light pipes (≈50 %). In addition, the quantum efficiency of the PMT's was specified to be ≈10 %, which further reduced the signal obtained.

Data Acquisition
The data acqusition system performed the following three functions: (i) Digital modulation of the experimental parameters necessary for acquiring the complete set of phase and contrast measurements required to extract the eEDM, as described in section 3.1.2.
(ii) Rapid (5 MSa/s) acquisition and storage of high-bandwidth fluorescence waveforms for the spin precession measurement.
(iii) Monitoring and logging of experimental parameters useful for checking the experimental state and for searching for systematic errors (e.g. magnetic fields, beam source temperatures). Data acquisition timing was controlled by a digital delay generator. † Every 20 ms, a TTL signal was produced which triggered the ablation laser Q-switch, in turn creating a pulse of molecules. Molecule fluorescence signals, measured as a PMT photocurrent, were captured on a 20-bit digital oscilloscope ‡. The oscilloscope was triggered 6-7 ms after the ablation pulse, depending on the current molecule beam forward velocity, and recorded a 9 ms window of signal containing the entire molecule signal (1-2 ms) and several ms of background. The 100 kHz square wave that drove the fast polarisation switching of the state readout laser was synchronised with the 50 Hz Q-switch trigger so that the relative phase was fixed. The 5 MSa/s data rate of the oscilloscope enabled resolution of the time-dependent structure within each 5 µs polarisation bin; this structure could vary on timescales as short as the C-state lifetime 1/γ C ≈ 500 ns [65].
Signal waveforms, S(t), were captured from two PMTs -note that we were not counting individual photoelectrons, but instead amplified and read out a voltage proportional to the count rate. These waveforms were then transferred to the control PC where they were digitally averaged over 25 pulses to form one 'trace'. The traces were then written to a hard drive. A file containing auxiliary measurements was recorded synchronously with each fluorescence trace. This file included the states of the experimental switches and other auxiliary measurements such as E-field voltages, B-field currents, laser power and polarisation, magnetic field measurements, molecular beam buffer gas flow rate, buffer gas cell temperature, and the temperature, pressure and humidity in our lab. This data proved useful in searching for systematic errors as described in section 5.

Data Analysis
In this section we describe the data analysis routine used to extract the eEDM value, and other quantities, from our dataset of nearly 10 6 PMT fluorescence traces. The entire analysis was implemented with a 'blind' offset on the eEDM channel such that the channel's mean value was not known until after all the data had been acquired and the systematic error in the measurement had been determined. No analysis changes were made after the blind was revealed. Several data cuts were applied (before removal of the blind) to ensure that the resulting eEDM measurements would be nearly normally distributed and to filter data that was not taken under normal operating conditions.

Signal Asymmetry
As described in section 3.1, the accumulated phase Φ was read out by resonantly addressing the H → C transition with linearly polarised light and monitoring the resulting fluorescence. The state readout laser was switched between orthogonal polarisations,X andŶ , at 100 kHz (with 1.2 µs of dead time between polarisations) in order to normalize against molecular flux variations. By switching at a rate fast enough that each molecule experienced both polarisations, we achieved nearly photon-shot-noise-limited phase measurements [62]. With a sufficiently wide laser beam, all molecules were completely optically pumped by both laser polarisations during their ∼20 µs fly-through time. We induced approximately one fluorescence photon from each molecule by projecting the molecule state onto the two orthogonal spin states excited by laser beams with orthogonal polarisations.
The rapid switching of the laser polarisation resulted in a modulated PMT signal, S(t), as shown in figure 22. For the following discussion we consider the polarisation state to switch at a time t = 0. Immediately after, there is a rapid increase in fluorescence as the molecules in the laser beam are quickly excited; while Ω r t 1, where Ω r ∼ 2π × 1 MHz is the Rabi frequency on the H to C transition, the fluorescence increases as S(t) ∝ Ω 2 r × t 2 . At later times, when Ω r t 1, population is about evenly mixed between the H and C states (since Ω r γ C ); hence, S(t) decays exponentially with a time constant of roughly 1/(2γ C ) ≈ 1 µs. Molecules that were not present at t = 0 continue to enter the laser beam, causing S(t) to approach a steady state. The laser is then turned off and the signal decays exponentially with time constant 1/γ C ≈ 0.5 µs. The next laser pulse, with orthogonal polarisation, is turned on 1.2 µs ≈ 2.5/γ C after the end of the previous one to prevent significant overlap of contributions to S(t) induced by different polarisations. A low-pass filter in the PMT voltage amplifier with a cut-off frequency of 2π × 2 MHz removed † SRS DG645. ‡ National Instruments PXI-5922. any short timescale dynamics from S(t), and prevented aliasing of high frequency components in the signal given our fixed digitization rate of 5 MSa/s.
To determine the fluorescence F (t) produced by each polarisation state, we subtracted a time-dependent background, B(t ), taken from data with no molecule fluorescence present, i.e. F (t) = S(t)−B(t ). Examples of the extracted F (t) and B(t ) time series are shown in figure 22A and B, respectively. B(t ) was modulated in time due to scattered light from the state readout laser beam and has a DC electronic offset intrinsic to the PMTs. The first millisecond of data, which contains no fluorescence, was used to determine B(t ). We assumed that B(t ) was periodic with the switching of the laser polarisation but did not depend on the polarisation; we inferred its value by averaging together the recorded PMT signal across all polarisation bins for ≈1 ms of data taken before the arrival of the molecule pulse. Since molecule beam velocity variations caused jitter in the temporal position of the molecule pulse within the trace, 9 ms of data were collected per pulse, despite the fact that only the ≈2 ms of strong signal with F (t) B(t ) and ≈1 ms of background contained useful information for the spin precession measurement.
Integrating F (t) over times associated with pairs of orthogonally polarised laser pulses resulted in signals F X , F Y . The integration was performed over a specified time window that we denoted as a 'polarisation bin'. Figure 26B shows two typical choices of polarisation bin and illustrates that the extracted eEDM is not significantly affected by this choice. Figure 35 shows that most of the extracted quantities did not vary linearly within the polarisation bin (Pol. Cycle Time Dependence column).
After polarisation binning, the data displayed a fluorescence signal modulated by the envelope of the molecule pulse, as in figure 22C. Figure 22D shows the asymmetry, A, computed from these data. The asymmetry is computed for each 10 µs polarisation cycle, so that for the i th cycle we have The molecule phase, and hence asymmetry (see equation 11), had a linear dependence on the time after ablation because the molecules precessed in a magnetic field over a fixed distance; the slower molecules, which arrive later, precessed more than the faster molecules, which arrived earlier. We applied a fluorescence signal threshold cut of around F = (F X + F Y )/2 ≥ 3 × 10 5 s −1 , indicated by dashed lines in figure 22C,D. Section 4.3 describes the threshold choice in detail.
To determine the statistical uncertainty in A, n ≈ 20-30 adjacent asymmetry points were grouped together. For each group, j, centred around a time after ablation t j , we calculated the mean,Ā j , and the uncertainty in the mean, δĀ j , depicted as red points and error bars in figure 22D. For smaller n, the variance in the sample variance in the mean grows, in which case, error propagation that utilises a weighted mean of data ultimately leads to an understimate of the final statistical uncertainty [102]. For larger n, the mean significantly varies within the group due to velocity dispersion, and the variance in the mean grows in a manner not determined by random statistical fluctuations. For the range n = 20-30 we observed no significant change in any quantities which were deduced from the measured asymmetry.
As described earlier in this section, the background, B(t ), which we subtracted from the PMT signal, S(t), was observed to be correlated with the fast switching of the readout laser beam polarisation. This can arise, for example, if the two polarisations have different laser beam intensities or pointings. We chose to use a polarisation independent B(t ) by averaging over the two polarisation states. This produced an asymmetry offset as per equation 31 and hence a significant Φ nr associated with the polarisation-dependent background. However, we did not consider Φ nr to be a crucial physical or diagnostic quantity. We found that this methodology produced accurate estimates of the uncertainties of quantities computed from the measured asymmetry, as verified by χ 2 analysis of measurements of Φ N E . We also found that none of the phase channels of interest changed significantly dependent on whether a polarisation-dependent B(t ) was used.

Computing Contrast and Phase
To compute the measured phase Φ we must also measure the fringe contrast C and relative laser polarisation angle θ = θ read −θ prep , as described in section 3.1. TheX andŶ laser polarisations were set by a λ/2 waveplate and were determined absolutely by auxiliary polarimetry measurements [65]. The contrast, defined as either   2C = −∂A/∂θ or 2C = ∂A/∂φ †, can be determined by dithering either the accumulated phase φ (by varying B z ) or the relative laser polarisation angle θ. We chose the latter as it could be changed quickly (< 1 s) by rotating a half-wave plate with a stepper-motor-driven rotation stage. Figure 23 shows the asymmetry as a function of θ, for a range of values of applied magnetic field. We ran the experiment at the steepest part of the asymmetry fringe (where θ = θ nr ) and measured the contrast, C j , for each asymmetry group,Ā j , by switching θ between two angles, θ = θ nr + ∆θθ, forθ = ±1 and ∆θ = 0.05 rad: To stay on the steepest part of the fringe, we chose θ nr = 0 rad for B = ±20 mG and θ nr = π/4 rad for B = ±1, ±40 mG. For these data |C| < 90% due to low preparation laser power; typically, however, |C| ≈ 95%. Solid lines represent the expected behaviour for a given magnetic field and contrast.
Because the fringe contrast was fairly constant over the duration of the molecule pulse (figure 24A), we used a weighted average ‡ of all C j measurements within the cut region for that trace to extract the accumulated phase. We also performed the analysis by fitting C j to a 2nd-order polynomial as a function of time during the ablation pulse; this led a better fit to the data, but had no significant effect on the results. We typically found |C| ≈ 95%. We believe that this was limited by a number of effects including: imperfect state preparation/readout, decay from the C state back to the H state and dispersion in the spin precession. We also observed that this value was constant over a ±2π × 1 MHz detuning range of the state preparation laser (figure 24B), indicating complete optical pumping over this frequency range. Recall that, as defined, C can be positive or negative, depending on the sign of the asymmetry fringe slope (see figure 23, or equation 32). Given that we worked near zero asymmetry where the fringe slope was steepest, and that θ nr was always chosen to be 0 or π/4, we computed the total accumulated phase as Here, q = 0, ±1 or ±2, corresponds to applied magnetic fields of ±1, ±20, and ±40 mG, respectively. We chose to apply a small magnetic field, B = 1 mG when operating at q = 0 rather than turning off the † Recall that in practice we consider C as an unsigned quantity for the purposes of data analysis. ‡ Each C j measurement is weighted by its computed uncertainty.
magnetic field completely so that we would not need to change the experimental switch sequence or data analysis routine for data taken under this condition. Figure 23 illustrates the correspondence between θ nr and applied magnetic field needed to remain on the steepest part of the asymmetry fringe.

Accounting for Correlated Contrast
It was possible for the magnitude of the contrast |C| to vary between different experimental states. For example, if the state preparation laser detuning or fluorescence signal background were correlated with any of the block switchesÑ ,Ẽ, orB, then contrast would also be correlated with those switches. As described in section 5, we observed bothÑ -andÑẼ-correlated contrast. The latter was particularly troubling since it could lead to a systematic offset in the measured eEDM if not properly accounted for: since A = −C cos[2(φ − θ)], a nonzero A N E could occur due to either C N E or φ N E . We accounted for contrast correlations by calculating C separately for each combination ofÑ ,Ẽ, andB experimental states ('state-averaged' contrast †): As previously discussed, we averaged or applied a quadratic fit to all C j (Ñ ,Ẽ,B) within a molecule pulse to computeC j (Ñ ,Ẽ,B). The precession phase was calculated from each state-specific asymmetry and contrast measurement (cf. equation 32): whereĀ is the average asymmetry over the twoθ states in a data block that share identical values ofÑ ,Ẽ andB. By construction, phases computed from state-averaged contrast are immune to contrast correlations. We also computed phases by ignoring contrast correlations (i.e. treating contrast as independent ofÑ ,Ẽ,B) and the result did not change significantly.

Computing Phase and Frequency Correlations
After extracting the measured phase Φ j (Ñ ,Ẽ,B), we performed the basis change described in equation 23, from this experiment switch state basis to the experiment switch parity basis, denoted by Φ p j , where p is a placeholder for a given experiment switch parity.
We observed that the molecule beam forward velocity, and hence the spin precession time τ j , fluctuated by up to 10% over a 10 minute time period. Since B z and g 1 are known from auxiliary measurements to a precision of around 1%, we were able to extract τ j from each block from the Zeeman precession phase measurement, Φ B j = −µ B g 1 B z τ j (see section 3.1.2). Velocity dispersion caused τ j to vary across the molecule pulse with a nominally linear dependence on time after ablation, t, however we observed significant deviations from linearity. Thus, we fit τ j to a 3rd order polynomial in t in order to evaluateτ j . Then, we evaluated the measured spin precession frequencies defined as for all phase channels p (see equation 22 for definition). We extracted the eEDM from ω N E j , which in the absence of systematic errors would be given by ω N E = −d e E eff independent of j.
From here on we will drop the j subscript that denotes a grouping of n adjacent asymmetry points about a particular time after ablation t j ; it is implicit that independent phase measurements were computed from many separate groups of data, each with different values of t j across the duration of the molecule pulse. At the end of the analysis, and whenever it was convenient to do so, we implicitly performed weighted averaging across the j subscript. Reproduced with permission from [57] Other phase channels couuld be used to search for and monitor systematic errors, discussed in detail in section 5, or to measure properties of ThO, as is the case with ω N B . We discuss the latter case here. This channel provided a measure of ∆g, the magnetic moment difference between upper and lowerÑ -levels, arising from perturbations due to other electronic and rotational states [57,76]. Because this difference limits the extent to which theÑ reversal can suppress certain systematic errors [67], it is an important quantity both in our experiment and in other experiments measuring eEDMs in molecules with Ω-doublet structure [103]. Figure 25 illustrates an observed linear dependence ∆g/2 = ηE, as predicted [57,76]. Since E and B z are precisely known from auxiliary measurements, the constant η can be directly calculated from our angular frequency measurements: Our measured value of η = −0.79 ± 0.01 nm/V was approximately half of what one would compute using the methods developed to understand the effect in the PbO molecule [76,104]. This discrepancy was subsequently understood as being primarily due to coupling to other fine-structure components in the 3 ∆ manifold [57,64].
The ω N B channel illustrates the importance of understanding phase channels besides that corresponding to the eEDM.

Data Cuts
Three data cuts were applied as part of the analysis: fluorescence rate threshold (see section 4.1), polarisation bin (see below), and contrast threshold (see below). These cuts made sure that we only used data taken under appropriate experimental conditions (e.g. only when lasers remained locked etc.) and thus ensured a high signal to noise ratio for the data used to extract the eEDM value. We thoroughly investigated how each of these cuts affected the calculated eEDM mean and uncertainty.
As previously mentioned, a fluorescence threshold cut of about F cut = 3 × 10 5 s −1 , was applied to each trace (average of 25 molecule pulses) to ensure that the fluorescence rate would always be larger than the background rate. This threshold was chosen to include the maximum number of asymmetry points in our measurement while also excluding low signal-to-noise asymmetry measurements that would increase the overall eEDM uncertainty, as described below. We also removed entire blocks (complete sets ofÑ ,Ẽ,B,θ) of data from the analysis if any of the block's experiment states had 0.5 ms of fluorescence data above F cut .
The count rates of uncorrelated fluorescence photoelectrons exhibit Poissonian statistics. In each block we averaged together four traces with the same experimental configuration. After such averaging, the number of detected photoelectrons within a pair of laser polarisation bins was 50, which was large enough that the photoelectron number distribution closely resembled a normal distribution. Because the asymmetry was defined as a ratio of two approximately normally distributed random variables (F X − F Y and F X + F Y ), its distribution was not necessarily normal. Rather, it approached a normal distribution in the limit of large F X + F Y [64]. The same followed for all quantities computed from the asymmetry, including the eEDM. The fluorescence threshold cut therefore ensured that the distribution of eEDM measurements was very nearly normally distributed. Including low-signal data would have caused the distribution to deviate from normal and increase the overall uncertainty. To check that this signal size cut did not lead to a systematic error in our determination of d e , the eEDM mean and uncertainty were calculated for multiple F cut values, as shown in figure 26. If the cut was increased above 6 × 10 5 s −1 the mean value was seen to move slightly (but within the computed uncertainties), and the uncertainty to increase. However, for all plausible values of the cuts the resulting value of d e was consistent, within uncertainties, with our final stated value.
As described in section 4.1, data points within a polarisation bin were averaged together when calculating the asymmetry (cf. figure 22A). These data points were separated by 200 ns. Numbering these points from when the readout laser beam polarisation is switched, we binned points 5-20 or points 0-25, depending on the analysis routine (see section 4.4 below) when reporting our final result. The former choice was made to cut out background signal and overlapping fluorescence between polarisation states while retaining as much of the fluorescence signal as possible whereas the latter was chosen to minimize the statistical uncertainty given the lack of evidence for systematic errors that depended on time within the polarisation switching cycle. As shown in figure 26, we checked for systematic errors associated with this choice by also using several different polarisation bins to compute the eEDM. The eEDM uncertainty increased, as expected, for polarisation bins that cut out data with significant fluorescence levels, but the mean values were all consistent with each other within their respective uncertainties.
In order for a block of data to be included in our final measurement, we also required that each of the 8 (Ñ ,Ẽ,B) experiment states had a measured fringe contrast above 80%. The primary cause of blocks failing to meet this requirement was the state preparation laser becoming unlocked. This cut resulted in less than 1% of blocks being discarded. If the contrast cut was lowered, or not applied at all, the eEDM mean and uncertainty change by less than 3% of our statistical uncertainty. As with the signal threshold, if this cut threshold was increased to 90%, close to the average value of contrast, C, then a larger fraction of data was neglected and the eEDM uncertainty was seen to increase.
For all the cuts discussed, we significantly varied the associated cut and in some cases removed it entirely. The eEDM mean and uncertainty were very robust against significant variation of each of these cuts, and the cuts were chosen before the blind offset applied to the eEDM channel was removed.

Differences Between Data Analysis Routines
As a systematic error check, we performed three independent analyses of the data. Each routine followed the general analysis method described above, but varied in many small details such as background subtraction method, cut thresholds, numbers of points grouped together to compute asymmetry, polarisation bin choice, etc. The analyses differed in the polynomial order of the fits applied to both the contrast C and the precession time τ vs. time after ablation t. The analyses also differed in the inclusion of a subset of the eEDM data that featured a particularly large unexplained signal in the ω N channel. Each of the three analyses independently computed the eEDM channel and the systematic error in the eEDM channel. The uncertainties for all three routines were nearly identical, and the means agreed to within ∆ω N E < 3 mrad/s, which is within the statistical uncertainty of the measurement δω N E = ±4.8 mrad/s. The eEDM mean and uncertainty were averaged over the three analyses to produce the final result.

EDM Mean and Statistical Uncertainty
The final data set used to report our result is shown in figure 27. It consisted of ∼10 4 blocks of data taken over the course of ∼2 weeks (figure 27B); each block contains ≈20 separate eEDM measurements distributed over the duration of the molecule pulse ( Figure 27A). All ≈2 × 10 5 measurements were combined with standard Gaussian error propagation to obtain the reported mean and uncertainty. Figure 27C,D shows histograms of all measurements on a linear (C) and log (D) scale, showing the distribution agrees extremely well with a Gaussian fit. The resulting uncertainty was about 1.2 times that expected from the photoelectron shot-noise limit, taking into account the photoelectron rate from molecule fluorescence, background light, and PMT dark current. When the eEDM measurements were fit to a constant value, the reduced χ 2 was 0.996 ± 0.006 where this uncertainty represents the 1σ width of the χ 2 distribution for the appropriate number of degrees of freedom.  The distribution of ∼200,000 separate eEDM measurements (black) matches very well with a Gaussian fit (red). The same data is plotted with both a linear and a log scale. In these histograms the mean of each individual measurement was normalized to its corresponding error bar.  When computing the eEDM result, data from superblocks were averaged together. The mean could be either weighted or unweighted by the statistical uncertainty in each superblock state. Weighted averaging minimized the resulting statistical uncertainty, but unweighted averaging could suppress systematic errors that have well-defined superblock parity from entering into the extracted value for ω N E .
Due to molecule number fluctuations, each block of data had a different associated uncertainty. However, roughly equal amounts of data were gathered for the 2 4 superblock states defined by the state readout paritỹ P, field plate lead configurationL, state readout laser polarisationR, and global laser polarisationG. For the reported eEDM value, unweighted averaging (or to be precise, performing the basis change prescribed by equation 23) was used to combine data from the differentP,R,L,G experiment states, since there were known systematic errors with well-defined superblock parity that were suppressed by these switches (see, for example, sections 5.2.2 and 5.3). Note, however, that figure 28 shows that these systematic errors produced no significant eEDM shift, and that the overall uncertainty was comparable (within 10%) when the data was combined with weighted or unweighted averaging.
Unequal amounts of data were collected for the B z , E, andk ·ẑ experimental states. For example, 40% (60%) of data were gathered with the state preparation and readout laser beams pointing east (west), k·ẑ = −1(+1). To account for this, we performed state-by-state analysis of the systematic errors: the primary systematic errors (described in section 5.2.6) were allowed to depend on the magnitude of the magnetic field (though B z = 1, 40 mG were grouped together), and the pointing direction, and separate systematic error subtractions were performed for each (B z ,k ·ẑ) state. After this subtraction, the systematic uncertainties were added in quadrature with the statistical uncertainties for each state, and the data from each state was averaged together weighted by the resulting combined statistical and systematic uncertainties.
The reported statistical uncertainty was obtained via the method above assuming no systematic uncertainty. The reported systematic uncertainty was defined such that the quadrature sum of the reported statistical and systematic uncertainties gives the same value as when incorporating the state-by-state analysis. A description of the methods used to evaluate the systematic error and the systematic uncertainty in the measurement is provided in section 5.10. To prevent experimental bias we performed a blind analysis by adding an unknown offset to the mean of the eEDM channel, ω N E . The offset was randomly generated in software from a Gaussian distribution with standard deviation σ = 150 mrad/s and mean zero. The mean, statistical error, procedure for calculating the systematic error, and procedure for computing the reported confidence interval were all determined before revealing and subtracting the blind offset.

Systematic Errors
According to Sozzi [105], 'The best way to handle systematic effects is not to have any...'. We approached our experimental design in a way to be very resilient to the systematic errors that impacted eEDM experiments in the past. We performed searches for unforeseen systematic errors, observed some, developed models to understand them, and carefully quantified them in auxiliary measurements, as described in this section.
A true eEDM should contribute to theÑẼ-correlated spin precession frequency, ω N E , with a signal ω N E T that does not vary with any experimental parameter. To discriminate between a systematic error in ω N E and ω N E T , we pursued a strategy to vary a large number of experimental parameters and imperfections while closely monitoring ω N E . If ω N E changes then there must be a systematic error correlated with that varied parameter. During our search for systematic errors we varied parameters including: applied electric and magnetic fields; magnetic field gradients; molecule beam pointing; and laser beam shape, pointing, detuning, and polarisation. In addition to monitoring ω N E , we monitored the spin precession frequency, contrast, fluorescence signal, and a number of additional experimental conditions such as molecule beam velocity, vacuum pressure and room temperature. We examined the correlations of these quantities with the experiment switches to determine whether there are any spurious signals that might point to unforeseen systematic errors, or a gap in our understanding of the experiment [65].

Determining Systematic Errors and Uncertainties
In total, we varied more than 40 separate parameters during our search for systematic errors (see Table 2). These fall into two categories. Category I contains parameters P which are optimally zero; P = 0 represents an experimental imperfection. We were able to use experimental data to put a direct limit on the size of Table 2: Parameters varied during our systematic error search. Left: Category I Parameters -These were ideally zero under normal experimental running conditions and we were able to vary them significantly from zero. For each of these parameters direct measurements or limits were placed on possible systematic errors. Right: Category II Parameters -These had no single ideal value. Although direct limits on these systematic errors could not be derived, they served as checks for the presence of unanticipated systematic errors. See the main text for more details on all the systematic errors referenced. -Search for correlations between all channels of phase, contrast and fluorescence signal -Correlations with auxiliary measurements of B-fields, laser powers, vacuum pressure and temperature -3 independent data analysis routines possible systematic errors proportional to these parameters. Category II contains parameters that have no optimum value and which we could vary significantly without affecting the nature of the spin precession measurement. The variation of these parameters could reveal systematic errors and serve as a check that we understood the response of our system to those parameters, but no quantitative bounds on the associated systematic errors were derived. For each Category I parameter P , we exaggerated the size of the imperfection by a factor greater than 10, if possible, relative to the maximum size of the imperfection under normal operating conditions,P , which was obtained from auxiliary measurements. Following previous work [41,42,106], we assumed a linear relationship between ω N E and P , and extracted the sensitivity of the ω N E to parameter P , ∂ω N E /∂P . The systematic error under normal operating conditions was computed as ω N E P = (∂ω N E /∂P )P . The statistical uncertainty in the systematic error (henceforth referred to as the systematic uncertainty) δω N E P was obtained from linear error propagation of uncorrelated random variables,

Magnetic Fields
where δP is the uncertainty inP and δ∂ω N E /∂P is the uncertainty in ∂ω N E /∂P . For parameters that had been observed to produce statistically significant shifts in ω N E , such as the non-reversing electric field, E nr , we monitored the size of the systematic error throughout the reported data set during Intentional Parameter Variations (described in section 3.1.2) and deducted this quantity from ω N E to give a value of the spin precession frequency due to T-odd interactions in the H state of ThO, Most Category I parameters did not cause a statistically significant ω N E P and were not monitored. For these parameters, we did not subtract ω N E P from ω N E , but rather included an upper limit of (ω N E P ) 2 + (δω N E P ) 2 1/2 in the systematic uncertainty on ω N E T , or chose to omit this parameter from the systematic error budget altogether based on the criteria described in section 5.10.
Where applicable, we also fit higher-order polynomial functions to ω N E with respect to P during the systematic error searches. No significant increase in the systematic uncertainty was observed using such fits and hence the contributions to the systematic error budget in Table 4 were all estimated from linear fits. We note, however, that certain non-linear dependences of ω N E on P could lead to underestimates of the systematic uncertainty, for example if ω N E has a small (large) nonzero value for large (small) values of P . In efforts to avoid this, data were taken over as wide a range as possible, it is, however, always possible that such non-linear dependence is present between the parameter values for which we took data. We had no models by which non-linear dependence could manifest by variation of the parameters investigated, so we believe the procedure outlined above produced accurate estimates of the systematic errors.

Systematic Errors Due to Imperfect Laser Polarisations
The dominant systematic errors in our experiment were due to imperfections in the laser beams used to prepare the molecular and read out the molecular state. Non-ideal laser polarisations combined with laser parameters correlated with the expected eEDM signal resulted in three distinct systematic errors which we refer to as the E nr , Ω N E r , and Stark Interference (S.I.) systematic errors. In this section, we model the effects of several types of polarisation imperfections on the measured phase Φ (sections 5.2.2 and 5.2.3) and discuss the correlated laser parameters that couple to these polarisation imperfections to result in systematic errors (section 5.2.6). We then discuss how we were able to suppress and quantify the residual systematic errors in the eEDM experiment (sections 5.2.5 and 5.2.6).

Idealized Measurement Scheme with Polarisation Offsets
As described in section 3.1, the molecules initially enter the state preparation laser beam in an incoherent mixture of the two states |±,Ñ . The bright state |B(ˆ prep ,Ñ ,P prep ) is then optically pumped away through |C,P prep leaving behind the dark state |D(ˆ prep ,Ñ ,P prep ) as the initial state for the spin precession. The molecules then undergo spin precession by angle φ evolving to a final state |ψ f = U (φ)|D(ˆ prep ,Ñ ,P prep ) where U (φ) = ± e ∓iφ |±,Ñ ±,Ñ | is the spin precession operator. The molecules then enter the state readout laser that optically pumps the molecules with alternating polarisationsˆ X andˆ Y (which are nominally linearly polarised and orthogonal) between |±,Ñ and |C,P read . For each polarisation, the optical pumping results in a fluorescence count rate proportional to the projection of the state onto the bright state, F X,Y = f N 0 | B(ˆ X,Y ,Ñ ,P read )|ψ f | 2 where f is the photon detection efficiency, and N 0 is the number of molecules in the addressedÑ level. We then compute the asymmetry, A = (F X −F Y )/(F X +F Y ), dither the linear polarisation angles in the state readout laser beams to evaluate the fringe contrast, C = (∂A/∂φ)/2 ≈ −(∂A/∂θ read )/2, and extract the measured phase, Φ = A/(2C) + qπ/4. † We then report the result of the measurement in terms of an equivalent phase precession frequency ω = Φ/τ where τ ≈ 1 ms is the spin precession time, which was measured for each block as described in section 4.2.2.
Let us first consider the idealized case in which all laser polarisations are exactly linear, Θ i = π/4 for each laser i ∈ {prep, X, Y }, the angle between the state preparation laser polarisation (prep) and state readout basis (X, Y ) is π/4, θ read − θ prep = −π/4, and the accumulated phase is small, |φ| 1 (i.e. no magnetic field is applied). Under these conditions, the measured phase Φ is equal to the accumulated phase φ. Now consider the effect of adding polarisation offsets d i to each of the three laser beams such thatˆ i →ˆ i + κd i , where κ = 1 is a perturbation parameter. It is useful to cast the polarisation imperfections in terms of linear angle imperfections, θ i → θ i + κdθ i and ellipticity imperfections, Θ i → Θ i + κdΘ i where S i = −2dΘ i is the laser ellipticity Stokes parameter; these are related bŷ Note that laser polarisations can have a nonzero projection in theẑ direction, but we assume in the discussion above thatˆ i represents a normalized projection of the laser polarisation onto the xy plane. † With these polarisation imperfections in place, the measured phase Φ gains additional terms: up to second order in κ. In the eEDM measurement, we switch between two values ofP ≡P read , the parity of the excited state addressed during state readout, and we setP prep = +1, the parity of the excited state addressed during state preparation. It is worth dwelling on equation 41 for a moment. A rotation of all polarisations by the same angle leaves the measured phase unchanged: dθ i → dθ i + dθ =⇒ Φ → Φ, as expected. A deviation in the relative angle between the state preparation and readout beams, dθ prep → dθ prep +dθ and dθ X,Y → dθ X,Y −dθ, enters into the phase measurement as Φ → Φ+2dθ, but is benign so long as dθ is uncorrelated with the expected eEDM signal. The laser ellipticities affect the phase measurement only when the state readout beams differ in ellipticity, and this contribution to the phase can be distinguished from the others by switching the excited state parity,P. This last term is particularly interesting because it allows for multiplicative couplings between polarisation imperfections in the state preparation and state readout beams to contribute to the measured phase.
Although the polarisation imperfection terms in equation 41 are uncorrelated with theÑẼ and hence do not contribute to the systematic error, we will see in later sections that additional imperfections can lead to changes in the molecule state that is prepared or read-out that are equivalent to correlations dθ N E i and dΘ N E i . The framework of equation 41 is useful for understanding how these correlations result in systematic errors in the eEDM measurement extracted from Φ N E .

Stark interference between E1 and M1 transition amplitudes
In this section we describe in detail how interference between multipole transition amplitudes can lead to a measured phase that mimicks an eEDM spin precession phase. We develop a general framework illustrating how such phases depend on laser polarisation and pointing.
In an applied electric field, opposite parity levels are mixed, allowing both odd parity (E1, M2,...) and even parity (M1, E2,...) electromagnetic multipole amplitudes to contribute when driving an optical transition. These amplitudes depend on the orientation of the electric field relative to the light polarisationˆ and the laser pointing directionk. This Stark interference (S.I.) effect forms the basis of precise measurements of weak interactions through parity non-conserving amplitudes in atoms and molecules [107][108][109]. However, it can also generate a systematic error in searches for permanent electric dipole moments which look for spin precession correlated with the orientation of an applied electric field. These Stark interference amplitudes have been calculated and measured for optical transitions in Rb [110,111] and Hg [112,113], and have been included in the systematic error analysis in the Hg EDM experiment [106,114].
In this section, we consider Stark interference as a source of systematic errors in the ACME experiment. There are two important differences between molecular and atomic systems. First, molecular states such as the H 3 ∆ 1 state in ThO can be highly polarisable and opposite parity states can be completely mixed by the application of a modest laboratory electric field. Second, molecular selection rules can be much weaker than atomic selection rules: the H 3 ∆ 1 → C 3 Π 1 transition that we drive is nominally an E1 forbidden spin-flip transition (∆Σ = 1, where Σ is the projection of the total electron spin S = 1 onto the intermolecular axis), but these states have significant subdominant contributions from other spin-orbit terms [74], between some of which the E1 transition is allowed. Both of these effects significantly amplify the effect of Stark interference in molecules relative to atoms. In this section we will derive the effect of Stark interference on the measured phase Φ.
Consider a plane wave vector potential A with real amplitude A 0 , oscillating at frequency ω, that is resonant with a molecular optical transition |g → |e ,with wave vector k = (ω/c)k, and complex polarisation : The interaction Hamiltonian H int between this classical light field and the molecular system is given by: where a indexes a sum over all of the particles in the system with charge e a , mass m a , position r a and momentum p a . Typically we apply the multipole expansion on the transition matrix element between states |g and |e ; the matrix element can then be written as where E λ describes the electric interaction of order O(( k · r) λ−1 ) and M λ describes the magnetic interaction of order O(α( k · r) λ−1 ) (where α is the fine structure constant) such that where L a is the orbital angular momentum, S a is the spin angular momentum, and g a is the spin g-factor for particle of index a (see e.g. [115]). For typical atomic or molecular optical transitions, if all moments are allowed, we expect the dominant corrections to the leading order E1 transition moment to be on the order of M1/E1 ∼ α ∼ 10 −2 -10 −3 and E2/E1∼ ka 0 ∼ 10 −3 -10 −4 , where a 0 is the Bohr radius. In this work we neglect the higher order contributions beyond E2, though the effects may by evaluated by using the expansion above.
During the state preparation and readout of the molecule state, transitions are driven between the state |g = ± d ± |±,Ñ and |e = |C,P , where d ± are state amplitudes that denote the particular superposition in |H that is being interrogated. The particular d ± combination that results in M = 0 describes the state that is dark, and the orthogonal state is bright and is optically pumped away.
It is convenient to expand the Hamiltonian H int in terms of spherical tensor operators. Furthermore, the laser is only resonant with ∆M = ±1 transitions, so the spherical tensor operators with angular momentum projections other than ±1 can be reasonably omitted. In table 3, we factor the first 4 multipole operators into products of molecule and light field operators and express the molecular operators in terms of spherical tensors T λ ±1 of rank λ = 1, 2. The E1 and M1 terms consist of vector operators with λ = 1. The E2 and M2 operators are rank 2 cartesian operators which can have spherical tensor operator contributions for λ = 0, 1, 2. The rank λ = 0 components of the E2 and M2 operators, and the λ = 1 component of the E2 operator, vanish. The rank λ = 1 component of the M2 operator does not vanish, but the light field angular dependence of this operator is equivalent to E1, so we may treat it as such.
Using well-known properties of angular momentum matrix elements [60], we may write the transition matrix element in the following form, such that ε eff is the 'effective E1 polarisation' (i.e. including the effects of interference between multipole transition matrix elements is equivalent to an E1 transition with this polarisation) with the form wheren =ÑẼẑ is the orientation of the internuclear axis in the laboratory frame, a E2 (P) = c E2 (P)/( √ 2c E1 ) and a M1 = c M1 /c E1 are real dimensionless ratios describing the strength of the M1 and E2 matrix elements relative to E1, and the c coefficients are matrix elements, which are defined using the state notation |A, J, M, Ω for electronic state A, and 'E1, M1, E2' refer to the corresponding molecular operators in table 3. It is useful to define the Rabi frequency Ω r = |M| as the magnitude of the amplitude connecting to the bright state, and the unit vectorε eff corresponding to the projection of ε eff onto the xy plane,ε This completely determines the bright and dark states, which have been previously defined in equations 12 and 13 for solely E1 transition matrix elements. The odd parity E1 and even parity M1 and E2 contributions to the effective polarisation differ by a factor ofÑẼ, which is correlated with the expected eEDM signal. Expanding the effective E1 polarisation in terms of switch parity components,ε eff =ε nr eff +ÑẼd ε N E eff , and evaluating the effectiveÑẼ correlated polarisation imperfections using equation 40, we find that the bright and dark states have effective polarisation correlations given by:ẑ It is useful to use a particular parameterization of the laser pointingk and polarisationˆ to expand the expression in equation 53 in terms of pointing and polarisation imperfections. The state preparation laser k-vector is aligned along (or against) theẑ direction in the laboratory, so it is convenient to parameterize the pointing deviation from normal by spherical angle ϑ k , and the direction of this pointing imperfection by polar angle ϕ k in the xy plane, such that: k = cos ϕ k sin ϑ kx + sin ϕ k sin ϑ kŷ + cos ϑ kẑ .
We may use a parameterization for the polarisationˆ that is similar to that in equation 14, but a slight modification is required to ensure thatk ·ˆ = 0: =N −e −iθ cos Θˆ +1 + e iθ sin Θˆ −1 + zẑ (55) where N is a normalization constant that ensures thatˆ * ·ˆ = 1. With these parameterizations in place, and expanding about small ellipticities dΘ such that Θ = π/4 + dΘ, and small laser pointing deviation, ϑ k 1, we find that theÑẼ-correlated effective laser polarisation imperfections are given by: where S i = −2dΘ i describe the laser ellipticities. Hence, following equation 41, there is a systematic error in ω N E : where c i ≡ cos (2(θ i − ϕ i,k )) and s i ≡ sin (2(θ i − ϕ i,k )) describe the dependence of the systematic error on the difference between the linear polarisation angle θ i and the pointing angle ϕ i,k in the xy plane. There is another contribution to this systematic error that arises when the coupling to the off-resonant opposite parity excited state |C, −P is also taken into account. This additional contribution becomes significant when the ellipticities are comparable to or smaller than γ C /∆ Ω,C,J=1 ≈ 0.5%.
The eEDM channel, ω N E , was defined to be even under the superblock switches (includingP), hence those terms proportional toP in the equation above do not contribute to our reported result. Additionally, theG andR switches rotate the polarisation angles for each laser by roughly θ i → θ i + π/2 periodically and the resulting ω N E signal is averaged over these states. Provided that the pointing drift is much slower than the timescale of these switches, and to the extent that the laser polarisations constituting theR andG states are orthogonal, then these systematic errors should dominantly contribute to the ω N EG and ω N ER channels which were found to be consistent with zero (see Figure 35).
An indirect limit on the size of the systematic error due to Stark interference, ω N E S.I. , may be estimated by assuming a reasonable suppression factor by which the effects in ω N ER and ω N EG may 'leak' into ω N E . We monitored the pointing drift on a beam profiler and observed pointing drifts up to dϑ k ∼ 50 µrad throughout a full set of superblock states. The absolute pointing misalignment angle was not well known but was estimated to be larger than ϑ k 0.5 mrad. Hence we may estimate a conservative suppression factor dϑ k /ϑ k 1/10 by which pointing drift may contaminate ω N E from ω N ER and ω N EG . The twoR states are very nearly orthogonal, but theG states deviate sufficiently from orthogonal (see section 5.2.5) such that the leakage from ω N EG → ω N E will dominate the systematic error; we estimate a suppression factor of about c nr p /c G p ∼ s nr p /s G p ∼ 1/5. Based on the upper limits on the measured values for ω N ER and ω N EG combined with leakage from ω N ER and ω N EG into ω N E due to pointing drift, and leakage from ω N EG into ω N E due to non-orthogonality of the twoG states, we estimate the possible size of the systematic error to be ω N E S.I. 1 mrad/s. Note that the mechanism for this systematic error was not discovered until after the publication of our result [27] and hence was not included in our systematic error analysis there. Furthermore, since we did not observe this effect, this systematic error does not match any of the inclusion criteria outlined in section 5.10 and hence is not included in the systematic error budget in this paper. Since we did not understand the mechanism for this systematic error while running the apparatus, we were not able to place direct limits on the size of this systematic error. We estimate that the absolute pointing deviation from ideal was at most 5 mrad and the ellipticity of each laser was no more than S i ≈ 5%. The E1/M1 interference coefficient is a M 1 ≈ 0.1 for the H → C transition. This gives an estimate of ω N E S.I. ∼ 0.1 mrad/s before suppression due to theR andG switches. Hence, we do not believe that this systematic error significantly shifted the result of our measurement.

AC Stark shift phases
In this section we describe contributions to the measured phase Φ that depend on the AC Stark shifts induced by the state preparation and readout lasers. We describe mechanisms by which such phase contributions may arise, and we describe mechanisms by whichÑẼ correlated experimental imperfections may couple to these phases to result in eEDM-mimicking phases. Concise descriptions of some of the effects described here can be found in [64,65,92].
During our search for systematic errors as described in section 5.1, we empirically found that there was a contribution to the measured phase dΦ(∆, Ω r ) that had an unexpected linear dependence on the laser detuning, ∆, a quadratic dependence on laser detuning ∆ in the presence of a nonzero magnetic field, and a linear dependence on small changes to the magnitude of the Rabi frequency, dΩ r /Ω r , in the presence of a nonzero magnetic field, where i ∈ {prep, X, Y } indexes the state preparation and readout lasers. The coefficients we measured were α ∆ ∼ 1 mrad/(2π ×MHz), α ∆ 2 ∼ 1 mrad/ (2π × MHz) 2 and β dΩr ∼ 10 −3 . We performed these measurements by independently varying the laser detunings ∆ i across resonance using AOMs or modulating the laser power using AOMs with the set-up depicted in figure 10 and extracting the measured phase Φ. Examples of such measurements are given in figure 31.
We determined that this behaviour can be caused by mixing between bright and dark states, due to a small non-adiabatic laser polarisation rotation or Zeeman interaction present during the optical pumping used to prepare and read out the spin state. The mixed bright and dark states differ in energy by the AC Stark shift, which leads to a relative phase accumulation between the bright and dark state components that depends on the laser parameters ∆ and Ω r . We shall now derive the AC Stark shift phase that results in equation 62, under simplifying assumptions amenable to analytic calculations.
Consider a three level system consisting of the bright |B(ε,Ñ ,P) and dark |D(ε,Ñ ,P) states and the lossy excited state |C,P with decay rate γ C . For simplicity, assume that there is no applied magnetic field for the time being. In this system, the instantaneous eigenvectors (depicted in figure 29C) are |B ± ≡ ± κ ± |C,P + κ ∓ |B(ε,Ñ ,P) , |D ≡ |D(ε,Ñ ,P) , (63) and the instantaneous eigenvalues are such that the mixing amplitudes κ ± are given by The effect of the decay of the excited state (which occurs almost entirely to states outside of the three level system) may be taken into account by adding an anti-Hermitian operator term in the Schrodinger equation, |ψ = −i(H − i 1 2 Γ)|ψ , where Γ = γ C |C,P C,P| is the decay operator. This formulation is equivalent to the Lindblad master equation, that governs the time evolution of the density matrix ρ = |ψ ψ|. In practice, we implement this decay term by calculating the time evolution of the system according to H, and then making the substitution ∆ → ∆ − iγ C /2 before calculating squares of amplitudes. It is useful to work in the dressed state basis, |D , |B ± , (basis C in figure 29) because these are nearly stationary states and have simple time evolution in the case that laser polarisation and Rabi frequency are stationary. If we allow the laser polarisation to vary in time, then the dressed state basis varies in time, and the system evolves according to the Hamiltonian, where U is the transformation from time independent basis A to time dependent basis C (from figure 29), U HU † is diagonal, and −iUU † is a fictitious force term that arises because we are working in a non-inertial frame when the laser polarisation is time dependent [75].
Assuming that the polarisation is nearly linear, Θ ≈ π/4, but allowing the polarisation to rotate slightly, and allowing for a nonzero two photon detuning due to the Zeeman shift δ = −g 1 µ B B zB , the Hamiltonian in the dressed state picture is: whereχ =Θ − i(θ + δ) can be considered to be a complex polarisation rotation rate,Ω r is the rate of change of the Rabi frequency, and∆ is the rate of change of the detuning. Note that this Hamiltonian implies that the effect of a two photon detuning arising from the Zeeman shift is equivalent to that of a linear polarisation rotating at a constant rate. Figure 29: Energy level diagrams depicting the Hamiltonian when the three-level H → C transition is addressed by the state preparation or readout lasers in three different bases. Solid double-sided blue arrows denote strong laser couplings between H and C. Wiggly red arrows denote spontaneous emission from C to states outside of the three level system. Dashed orange lines denote weak couplings induced by laser polarisation rotation. Basis A is useful for describing the spin precession phase induced by the Zeeman and eEDM Hamiltonians. Basis B is useful for describing the states that are prepared and read-out in the spin precession measurement. Basis C is useful for evaluating the AC Stark Shift phases induced by laser polarisation rotations.
We may then apply first order time-dependent perturbation theory in this picture to determine the extent of bright/dark state mixing due toχ in the time evolution of the system. If we parameterize the time-dependent state as then in the case of a uniform laser fieldΩ r = 0, of duration t and with a constant detuning∆ = 0, the time evolution of the coefficients is given at first order by: In the state preparation region, the molecules begin in an incoherent mixture of the states |B(ε prep ,Ñ ,P) and |D(ε prep ,Ñ ,P) and then enter the state preparation laser beam. In the ideal case of uniform laser polarisation, molecules that were in the bright state are optically pumped out of the three level system, and molecules that are in the dark state remain there; this results in a pure state, |D(ε prep ,Ñ ,P) . However, if there is a small polarisation rotation by amount dχ ≡ t 0χ (t ) dt ≡ dΘ − i(dθ − g 1 µ B B zB t), such that |dχ| 1, then the dark state obtains a bright state admixture that may not be completely optically pumped away before leaving the laser beam. † In this case, the final state can be written as whereε prep is the effective polarisation that parameterizes the initial state in the spin precession region and Π is an amplitude that accounts for the AC Stark shift phase and the time dependent dynamics of the non-adiabatic mixing due to the polarisation rotation, The deviations between the effective polarisation and the actual laser polarisation can be viewed as effective polarisation imperfections, that lead to shifts in the measured phase Φ as described in equation 41. For definiteness, consider the case in which the polarisation rotation rateχ(t ) = dχ/t is a constant for the duration of the optical pumping pulse t. In this case, This function has the property that ImΠ is an odd function in ∆ that can take on values up to order unity across resonance (a frequency range on the order of γ C ) and is exactly zero on resonance. ReΠ is an even function quadratic in ∆ about resonance, and depends on Rabi frequency on resonance. If the laser beam intensity reduces quickly as the molecule leaves it then most of the AC Stark shift phase arises from the last Rabi flopping period before the molecule exits the laser beam (providedχ is nonzero during that time). If the intensity reduces slowly, the AC Stark shift phase can be exacerbated since the bright state amplitude is not as effectively optically pumped away while Ω r < γ C . Nevertheless, beamshaping tests shown in figure  31 and numerical simulations indicate that Π is not very sensitive to the shape of the spatial intensity profile of the laser beam or the shape of the spatial variation of the polarisation.
If we consider only the first order contribution to the shift in the measured phase in equation 41, dθ prep,eff , and neglect the second order shift that arises due to dΘ prep,eff , then we can relate the parameters in equation 62 to the amplitude Π accounting for the AC Stark shift phase and the complex polarisation rotation dχ, by We can interpret these results as follows. The linear dependence of the measured phase on detuning, α ∆,prep , comes from a spatially varying ellipticity in the x direction coupling to the AC Stark shift phase. Similarly, the quadratic dependence of Φ on ∆, α ∆ 2 ,prep , and the dependence of Φ on a relative change in Ω r , β dΩr,prep , come from either a spatially varying linear polarisation in the x direction or a Zeeman shift, each coupling to the AC Stark shift phase. Here, we only analyzed the phase shift that results from AC Stark shift effects in the state preparation laser beam, but there is an analogous phase shift in the state readout beam. There are several other subdominant effects that also contribute to the AC Stark shift phase behavior described in equation 41 in the presence of polarisation imperfections. The opposite parity excited state |C, −P couples strongly to the dark state, but the mixing between these two states is weak because the transition frequency is off-resonant by a detuning ∆ Ω,C,J=1 ≈ 2π × 51 MHz γ C . In the case that an optical pumping laser has nonzero ellipticity, the bright state gains a weak coupling to the opposite-parity excited state proportional to this ellipticity. Then, two-photon bright-dark state mixing ensues in such a way that the mixing amplitude, and hence the measured phase, depends on the laser detuning.
The rapid polarisation switching of the state readout beam can also introduce AC Stark shift-induced phases in the absence of a polarisation gradient, if the average ellipticity between the two polarisations is nonzero. Suppose a particular molecule is first excited by theˆ X polarised beam. The two bright eigenstates |B ± are mostly optically pumped away, resulting in a fluorescence signal F X . The population remaining in the bright eigenstates acquires a phase relative to the dark state, due to the AC Stark shift. Then the molecules are optically pumped by theˆ Y polarised beam. If there is a nonzero average ellipticity,ˆ Y is not quite orthogonal toˆ X and the new bright eigenstates that give rise to the fluorescence signal F Y are superpositions of the former bright and dark states that acquired a relative AC Stark shift phase. This results in a fluorescence signal, and hence measured phase component, that depends linearly on laser detuning ∆.

Polarisation Gradients from Thermal Stress-Induced Birefringence
The AC Stark shift phases described in the previous section can be induced by polarisation gradients inx across the state preparation and readout laser beams. In this section we describe a known mechanism by which these arose. Recall that these laser beams passed through transparent, ITO-coated electric field plates. For an absorbance α and laser intensity I, the rate of heat deposition into the plates isQ (x, y) = α I (x, y). The laser beam profile is stretched in the y direction to ensure that all molecules are addressed. For simplicity we assume that the heating distribution,Q (x, y) =Q (x), is completely uniform in the y direction. We also assume that there are no shear stresses, i.e. local expansion of the glass is isotropic. Under these assumptions, the relationship between the heating rate,Q, and the internal stress tensor σ ij (where i, j are Cartesian indices) is where E, α V and κ are the Young's modulus, coefficient of thermal expansion, and thermal conductivity, respectively [116]. Unit vectorsx andŷ correspond to the principal axes of the stress tensor due to the symmetry of the heating function, hence the off-diagonal (shear) elements are zero, σ xy = 0. The other diagonal component, σ xx , is uniform across the plates, and equal to σ yy far away from the laser. The stress-optical law states that the birefringence and stress are linearly proportional along the principal axes of the stress tensor [117]. The difference between the indices of refraction in the x and y directions is then ∆n = K (σ xx − σ yy ), where K ≈ 4 × 10 −6 MPa −1 is the stress-optical coefficient for Borofloat glass [118]. The retardance of an incident laser beam of index i is Γ i = 2π∆n (t/λ), where t is the thickness of the field plates (in the z direction), and λ is the wavelength of light. Hence, in this limit, the retardance due to thermal stress-induced birefringence is related to the laser intensity by: where η = 2πKEα V α/κ ≈ 26 × 10 −6 W −1 is a material constant of Borofloat glass [118]. The ellipticity imprinted on the nominally linearly polarised laser beam is given by where θ i is the linear polarisation angle and φ Γ,i is the orientation of the fast axis of the birefringent material (nominallyx in our case).
Assuming the laser has total power P , a Gaussian profile in x with standard deviation w x , and a top-hat profile in y with half width w y , the intensity is given by where 2w y w x . There is then an analytic solution to equation 82 from which we extracted a retardance gradient in the laser tail, x = w x , of for a nominal laser power of ≈2 W. Similar results were obtained from numerical finite element analysis. Thermal stress-induced birefringence has been observed in similar systems such as in UHV vacuum windows [119], laser output windows [120], and Nd:YAG rods [121]. The estimates of the ellipticity gradient agree well with measurements of the polarisation of the beam, as shown in figure 30. These polarimetry measurements were adapted from the procedure described in [122]; a polarimeter was constructed consisting of a rotating quarter-wave plate, fixed polariser, and fast photodiode. The use of a fast photodetector allows for polarimetry of the probe beam during the 100 kHz polarisation switching. The resolution of the system was such that we could quickly measure the normalized circular Stokes parameter, S, to a few percent, which is sufficient to measure typical birefringence gradients of ∼10% across the beam.

Suppression of AC Stark Shift Phases
We were able to suppress the magnitude of the AC Stark shift phases in several different ways that are illustrated in figure 31. The ellipticity gradient across the state preparation laser beam was suppressed by tuning the linear polarisation angle: as per equation 83, the ellipticity gradient is proportional to sin(2θ prep − 2φ Γ,prep ), which vanishes when the polarisation is aligned along a birefringence axis, i.e. θ prep = φ Γ,prep , φ Γ,prep + π/2. To determine φ Γ,prep we measured the total accumulated phase as a function of laser detuning for various θ prep and then extracted the slope α nr ∆,prep = ∂Φ nr /∂∆ prep for small detuning values. Note that when fitting the phase vs. detuning data we found that cubic functions provided significantly better fits over the detuning ranges used (see Figure 31(B)). We then selected θ prep to minimize α nr ∆,prep . This suppressed α nr ∆,prep by about a factor of 50 relative to its original value, to α nr ∆,prep 0.1 mrad/(2π × MHz). Another method implemented to suppress AC Stark shift phases was to reduce the time-averaged power of the state preparation laser incident on the field plates. We used a chopper wheel to modulate the laser at 50 Hz, synchronous with the molecular beam pulses, with a 50% duty cycle. We estimated the time scale for thermal changes to be on the order of Q/Q ∼ 2ρCw 2 x /κ ∼ 10 s, where ρ and C are the density and heat capacity of Borofloat respectively, so did not anticipate any significant transient effects to be introduced. This modification reduced the retardance gradient, and hence the value of α nr ∆,prep , by about a factor of two, as shown in Figure 31(C).
Finally, α nr ∆,prep was suppressed by shaping of the laser beam intensity profile. AC Stark shift phases were most significant at the downstream edge of the state preparation laser beam. Here, the intensity is such that bright-dark state mixing is still occurring but the bright state is not efficiently optically pumped away. By making the spatial intensity profile drop off more rapidly, we reduced the time that molecules spent in this intermediate intensity regime. This was achieved by taking advantage of the aspherical distortion introduced by misaligning a telescope immediately before the laser beam entered the spin-precession region. This suppressed α ∆,prep and β ∆ 2 ,prep by ≈2, as shown in Figure 31(C). In addition to a phase suppression, we noticed that the optimal laser polarisation angle changed after implementing the steps described, as can be seen in Figure 31(C). The reason for this change is not definitively known, but we suspect that as we suppressed the birefringent contribution to the AC Stark shift phase, the non-birefringent contributions (i.e. the phase due to nonzero ellipticity causing bright-dark state mixing via the off-resonant opposite parity excited state) became fractionally larger, and we needed to tune the polarisation angle to obtain cancellation between these two classes of effects.
We observed much smaller AC Stark shift phases in the state readout laser beam than in the state preparation laser beam. This is not surprising since the effect is largely birefringent; the contributions to the effective polarisation imperfections for theˆ X andˆ Y polarised lasers should be opposite in sign, dθ X ∝ sin(2(θ read − φ Γ,read )), dθ Y ∝ sin(2(θ read − φ Γ,read + π/2)), such that they cancel each other in the measured phase (cf. equation 41). The residual AC Stark shift phases measured in the state readout beam gave α nr ∆,read ≈ 0.5 mrad/(2π × MHz). This was sufficiently small that the methods of suppression described above were only implemented in the state preparation region.

Systematic Errors due to Correlated Laser Parameters
In the discussion above, we described how polarisation imperfections can lead to contributions to the measured phase that depend on the AC Stark shifts and hence on the laser detunings ∆ i and Rabi frequencies Ω r,i . However, these phases only produce a systematic error in ω N E if there is a nonzero correlation ∆ N E i or Ω N E r,i of the laser detuning or Rabi frequency. We observed such correlations and discuss them in this section. We will also describe how we evaluated the associated systematic errors.
In section 3.2.4 (see figure 11) we discussed how a non-reversing component of the applied electric field, E nr , could produce a ∆ N E . In an entirely analogous manner, the Rabi frequency magnitude Ω r of the H → C transition can exhibit the following correlations: Here, Ω nr r,i is the dominant component of the Rabi frequency for laser i ∈ {prep, X, Y }, which could fluctuate in time on the order of 5% due to laser power instability. Ω N r,i is generated by a laser power difference between theÑ states. This arose because we routed the laser light along different paths through a series of AOMs for each state. We measured this effect with photodiodes and found that the largest fractional power correlation was Ω N r /Ω nr r ≈ 2.5 × 10 −3 . An additional contribution to Ω N r,i and a contribribution to Ω N P r,i on the same order arises due to Stark mixing between rotational levels in H and C, leading toÑ -andÑP-correlated transition amplitudes on the H → C transition.
Although we did not observe a laser power correlation withÑẼ we did observe signals consistent with a Rabi frequency correlation, Ω N E r . A nonzeroÑẼ-correlated fluorescence signal (as defined in section 4.1) that also reversed with the laser propagation directionk ·ẑ, F N E /F nr ≈ −(2.4 × 10 −3 )(k ·ẑ), together with a nonzero ω N EB ≈ (2.5 mrad/s)(B z /mG)(k ·ẑ), provided the first evidence that a nonzero Ω N E r existed in our system. We believe that this fluorescence correlation arises from a linear dependence of the fluorescence signal size on Rabi frequency, F N E = (∂F/∂Ω nr r )Ω N E r , which is nonzero since the state readout transitions were not fully saturated. We believe that the signal in ω N EB was caused by a coupling between the Rabi-frequency correlation and the B-odd AC Stark shift phase, ω N EB = 1 τ β B dΩr B z (Ω N E r /Ω nr r ). We were able to verify a linear dependence of both of these channels on Ω N E r by intentionally correlating the laser intensity withÑẼ using AOMs; this is shown for the Φ N EB channel in Figure 32. Varying the size of this artificial Ω N E r allowed us to measure the value present in the experiment under normal operating conditions, Ω N E r /Ω nr r = (−8.0 ± 0.8) × 10 −3 (k ·ẑ). Ω N E r can couple to β nr dΩr,i as per equations 62 and 80 to result in a systematic error in ω N E . A nonzero β nr dΩr,i can be produced by a linear polarisation angle gradient (not observed in the experiment) or by a non-reversing Zeeman shift component g 1 µ B B nr z . While searching for a model to explain the intrinsic Ω N E r , we developed the Stark interference model presented in section 5.2.2. For unnormalized effective polarisation ε eff = ε nr eff +ÑẼd ε N E eff , this model predicts , which correctly predicts the dependence of Ω N E r on the laser propagation directionk ·ẑ. However, the factors a M1 and a E2 , which correspond to the ratio of M1 and E2 amplitudes to the E1 amplitude, must be real for a plane wave, so Im [(a M 1 + a E2 )] = 0. Hence this model fails to explain this Rabi frequency correlation unless there is some additional effect that introduces a phase shift between the E1 and M1 amplitudes. For example, interference between the E1 amplitude due to the incident laser beam, and a phase shifted M1 amplitude due to a (low intensity) reflected beam can lead to a nonzero Ω N E r by this model. However, this phase factor oscillates spatially on the scale of the light wavelength, which is very small compared to the size of the molecule cloud and hence should average out over the entire molecular beam cloud. The origin of the intrinsic Ω N E r is still not fully understood, and we are continuing to explore models to understand this effect. resulting from correlated power P N E systematically shifts ω N EB in accordance with equation 62. Φ N EB is zero when the applied P N E is such that there is no netÑẼ-correlated Rabi frequency. The intrinsic Ω N E r (i.e. that inferred when P N E = 0) changed sign withk ·ẑ within the resolution of the measurement. The slopes between the two measurements differ due to differences in the AC Stark shift phase, believed to be due to differences in the spatial intensity profile and polarisation structure between the two measurements.
Given the empirical AC Stark shift phase model in equation 62, the resulting systematic errors in the frequency measurement are given by Early in the experiment, we observed a nonzero systematic shift ω N E E nr and took the steps outlined in section 5.2.5 to suppress it. To verify that the steps taken were effective, we examined ω N E as a function of an intentionally applied non-reversing electric field. The resulting data are shown in figure 33. The original slope, ∂ω N E /∂E nr = (6.7 ± 0.4)(rad/s)/(V/cm), corresponded to a systematic shift of ω N E E nr ≈ −34 mrad/s when combined with the measured E nr ≈ −5 mV/cm. Following the modifications described above, the ∂ω N E /∂E nr slope was greatly suppressed, reducing the systematic error to ω N E E nr < 1 mrad/s, well below the statistical uncertainty in the measurement of ω N E .
Non-reversing electric -eld, E nr (mV/cm)  Figure 33: Linear dependence of the ω N E channel on an applied non-reversing electric field observed before (red) and after (black) we suppressed the known AC Stark shift phase by optimizing the preparation laser beam shape, time-averaged power and polarisation.
Because we observed that the parameters E nr and Ω N E r caused systematic errors in ω N E , we intermittently measured the size of the associated systematic errors throughout the datasets that were used for our reported result. We measured ∂ω N E /∂E nr by applying a range of large non-reversing electric fields, up to around 70 times that present under normal running conditions. The value of ∂ω N E /∂Ω N E r was measured by applying a correlated laser power P N E in the state preparation and state readout beams with a magnitude corresponding to an applied Ω N E r that was up to 20 times that measured under normal operating conditions. These parameters were measured for multiple values of the magnetic field magnitude, B z , for which different state readout laser beam polarisations were required. Due to known birefringent behavior of the AC stark shift phases, we allowed for this possibility for all AC stark shift phase systematic errors. We measured ∂ω N E /∂E nr for bothk ·ẑ = ±1, but the Ω N E r systematic error was only discovered after thek ·ẑ = +1 dataset and hence ∂ω N E /∂Ω N E r was only monitored during thek ·ẑ = −1 dataset. The Ω N E r systematic error during thê k ·ẑ = +1 dataset was determined from auxiliary measurements of the AC Stark shift phase. As described in section 3.2.5, E nr (x) exhibits significant spatial variation along the beam-line axis, x. However the E nr that was intentionally applied to determine ∂ω N E /∂E nr was spatially uniform, and hence these measurements were insensitive to the difference (E nr (x prep ) − E nr (x read )) between the state preparation laser beam at x prep and the state readout beam at x read . For this reason, we deduced the systematic error proportional to the difference (E nr (x prep ) − E nr (x read )) from auxiliary measurements of the AC Stark shift phase parameters, α nr ∆,i . In summary, the systematic errors proportional to E nr and Ω N E r that were evaluated and subtracted from ω N E to report a measured value of ω N E T can be expressed as where (∂ω N E /∂E nr ) and (∂ω N E /∂Ω N E r ) were monitored by Intentional Parameter Variations (see section 3.1.2) throughout the dataset used for our reported result, and E nr (x prep ), E nr (x read ), Ω N E r , α nr ∆,i , and β nr dΩr,i were obtained from auxiliary measurements. These two systematic errors account for almost all of the systematic offset that was subtracted from ω N E to obtain ω N E T as described in section 5.10.

A N E asymmetry effects
In addition to the dependence of the measured phase on laser detuning and Rabi frequency, we observed dependence of the asymmetry A (as defined in section 4.1) on the laser parameters ∆ read and Ω r,read , due to differences between the properties of the X and Y readout laser beams. The laser-induced fluorescence signal F (∆, Ω r ) varies quadratically with detuning (for small detuning) and linearly with Rabi frequency. Under normal conditions, the signal sizes from X and Y are comparable, F X ≈ F Y ≈ F . If the X and Y beams have different wavevectors, k X,Y = k nr ± k XY , and k XY has some component alongx, then the two beams will acquire different Doppler shifts. This leads to a linear dependence of the asymmetry on detuning, which in turn can couple to ∆ N E to result in a contribution to A N E , Similarly, if the two readout beams differ in Rabi frequency, Ω r,X/Y ≈ Ω nr r ± Ω XY r , the asymmetry becomes linearly dependent on Rabi frequency, which in turn can couple to Ω N E r to result in a contribution to A N E , However, these asymmetry effects are very distinguishable from spin precession phases and polarisation misalignments. Since theP andR switches effectively swap the role of the X and Y readout beams, the A N E effects described above do not contribute to ω N E when summed over these switches. Additionally, asymmetry effects, once converted to an equivalent frequency or phase, depend on the sign of the contrast, C, unlike true phases. In the B z ≈ 20 mG configuration, sgn(C) = sgn(B z ), but sgn(C) has no dependence on sgn(B z ) for B z ≈ 1, 40 mG. Hence asymmetry correlations A N E are mapped onto frequency correlations ω N EPR or ω N EBPR depending on the magnetic field magnitude.
If the pointing or Rabi frequency differences between the X and Y beams drift on timescales comparable to or shorter than theP orR switches, these effects can occasionally 'leak' into the 'adjacent' channels ω N EP , ω N ER , ω N EBP , ω N EBR ; however, we have not seen any evidence of these effects contributing to the ω N E channel itself, and hence did not include systematic error contributions due to these effects in our systematic error budget.

5.4.Ẽ-Correlated Phase
Previous eEDM measurements have often been limited by a variety of systematic errors that would have produced anẼ-correlated phase precession frequency in our experiment, ω E [17,38,41,50], such asẼcorrelated leakage currents, geometric phases, and motional magnetic fields. Our ability to spectroscopically reverse the molecular orientation through a choice ofÑ distinguished these effects from an eEDM-generated phase. In addition, the aforementioned effects scale with the magnitude of the applied electric field, which was orders of magnitude smaller in our experiment than previous similar eEDM experiments due to the high polarisability of ThO [41]. Also, because the molecular polarisation was saturated, the eEDM phase should have been independent of the magnitude of the applied field. We also note that any shifts from leakage currents and motional magnetic fields coupled through the magnetic dipole moment, which is near-zero in the H-state of ThO. Thus we expected ω E to be substantially suppressed, and that it should not enter ω N E at any significant level.
The reversal ofÑ did not, however, entirely eliminate an eEDM-like phase due to ω E . As discussed in section 4.2.2, there was a small and E-field dependent difference between the g-factors of the twoÑ levels [57,76], which meant that a systematic error in the ω E channel showing up in ω N E at a level given by ω N E ω E = (ηE/g 1 )ω E . We verified this relation by intentionally correlating a 1.4 mG component of our applied magnetic field withẼ. This deliberate B E resulted in a large shift in the value of ω E and a ∼1000-times smaller offset of ω N E , as illustrated in figure 34. The intentionally applied B E was the only experimental parameter that was observed to produce a measurable shift in ω E . Even large (∼20 mG) magnetic fields components alongx andŷ, which exaggerate the effect of motional magnetic fields, did not shift ω E (this is expected, since the large tensor Stark shift in |H, J = 1 dramatically suppresses the effect of motional magnetic fields [36]). For our eEDM data set, ω E was consistent with zero. We included a contribution from ω E in our error budget for ω N E by multiplying the mean and uncertainty of the extracted ω E by our measured |E|-dependent suppression factors ηE/g 1 .

5.5.Ñ -Correlated Laser Pointing
We discovered a nonzero, time-dependent signal in ω N which was associated with anÑ -correlated laser pointing,k N ≈ 5 µrad. An investigation into the mechanism behind this effect was inconclusive. We found that the pointing correlation appeared downstream of the AOMs that created the rapid polarisation switching and improved alignment was able to reduce the effect. We also found that the observed pointing was in some way correlated with the seed power and input angle of incidence into the high-power fiber amplifier immediately upstream of the polarisation switching, despite the fact that pointing out of the amplifier did not fluctuate. Since we used four different sets of AOMs to perform theÑ andP switches before the amplifier, we observed laser pointing correlated with both of these switches. By matching the characteristics of these four beam paths we were able to suppressk N to < 1 µrad.
The effect ofk N on ω N was studied by exaggerating the former with piezoelectrically actuated mirrors. Examining ∂ω N /∂k N showed significant fluctuations in its value. We were unable to identify the mechanism by whichk N affected ω N .
We had no evidence that the effect causing the observed variation in ω N also caused a systematic error in ω N E , but to be cautious we included an associated systematic uncertainty in our systematic error budget (section 5.10). Assuming a linear relationship between ω N E and ω N , we extracted ∂ω N E /∂ω N from a combination of data taken under normal conditions and with an exaggerated ω N induced by an exaggerated k N . We then placed an upper limit on a possible systematic error ω N E ω N based on the value of ω N obtained under normal running conditions. The resulting systematic uncertainty was four times smaller than our statistical uncertainty.

Laser Imperfections
Of the lasers used in our experiment, only the state preparation and readout lasers were known to produce possible systematic errors; imperfections in the rotational cooling, optical pumping or target ablation lasers simply resulted in a reduction in usable molecule flux. As part of our search for systematic errors, we intentionally exaggerated all known state preparation and readout laser imperfections possible without dismantling the apparatus (cf. table 2). In this section we describe this procedure and the resulting contributions to our systematic error budget.

Laser Detuning
The correlated components of the state preparation and readout laser beam detunings are described in detail in section 3.2.4. Each detuning component was separately exaggerated and in some cases multiple components were simultaneously exaggerated. Most of the detuning terms in equation 27 were exaggerated to ±2π × 1-2 MHz. No detuning or detuning correlation produced a significant shift in ω N E other than ∆ N E caused by E nr , discussed in section 5.2.6. In some cases, shifts in other phase channels were induced, but all shifts were consistent with well-understood AC Stark shift and asymmetry models described in sections 5.2.3 and 5.3. For example, the combination of nonzero ∆ N and ∆ nr coupled to the B-dependent component of the AC stark shift phase (equation 80) induces a significant shift in ω N B (cf. equation 62). Asymmetry correlations also resulted from these detuning correlations, but these were only manifested in channels odd with respect toP andR, and hence had no plausible effect on ω N E . Because the YbF eEDM experiment [99] observed unexplained dependence of the measured eEDM value on state preparation microwave detuning, we included a systematic error contribution from all detuning imperfections in our systematic error budget.

Laser Pointing and Intensity
Similar to detuning imperfections, the state preparation and readout lasers could have imperfect pointing and correlated intensities. Ideally the laser propagation direction,k, would have been parallel to the laboratory electric field. This would have diminished the amount ofẑ polarised light experienced by the molecules, which could drive unwanted off-resonant transitions, and prevented stray retroflection from the ITO field plate surfaces. Using this ITO retroflection as a guide, we alignedk perpendicular to the field plate surface, and therefore parallel toÊ, to within ∼ 3 mrad. To test for errors related to imperfect pointing, both the state preparation and readout pointing misalignments were exaggerated in the x-direction to ±10 mrad, as was the relative pointing of theX andŶ state readout beams. The vacuum windows and ∼3.8 cm wide holes in the magnetic shields prevented us from further misaligning the beams. To decouple pointing imperfections from detuning imperfections, the state preparation and readout laser frequencies were tuned to resonance after each pointing adjustment. No shift in ω N E was observed and no systematic error contribution from pointing imperfections was included. Pointing imperfections were only observed to affect the signal asymmetry, as previously discussed in section 5.3.
Unlike laser pointing and detuning, there was no 'ideal' value for laser intensity. The state preparation and readout laser intensities were chosen such that we were driving optical pumping to completion on the H → C transition without producing unnecessary thermal stress on the field plates. We decreased each laser intensity by a factor of four to check that there was no variation in ω N E . We observed a nonzero Ω N r caused by theÑ -correlated seed power into the high-power fiber amplifiers and by Stark mixing between rotational levels in H and C as discussed in section 5.2.6. We exaggerated this imperfection by a factor of 20. Only ω N B was shifted, consistent with our understanding of the B-correlated AC Stark shift phase. These intensity systematic error checks were not included in the systematic error budget.

Magnetic Field Imperfections
The H state is very insensitive to a magnetic field B z due to its small g-factor, as discussed in section 2.2. Sensitivity to the transverse fields is even further suppressed by the large size of the tensor Stark shift relative to the Zeeman interaction. Nevertheless, there are known mechanisms by which magnetic field imperfections can contribute to systematic errors: B nr z can contribute to the ω N E Ω N E r systematic error discussed in section 5.2.6, and transverse fields B nr x and B nr y can lead to the geometric phase systematic errors [67] discussed in section 5.4. We designed the experiment to allow a wide variety of magnetic field tilts and gradients to be applied as described in section 3.2.6 and we directly looked for systematic errors resulting from these magnetic field imperfections.
Both B-correlated and uncorrelated imperfections were applied. We did not precisely measure the residual values of each of these parameters along the molecule beam line until we had studied all systematic errors and collected our published data set. Based on the projected ∼10 5 magnetic shielding factor, we expected all stray magnetic fields and gradients to be on the order of 10 µG and 1 µG/cm, respectively. For this reason we only exaggerated these imperfections to ∼2 mG and ∼0.5 mG/cm. When we mapped out the magnetic field with a magnetometer inserted between the electric field plates as described in section 3.2.6, we discovered that several imperfections were much larger than we expected (e.g. B y ≈ 0.5 mG). This was caused by poor magnetic shielding due to insufficient shield degaussing. For this reason we gathered additional eEDM data with some magnetic field parameters exaggerated by an additional factor of five. ω N E and nearly all other frequency channels, apart from ω nr and ω B were not observed to be affected by any of these magnetic field parameters. Because uncorrelated stray magnetic fields and magnetic field gradients caused unexpected eEDM offsets in the PbO eEDM experiment [50], we included contributions from all uncorrelated magnetic field imperfections in our systematic error budget described in section 5.10.

Electric Field Imperfections
Unlike the magnetic field, we do not have the ability to control electric field gradients and stray electric fields, aside from the average value of E nr . The field plates were located at the center of the experiment, inside the vacuum chamber and magnetic shields and coils, with no direct access available. To search for systematic errors related to the electric field, equal amounts of eEDM data were gathered with two different electric field magnitudes. The ω N E values from both field magnitudes were consistent with each other. The YbF eEDM experiment observed unexplained eEDM dependence on the voltage offset common to both field plates. For this reason we exaggerated this offset by a factor of 1000 (relative to its residual value of ≈5 mV) and, even though it did not shift our eEDM measurement, included it in our systematic error budget.

Molecule Beam
The molecule beam should have ideally travelled parallel to the electric field plates and well-centred between the plates. This minimizes Doppler shifts, protects the plates from being coated with ThO, and ensures that the molecules experience the most uniform electric field. The entire beam source vacuum chamber sat on a two axis (yz) translation stage. The exit aperture of the buffer gas cell was aligned to within 1 mm of the centre of the fixed collimators and electric field plates, using a theodolite. Geometric constraints only allowed us to exaggerate the cell misalignment by roughly a factor of three (up to 3 mm) before the molecules would have hit the sides of the field plates. We also varied the transverse spatial and velocity distributions by using adjustable collimators between the beam source and spin-precession region to block half of the beam from the ±x, ±ẑ directions. The value of ω N E was not observed to shift with any molecule beam parameter adjustment.

Searching for Correlations in the eEDM Data Set
In addition to performing systematic error checks for possible variations of ω N E with various experimental parameters, we searched for statistically nonzero values within the set of 1536 possible correlations with the block and superblock switches. This analysis was performed for our primary measured quantities ω, C, and F and for a wide range of auxiliary measurements such as laser powers, magnetic field, room temperature, etc. We also examined the switch-parity channels of ω, C, and F as a function of time within the molecule beam pulse, and as a function of time within the polarisation switching cycle. We used the Pearson correlation coefficient to look for correlations between the aforementioned switch-parity channels and used the autocorrelation function to look for signs of time variation of the mean within those channels. Figure 35 illustrates data from such a search with a subset of the previously described quantities. In this search, we looked at 4390 quantities and we set the significance threshold at 4σ which correponds to a probability of p ≈ 0.25 that there will be one or more false positives above that threshold. We represented the significance of each of these quantities with a grayscale pixel. Each pixel that was significant at the 4σ level is marked with a symbol corresponding to a known explanatory physical model, or a red dot if the signal is not yet explained. The fact that we understand most of the significant signals present in our experiment, combined with the fact that the statistical distribution of the remaining signals below the significance threshold is consistent with a normal distribution, gives us added confidence in our models of the experiment and our reported eEDM result.
Channels/correlations marked with symbols are significantly nonzero due to known mechanisms as follows: • Green stars: Correlations due to the nonzero and drifting signal in the ω N channel described in section 5.5. • Light blue squares: Signals in ω N EB channels due to theB-odd AC stark shift phase coupling to Ω N E r as described in section 5.2.6. • Orange triangles: Correlations due to contrast or asymmetry coupling to Ω N E r . Contrast correlations arise simply because there is a linear dependence of total contrast on Rabi frequency, and the asymmetry correlation is described in section 5.3. • Brown diamond: Correlations in C N and related contrast channels due to nonzero Rabi frequency correlations Ω N r and Ω N P r . These arise due to laser power correlations with theÑ andP switches and due to Stark mixing between rotational levels in H and C, which createÑ -andP-correlated transition amplitudes on the H → C transition as described in section 5.6.2. • Red dot: Signals above our significance threshold for which we have been unable to find a plausible explanation. Even if these quantities arise from real physical effects, they would need to couple to other correlated quantities to contribute to ω N E and there is no evidence for this in the eEDM dataset.

Systematic Error Budget
The method used for construction of a systematic uncertainty varies from experiment to experiment (see for example [123,124]), and it is ultimately a subjective quantity. Even if individual contributions are derived from objective measurements, their inclusion or exclusion in the systematic uncertainty is subjective. Furthermore, the systematic uncertainty cannot possibly be a measure of the uncertainty in all systematic errors in the experiment, but rather only those which were identified and searched for. Although we work hard to identify all significant systematic errors in the measurement, we cannot rule out the possibility that some were missed. Our criteria for including a given quantity in the systematic uncertainty consist of three classes of systematic errors in order of decreasing importance of inclusion: (A) If we measured a nonzero correlation between ω N E and some parameter which had an ideal value in the experiment, we performed auxiliary measurements to evaluate the corresponding systematic error and subtract that error from ω N E to obtain ω N E T . The statistical uncertainty in the shift made to ω N E contributed to the systematic uncertainty.
(B) If we observed a signal in a channel that we deemed important to understand, and it was not understood, but was not observed to be correlated with ω N E , we set an upper limit on the shift in ω N E due to a Figure 35: Over 4,000 switch-parity channels (left) and correlations between switch parity channels (upper right) computed from the eEDM data set. The deviation of each quantity from zero in units of the statistical uncertainty is indicated by the grayscale shading. We set a significance threshold of 4σ above which there is a probability of p = 25% of finding at least 1 false positive. We mark each significant channel/correlation with a symbol corresponding to a model known to produce a signal in that channel. The quantities below this threshold exhibit a normal distribution, shown in the lower right. possible correlation between the two channels. Since such a signal represented a gap in our understanding of the experiment, we added this upper limit as a contribution to the systematic uncertainty.
(C) If a similar experiment saw a nonzero, not understood correlation between their measurement channel and some parameter with an ideal experimental value, but we did not observe an analogous correlation, we set an upper limit on the shift in ω N E due to this imperfection. Since this signal may have signified a gap in our understanding of our experiment, we added this upper limit as a contribution to the systematic uncertainty. Table 4 contains a list of the contributions to our systematic error, grouped by inclusion class, with the corresponding shifts and/or uncertainties. Accounting for class A systematic errors was obligatory, and the removal of these errors from ω N E can be viewed as a redefinition of the measurement channel to ω N E T which does not contain those unwanted effects. These systematic errors consisted of those that depended on the parameters E nr , Ω N E r , and ω E as described in sections 5.2.6 and 5.4, and as such our reported measurement of the T -odd spin precession frequency is defined as ω N E The class B and class C systematic errors were included in the systematic uncertainty to lend credance to our result despite unexplained signals and unexplained systematic errors in experiments similar to ours. All uncertainties in the contributions to the systematic error were added in quadrature to obtain the systematic uncertainty.
With reference to the class B criterion, we deemed the following channels as important to understand: ω N , ω E , ω EB , and ω N EB . Signals were initially not expected in any of these channels and could be measured with the same precision as ω N E . The ω nr , ω B and ω N B channels were not included in our systematic error since the Zeeman spin precession signals present in these channels had non-stationary means and additional noise due to drift in the molecule beam velocity. Only one of these channels, ω N , described in section 5.5, met the class B inclusion criterion.
With reference to the class C criterion, we defined the set of experiments similar to ours to include other eEDM experiments performed in molecules: the YbF experiment [42] and the PbO experiment [50]. The PbO experiment observed unexplained systematic errors coupling to stray magnetic fields and magnetic field gradients (cf. section 5.7), and the YbF experiment observed unexplained systematic errors proportional to detunings (cf. section 5.6.1) and a field plate ground voltage offset (cf. section 5.8). Thus we included the systematic uncertainty associated with the aforementioned effects in our budget.
After having accounted for the systematic errors and systematic uncertainty, we reported ω N E T , the where the combined uncertainty is defined as the quadrature sum of the statistical and systematic uncertainties, σ 2 = σ 2 stat + σ 2 syst . This result is consistent with zero within 1σ. Since σ syst is to some extent a subjective quantity, its inclusion should be borne in mind when interpreting confidence intervals based on σ. Nevertheless, this inclusion decision does not have a large impact on the meaning of the resulting confidence intervals since σ is only about 20% larger than σ stat .

Confidence Intervals
A classical (i.e. frequentist) confidence interval [125] is a natural choice for reporting the result of an eEDM measurement. For repeated and possibly different experiments measuring the eEDM, the frequency with which the confidence intervals include or exclude the value d e = 0 suggests whether the results are consistent or inconsistent, respectively, with the Standard Model. Furthermore, the confidence level (C.L.) represents an objective measure of the a priori probability that the confidence interval assigned to any one of these measurements, selected at random, includes the unknown true value of the eEDM d e,true . Since no statistically significant eEDM has yet been observed, the recent custom has been for electron eEDM experiments to report an upper limit at the 90% C.L. [41,42]. The proper interpretation of such limits is that if the experiment were performed a large number of times, and the confidence interval were computed in the same way for each experimental trial, d e,true would fall within the interval 90% of the time.
Feldman and Cousins pointed out that in order for this interpretation to be valid, the confidence interval construction must be independent of the result of the measurement [126]. If the procedure for constructing 90% confidence intervals is chosen contingent upon the measurement outcome, the resulting intervals may 'undercover', i.e. fail to include the true value more than 10% of the time. This happens, for example, if an upper bound is reported whenever the measured result falls within a few standard deviations of zero, and a two-sided confidence interval is reported whenever the measured result is significant at more than a few-sigma level. Feldman and Cousins termed this inconsistent approach 'flip-flopping'.
In order to avoid flip-flopping, we chose a confidence interval construction, the Feldman-Cousins method described in reference [126], that consistently unifies these two limits. We applied this method to a model with Gaussian statistics, in which the measured magnitude of the eEDM channel, x = |ω N E T,meas |, is sampled from a folded Gaussian distribution where the location parameter is the unknown true magnitude of the eEDM channel, µ = |ω N E T,true |, and the scale parameter σ is equal to the quadrature sum of the statistical and systematic uncertainties given in equation 94 and at the bottom of table 4.
The central idea of the Feldman-Cousins approach is to use an ordering principle which, for each possible value of the parameter of interest µ, ranks each possible measurement outcome x by the 'strength' of the evidence it provides that µ is the true value. The values of x that provide the strongest evidence for each value of µ are included in the confidence band for that value. In the Feldman-Cousins method, the metric for the strength of evidence is the likelihood of µ given that x is measured [i.e. L(µ|x) = P (x|µ)], divided by the largest probability x can possibly achieve for any value of µ. The denominator in this prescription takes into account the fact that an experimental result that is somewhat improbable under a particular hypothesis can still provide good evidence for that hypothesis if the result is similarly improbable under even the most favorable hypothesis. This approach has its theoretical roots in likelihood ratio testing [127].
Our specific procedure for computing confidence intervals was a numerical calculation performed using the following recipe (cf. figure 36): (i) Construct the confidence bands on a Cartesian plane, of which the horizontal axis represents the possible values of x and the vertical axis the possible values of µ. Divide the plane into a fine grid with x-intervals of width ∆ x and µ-intervals of height ∆ µ . We will consider only the discrete possible values x i = i∆ x and µ j = j∆ µ , where the index i(j) runs from 0 to n x (n µ ).
(ii) For all values of i, maximize P (x i |µ j ) with respect to µ j . Label the maximum points µ max,i .
(iii) For some value of j, say j = 0, compute the likelihood ratio R(x i ) = P (x i |µ j )/P (x i |µ max,i ) for every value of i.
(iv) Construct the 'horizontal acceptance band' at µ j by including values of x i in descending order of R(x i ). Stop adding values when the cumulative probability reaches the desired C.L. of 90%, i.e., xi P (x i |µ j )∆ x = 0.9.
(v) Repeat steps (iii)-(iv) for all values of j.
(vi) To determine the reported confidence interval, draw a vertical line on the plot at x = |ω N E T,meas |. The 90% confidence interval is the region where the line intersects the constructed confidence band.  The left-hand plot in figure 36 was generated using the prescription above at several different C.L.'s. Note that the 90% confidence intervals switch from upper bounds to two-sided confidence intervals when the value of |ω N E T,meas | becomes larger than 1.64σ. This is the level of statistical significance required to exclude the value d e = 0 from a 90% C.L. central Gaussian confidence band.
From equation (94), we find |ω N E T,meas | = 0.46σ with σ = 5.79 mrad/s. In our confidence interval construction, this corresponds to an upper bound of |ω N E T | < 1.9σ = 11 mrad/s (90% C.L.). A comparison between three different 90% confidence interval constructions for small values of µ is shown in the right-hand plot of figure 36. The black dashed lines represent the central confidence band for the signed values (rather than the magnitude) of µ and x, where µ is the mean of a Gaussian probability distribution in x. The blue lines give an upper bound constructed by computing the the value of µ such that the cumulative distribution function for the folded Gaussian in equation 95 is equal to 0.9 for each value of x. It should be noted that this upper bound is more conservative than a true classical 90% confidence band, as it overcovers for small values of µ (e.g., if the true value were µ true < 1.64σ, the confidence intervals of 100% of experimental results would include µ true ). We nevertheless include this construction for comparison because we believe that previous experiments have reported EDM upper bounds using this method [41,42,106]. These intervals have a valid interpretation as Bayesian 'credible intervals' conditioned on a uniform prior for µ [126]. Finally, the red lines represent the Feldman-Cousins approach described here, which unifies upper limits and two-sided intervals. For our measurement outcome, indicated by the vertical dot-dashed line, the Feldman-Cousins intervals yield a 7% larger eEDM limit than the folded Gaussian upper bound would have.

Physical Quantities
Under the most general interpretation, our experiment is sensitive to any P -and T -violating interaction that produces an energy shift ω N E T . The eEDM is not the only such predicted interaction for diatomic molecules [128], and in particular a P -and T -odd nucleon-electron scalar-pseudoscalar interaction would also manifest as aÑẼ-odd phase in our experiment. Thus, we write where W S is a (calculated) energy scale specific to the species of study [28,29,[129][130][131] and C S is a dimensionless constant characterizing the strength of the T -violating nucleon-electron scalar-pseudoscalar coupling relative to the ordinary weak interaction. We can use our measurement to set an upper limit on d e by assuming that C S = 0 and that ω N E T is therefore entirely attributable to the eEDM. Taking the effective electric field to be the unweighted mean of the two most recent calculations of this quantity [28,29], E eff = 78 GV/cm, we can interpret our result in equation (94) as: ⇒ |d e | < 9.3 × 10 −29 e · cm (90% C.L.), where the second line is obtained by appropriately scaling the upper bound on ω N E T derived in section 6.1. If, instead, we assume that d e = 0, our measurement of ω N E T in ThO can be restated as a measurement of C S . Using an unweighted mean of the most recent calculations of the interaction coefficient, W S = −2π × 282 kHz [28,29], we obtain: ⇒ |C S | < 6.2 × 10 −9 (90% C.L.), which, at the time, was an order of magnitude smaller than the existing best limit set by the 199 Hg EDM experiment [114], and is still a factor of 2 smaller than the recently improved limit from the same group [132].

Summary and Outlook
Our new limit on the size of the electron's electric dipole moment [27]: d e ≤ 9.3 × 10 −29 e · cm. (101) † Note that the sign of the C S term is opposite to that used, incorrectly, in our original paper [27]. In addition, here W S differs in magnitude from the related quantity W T,P given explictly in [28,29]. A detailed discussion of the sign and notational conventions for this Hamiltonian is provided in Appendix A.
represented an order of magnitude improvement on previous bounds [41,42] and more strongly constrained the viable parameter space for many extensions to the Standard Model, while probing one-loop effects of new physics at a mass scale of ∼10 TeV.
We have presented our experimental method for measuring an eEDM-induced precession phase in the dipolar molecule ThO, detailing the way we utilise several experimental switches to isolate the component of accumulated phase with the correct symmetry properties. We described the apparatus that we used to carry out our measurement and have presented a thorough analysis of the systematic errors present in the experiment, showing in detail the approach to finding and quantifying shifts of the eEDM-associated phase and their corresponding uncertainties.
Despite the success of the experiment in reducing the limit on the value of the eEDM, there are several aspects of the experimental procedure that we are improving on which will significantly further enhance our statistical sensitivity. These upgrades include: • Thermochemical Source: Instead of relying on ablation to generate ThO molecules from a ThO 2 target, a relatively uncontrolled process, we are developing a new method using a thermochemical reaction-based beam source. This relies on the specific reaction [133,134] Th(s) + ThO 2 (s) → 2ThO(g) (102) occurring in a precursor target made of a Th/ThO 2 mixture. Preliminary tests have demonstrated a roughly factor of 10 increase in the time-averaged molecular flux produced via this method. • Beam Geometry: In the current experiment, the molecules in the spin-precession region subtend a solid angle of ∼60 µsr relative to the beam source, meaning only ∼10 −5 -10 −4 of molecules produced reach the state readout region. This useful fraction of molecules can be increased in two ways: by shortening the distance between beam source and spin-precession region, and by increasing the spacing between the electric field plates so as to accomodate a beam with a larger transverse size. By making both of these changes to the apparatus we can increase the usable molecule number by a factor of ≈8. • State Preparation: In the current experiment we transfer molecules into the H state by optically pumping via |A, J = 0 (see section 3.2.4). This procedure is inefficient; only ∼35% of molecules addressed by the excitation laser are transferred into the H-state manifold, within which 1/6 of the population is in the desired superposition state |B(ˆ prep ,Ñ ,P) . We can significantly increase the number of molecules prepared by using stimulated Raman adiabatic passage (STIRAP) to perform coherent population transfer from X to H via C. We have demonstrated an estimated efficiency of 75% which will increase the usable molecule number by a factor of ≈12 [135]. • State Readout: We will be changing the transition which we perform our state readout on. The I state of ThO is another Ω = 1 state which has a number of advantages over the C state, namely a ∼10 times higher transition dipole moment from the H state, a larger branching ratio to the X state and a shorter fluorescence wavelength to X [65,92,136]. The latter allows for higher quantum efficiency detection of photons. In addition, the efficiency with which we collect the light is improved by using light pipes instead of fiber bundles. Together we anticipate a factor of ≈6 improvement in signal [135]. • Other Improvements: The suppression of known systematic effects was limited only by statistics. To the best of our knowledge, the limit on ∂ω N E /∂E nr (see section 5.2) could have been 10 times smaller if we had collected the data required to tune out that slope with such precision. Therefore, there is no reason to believe that the systematic effects we have discovered in this first generation measurement will limit the next generation of the experiment. However, we are taking additional measures to suppress such systematics, such as new electric field plates designed to minimise the absorbed laser power and hence the birefringence.
The ACME search for the electric dipole moment of the electron is now entering its second generation and we anticipate a new measurement that will either find a nonzero value of d e , or constrain it to be 10 −29 e · cm, thus probing one-loop interactions at an energy scale of ∼30 TeV. again S = 1 to a fair approximation). This choice appears, at first glance, to contradict the discussion in section 2.1, where for a single electron we wrote d e = 2d e s (where s = 1/2). However, these two definitions are in fact consistent when taking into account that in the H 3 ∆ 1 state of ThO only one of the two valence electrons (the one in the σ orbital) contributes significantly to the EDM energy shift, while both electrons contribute to the total spin S = 1. Hence, in our formulation, the molecule-frame projection d eff e ·n can take extreme values ±d e , as expected for a single contributing electron. (This 'single contibuting electron' approximation is valid for all molecules used to date in searches for the eEDM.) We then write the effective eEDM Hamiltonian H eff EDM in the standard form for interaction of an electric dipole moment with the internal electric field: where the + sign in the final expression arises from the sign convention for E int . In eigenstates of Ω, the expectation value of H eff EDM -that is, the energy shift ∆E EDM due to the eEDM-can be written as Now, we finally re-introduce the effective electric field E eff used throughout the main text of this paper. This is related to the internal field introduced above, via We can then use this notation to describe the effective nonrelativistic eEDM interaction, within a given electronic state and eigenstate of Ω (and otherwise independent of molecular structure), as follows: where the sign in the last expressions arises from the defined definitions of Σ (component of S alongn) and E eff (antiparallel ton). All relevant quantities are summarised pictorially in figure A1. In most of the theoretical literature on this subject, this energy shift is written in the unambiguous form ∆E EDM = +d e W d Ω. However, there has been no consistent definition in the literature for the relation between W d and E eff . In particular, both their relative signs and the dependence of their relative magnitude on the value of |Ω| (encompassing both the case of one-and two-electron systems) are often defined differently, or imprecisely. In our notation, the expressions above imply a general relationship between E eff and W d : This relation is valid for systems with one or two valence electrons (in the 'single contributing electron' approximation for the latter case), and regardless of the relative directions of Σ and Ω. Now we apply these general considerations to the specific case of the H state of ThO. Here, since Σ ≈ −Ω, we find that E eff = −W d with our conventions. Thus, the energy shifts can be written for ThO as ∆E EDM = −d e E eff Ω. (A.10) In our experiment, this gives rise to energy shifts, for a given direction of the laboratory electric field E, given where G F is the Fermi coupling constant, γ i are Dirac matrices, ρ N ( r) is the normalised nuclear density, Z(N ) is the proton (neutron) number, and C S,p and C S,n are dimensionless constants which describe the interaction strength (relative to that of the ordinary weak interaction) specifically for protons and neutrons, respectively. Using the definition C S = Z A C S,p + N A C S,n = Z A C S,p + 1 − Z A C S,n , (A.14) where A = Z + N , C S represents a weighted average of the couplings to protons and neutrons, and is different for every nuclear species. However, since the ratio Z/A is nearly the same for all heavy nuclei used in molecular and atomic EDM experiments (ranging only from Z/A = 0.41 for 133 Cs to Z/A = 0.39 for 232 Th), typically a common value for C S is assumed for all experiments of this type. Thus we can write In a given molecular electronic state, this gives rise to a non-relativistic, single-electron effective Hamiltonian of the form H eff SP = 2 s ·nC S Y S ; the factor of 2 is included so that the maximal energy shifts due to this term have the simple form ∆E max SP = ±C S Y S . By analogy with our discussion of the eEDM Hamiltonian, in a molecular state with S = 1 and a 'single contributing electron', as in the H 3 ∆ 1 state of ThO, we rewrite this in the form Here, in our notation, W S ≡ Y S [ Σ /Ω] (≈ −Y S in ThO). However, quantities analogous to Y S (in terms of which the energy shifts depend explicitly on the spin direction) are rarely introduced in the literature; instead, only forms analogous to W S (where the energies depend only on Ω) are used.
Our definition for C S was historically a standard notation used in the literature. However, in some recent papers (e.g. references [28,138]) it is implicitly assumed that the neutron coupling C S,n vanishes. In these papers, the factor AC S in equation A.15 is replaced by ZC S,p (or its equivalent in a different notation), † and the energy shift is written in the analogous form ∆E SP = C S,p W S,p Ω. These papers report values of W S,p in the H state of ThO, based on sophisticated calculations of the molecular wavefunctions. However, since there is no particular reason to expect this interaction to couple more strongly to protons than to neutrons, we prefer to report our results in terms of C S . To do so, we use the relation W S = (A/Z)W S,p to determine W S from the reported values for W S,p .
Finally, the experimentally determined energy shift arising from the nucleon-electron SP interaction is 19) and the total T-violating energy shift is Note that the sign of the C S term is opposite to that used, incorrectly, in our original paper [27]. † In reference [138] our C S,p is denoted as k T,P and our W S,p as W T,P ; in Ref. [28], our W S,p is denoted simply as W S . References [131,139] denote our C S as C SP and our W S as Wc. Table A1: Summary of the different conventions used in some of the literature relating to eEDM measurements/theory. Where entries are left blank the convention is not stated in the reference provided. Quantities in square brackets are not explicitly stated in the references but are implied. In some cases, nomenclature has been modified for consistency. Footnotes provide specific references for the equations shown.

Appendix B. Glossary of Abbreviations and Symbols
Appendix B.

Experiment Switches
During the course of the experiment, we performed many parameter switches. Most of these switch parameter symbols are denoted by a superscript tildeX , which indicates that that parameter takes on two values, X = ±1.
N Used as a quantum number,Ñ ≈Ẽsgn (MΩ), for states of |M | > 0, |Ω| > 0, that refers to states with opposite molecular alignment with respect to the applied electric field. It is also used to refer to the experiment switch between spectroscopically addressing states in |H, J = 1 with opposite values ofÑ .
E Denotes the alignment of the applied electric field with respect to the laboratoryẑ axis,Ẽ = sgn Ẽ ·ẑ where E is the applied electric field.
B Denotes the alignment of the applied magnetic field with respect to the laboratoryẑ axis,B = sgn B ·ẑ where B is the applied magnetic field. θ Denotes the state of the polarisation dither that is used to extract the contrast in the spin precession measurement. It refers to the direction of the offset angle in the xy plane of the state readout polarisation basisX,Ŷ , relative to the average polarisation of these lasers. P Used as a quantum number to denote the parity (eigenvalue of the parity operator P ) of a given molecular state of well-defined parity. It is also used to refer to the experiment switch between spectroscopically addressing states in |C, J = 1 with opposite values ofP with the state readout lasers. L Denotes the state of the mapping between the two output channels of the electric field voltage supply, and the two electric field plates which can be either connected normally (+1), or inverted relative to normal (-1). R Denotes the state of an experimental switch of the state readout polarisation basis offset angle with respect to the x-axis by either 0 (+1) or π/2 (−1). G Denotes the state of an experimental switch of the global polarisation; the state preparation and state readout lasers are rotated synchronously by a common angle. This can be thought of as a redefinition of thex andŷ axes in the xy plane. B z Denotes the magnitude of the magnetic field along theẑ direction in the laboratory, B z = | B ·ẑ|. This parameter is switched between three values differing by about 20 mG. In figure 35, channels X that are 'odd' with respect to this parameter refer to the linear variation ∂X/∂B z .
E Denotes the magnitude of the electric field, E = | E|. This parameter is switched between two values.
k ·ẑ Denotes the orientation of both the state preparation and the state readout laser pointing directions with respect to the laboratoryẑ axis. This is a binary switch,k ·ẑ = ±1, but we do not denote this switch with a tilde as we do with the other binary switch parameters.

Appendix B.2. Laser Parameters
There are a variety of laser parameters which are used to describe the state preparation laser that is denoted with a subscript 'prep', or the state readout lasers that are denoted with a subscript 'read' if the property applies to both state readout lasers, or with subscripts X and Y , if the parameter can vary between the two readout lasers.
k Laser pointing direction. In this paper, the pointing direction is always nearly aligned or antialigned with respect to the laboratoryẑ axis such thatk ·ẑ ≈ ±1. ϑ k Defined in equation 54. Polar angle of deviation of the pointingk from aligned or anti-aligned with thê z axis. ϕ k Defined in equation 54. Azimuthal angle denoting the direction in the xy plane, relative to the x-axis, of the deviation of the pointingk from theẑ axis. Complex laser polarisation. The readout laser polarisations are also referred to asX andŶ as an alternative toˆ X andˆ Y at some points.
ε Effective polarisation. Used to parameterize the effect of experiment imperfections on the molecule state as the polarisation vector that would be required to obtain the same molecule state in the absence of those experiment imperfections. Ω r Rabi frequency for a particular laser beam and transition. Defined as the transition dipole matrix element multiplied by the amplitude of the electric field associated with the laser beam.
Γ Optical retardance for some birefrigent element along the laser beam path.
φ Γ Angle in the xy plane of the fast axis associated with an optical retardance Γ.

Appendix B.3. Molecular States and Parameters
These symbols are all used to describe the molecular energy level structure and the manner in which our laser light interacts with the molecules, in particular for the state preparation and readout processes.
J Total angular momentum.
M Projection of J onto the laboratoryẑ-axis.
Ω Projection of J onto the internuclear axis,n.
B H Rotational constant of the H state. E eff 'Effective electric field' to which we consider the eEDM to be subjected.
∆ Ω, 1 The Ω-doublet splitting of the |H, J = 1 state. γ Decay rate of the a given electronic state. The electronic state label is given in the subscript. In most of the paper, only γ C , the decay rate of the C state, is relevant.
Ω r Transition Rabi frequency, which is proportional to the square root of the laser intensity.
E B± , E D Instantaneous eigenenergies of the dressed three-level system, defined in equation 64.
χ Complex polarisation rotation rate defined in section 5.2.3.
Π Defined and discussed in section 5.2.3 and equation 74. This is a factor in the AC Stark shift phase that is independent of laser polarisation but depends on the laser detuning and Rabi frequency.
v The mean longitudinal velocity of the molecular beam.

Appendix B.4. Measurement Quantities
These symbols represent quantities related to the measurement of the accumulated phase and the way in which it is extracted during data analysis, as well as some related quantities pertaining to systematic studies.
N Total number of measurments performed, equivalent to the number of detected photoelectrons.
N 0 Number of molecules in the state readout region in the particularÑ level being addressed.
f Fraction of fluorescence photons emitted in the state readout region that are detected.
S Recorded photoelectron count rate measured on the photodetectors.
F Photoelectron count rate due to the molecule fluorescence. F X,Y is used to denote the molecular fluorescence induced by the X and Y state readout lasers, respectively. F cut is used to denote the fluorescence threshold above which data was included in the analysis.
B Background count rate primarily due to scattered light from the state readout lasers. This background signal is subtracted from the raw photoelectron signal S to obtain the fluorescence photoelectron count rate, F = S − B. A Signal asymmetry as defined in equation 20.
C Spin precession fringe contrast, as defined in equation 21, is the sensitivity of the asymmetry to molecular spin precession.
φ Actual spin precession phase of the molecules as defined in equation 16.
τ Measured spin precession time as described in sections 3.1 and 4.2.2. ω Measured spin precession frequency, as defined in equation 37, ω = Φ/τ . ω N E The measurement channel of interest, the spin precession frequency channel that is correlated withÑ andẼ. The expected eEDM signal should contribute to this channel.
ω N E

T
The contribution to spin precession frequency ω N E induced by T -odd spin precession effects in the H state in ThO.
ω N E P A systematic error in the ω N E channel that is proportional to some parameter P .