Revisiting $^{129}$Xe electric dipole moment measurements applying a new global phase fitting approach

By measuring the nuclear magnetic spin precession frequencies of polarized $^{129}$Xe and $^{3}$He, a new upper limit on the $^{129}$Xe atomic electric dipole moment (EDM) $ d_\mathrm{A} (^{129}\mathrm{Xe})$ was reported in Phys. Rev. Lett. 123, 143003 (2019). Here, we propose a new evaluation method based on global phase fitting (GPF) for analyzing the continuous phase development of the $^{3}$He-$^{129}$Xe comagnetometer signal. The Cramer-Rao Lower Bound on the $^{129}$Xe EDM for the GPF method is theoretically derived and shows the potential benefit of our new approach. The robustness of the GPF method is verified with Monte-Carlo studies. By optimizing the analysis parameters and adding data that could not be analyzed with the former method, we obtain a result of $d_\mathrm{A} (^{129}\mathrm{Xe}) = 1.1 \pm 3.6~\mathrm{(stat)} \pm 2.0~\mathrm{(syst)} \times 10^{-28}~ e~\mathrm{cm}$ in an unblinded analysis. For the systematic uncertainty analyses, we adopted all methods from the aforementioned PRL publication except the comagnetometer phase drift, which can be omitted using the GPF method. The updated null result can be interpreted as a new upper limit of $| d_\mathrm{A} (^{129}\mathrm{Xe}) |<8.3 \times 10^{-28}~e~\mathrm{cm}$ at the 95\% C.L.


Introduction
A quantum field theory that models the formation of the imbalance of matter over antimatter in our universe must fulfill the Sakharov conditions [1]. One of those conditions is the CP violation (C is charge conjugation and P is parity reversal). The best-tested standard model (SM) of particle physics provides two sources of CP violation, the phase of the Cabibbo-Kobayashi-Maskawa matrix and the termθ in the QCD Lagrangian [2]. However, the CP violation within the SM is too small to produce the observed rate of the matter to antimatter asymmetry, motivating searches for physics beyond-the-SM (BSM). BSM theories generally include additional sources of CP violation [2,3], such as a larger permanent electric dipole moment (EDM) of fundamental or composite particles [4,5]. So far, all measurement results of EDMs in more than ten diverse systems, with the first published in 1957 [6], are consistent with zero. These null results are interpreted as upper limits on EDMs and place constraints on various sources of CP violation and masses of BSM particles, thus directing the search of BSM scenarios [7].
Long spin-coherence time and obtainable high polarization leading to high signal-to-noise ratios (SNR) make several diamagnetic systems such as the 199 Hg and 129 Xe atom promising candidates for EDM experiments. Over the last 40 years, significant progress was made in the determination of upper limits for EDMs of diamagnetic systems (see Fig. 1). At present, the 199 Hg atomic EDM measurement is the most sensitive, and its upper limit sets constraints to multiple sources of CP violation [8]. Considering various potential contributions to an atomic EDM, an improved limit on other systems, like the 129 Xe EDM d A ( 129 Xe), will tighten these constraints. The theoretical results for 129 Xe EDM are more accurate and reliable than those obtained for 199 Hg EDM, therefore 129 Xe has the potential to probe new physics [9].
Recently, new upper bounds on the 129 Xe EDM using 3 He comagnetometry and SQUID detection have been reported by a joint collaboration between the University of Michigan, the Technical University of Munich and the Physikalisch-Technische Bundesanstalt(PTB) [11] as well as another independent group with comparable sensitivities [15], which are about five times smaller than the previous limit set in 2001 [16]. One of the challenges in both L. For all systems, the current upper bound has decreased more than an order of magnitude compared to their first published result [8,10,11,12,13,14,15,16].
experiments is the comagnetometer frequency drift, which is several magnitudes larger than the expected frequency shift due to a potential 129 Xe EDM [17]. One approach to correct for the impact of the comagnetometer drift on the measured d A ( 129 Xe) is using a deterministic physical model to fit the comagnetometer frequency drift [15,18]. However, the physical origin of the comagnetometer frequency instability is subject of a controversial debate [19,20], which was inspired by another theoretical model and motivated the performance of recent experiments to substantiate the former criticism [21,22]. Instead, in Ref. [11] a phenomenological method was used, which does not need any physical model on the comagnetometer frequency drift, but a distinct pattern of electric fields with switching polarity.
We will refer to that as the Pattern Combination (PC) method from here on.
Here, we propose a new analysis based on a Global Phase Fitting (GPF) method, where the EDM value is estimated by a single fit to the comagnetometer phase development within one complete measurement. Besides an experimentally deduced EDM function as used in Ref. [15], allowing to analyse any electric field pattern, our GPF method uses a polynomial function to account for the comagnetometer frequency drift. Sec. 2 gives a short description of the basic principle of measuring the 129 Xe EDM d A ( 129 Xe) using comagnetometry. In addition, the PC method is introduced for comparison with the GPF method. The GPF method is elucidated in detail in Sec. 3, including the derivation of the Cramer-Rao Lower Bound (CRLB). The CRLB of the variance on the EDM value estimation using the GPF method is a factor of four smaller than that of the PC method. In Sec. 4 we validate the GPF method with Monte-Carlo simulations and compare the results of the PC and GPF method using the experimental data obtained for Ref. [11]. Eventually we recalculate the systematic uncertainties based on Ref. [11] and derive a new upper limit for the permanent 129 Xe EDM. For 129 Xe atoms stored in a cell permeated by a uniform magnetic field B and an electric field E, that is parallel to B, their nuclear spin precesses at an angular frequency where F Xe = 1/2 is the total angular momentum number and γ Xe is the gyromagnetic ratio of 129 Xe. The magnetic field B in Eq. (1) becomes an interference term for directly calculating d A ( 129 Xe) from ω Xe . To overcome the experimental difficulties on controlling and measuring B, comagnetometry was introduced with two collated species measured at the same time [16,17,23]. 3 He is an ideal candidate for comagnetometry due to its potentially high SNR and a negligible EDM compared to d A ( 129 Xe) [24]. The weighted frequency difference between 129 Xe atoms and 3 He atoms is defined as and commonly named the comagnetometer frequency. Here ω He = |γ He B| is the spin precession frequency of 3 He atoms with γ He being its gyromagnetic ratio. Therefore, ω co can be written as showing that ω co is independent of the magnitude of the background magnetic field but depends on its orientation relative to the applied electric field. The current measurement sensitivity of ω co is in the nHz range for a single measurement, while the comagnetometer frequency drift is at the µHz level, which causes a non-negligible systematic error [21,22,25]. Multiple physical models to describe the comagnetometer drift were proposed. The dominant terms thereby vary in different models. Furthermore, several parameters, such as the longitudinal relaxation time T 1 of the nuclear spins, used in these models are unknown or difficult to measure, making the frequency drift correction with a deterministic model inaccurate. By using a phenomenological model such as proposed here and in [11], these currently unsolved difficulties can be omitted.

Parameters of two measurement campaigns
The data used in our analysis were collected in the joint collaboration at the Berlin Magnetically Shielded Room (BMSR-2) facility at PTB Berlin. Table 1 summarizes the main experimental parameters of the two measurement campaigns carried out in 2017 and 2018, respectively. More details on the setup and process are given in Ref. [26]. The spin precession signal of the transverse magnetization of 3 He and 129 Xe was recorded by a dc-SQUID system with two channels (Z1,Z2). The high voltage and leakage current between the two electrodes of the cell were monitored. A background magnetic field B 0 in the range of 2.6 µT -3 µT was applied to shift ω Xe and ω He to 30 Hz -36 Hz and 90 Hz -98 Hz, respectively, which are well above the vibrational interference signals (see Fig. 2). In order to further decrease the impact of the vibrational noise, a software SQUID gradiometer (Z1 − Z2) was used. The left panel of Fig. 2 shows the raw SQUID gradiometer signal in pT (gray) of one run from the 2018 campaign lasting 35000 s exemplarily. This run comprises two so called sub-runs with 36 segments each. A segment is defined as the time of constant electric field. For the two sub-runs shown in Fig. 2, the segments last 300 s and 600 s, respectively. The first sub-run ranging from 50 s to 12400 s is used as an example in the data analysis section.

PC method
As mentioned above, one approach to mitigate the effect of the comagnetometer frequency drift is repetitively reversing the direction of the electric field E. This allows to separate the impact of d A ( 129 Xe) on ω co from other interference terms. The E modulation method has been applied in diverse EDM experiments with varied modulation patterns [10,14]. For the PC method, the common E pattern for one sub-run consists of 36 segments with an equal time interval t s , and the sign of E changes according to the following sequence ±[0 + --+ -+ + --+ + -+ --+ 0, 0 -+ + -+ --+ + --+ -+ + -0]. The segments of zero voltage were added to allow for systematic error studies. The PC method determines the EDM value from averaging the comagnetometer frequencies ω co from 2 n (n ∈ N) consecutive segments omitting those with zero voltages. This pattern is constructed to cancel the effect of the comagnetometer frequency drift up to n − 1 order when parametrized in polynomials. The effect of the higher order (above n − 1) drift dependency imposing a false EDM on each sub-run is deduced by applying polynomial fits to all ω co within the sub-runs, leading to a correction for the EDM and an additional systematic uncertainty (for more details see Ref. [26]).

GLOBAL PHASE FITTING METHOD
The general data-processing procedure for the GPF method is illustrated in Fig. 3. For this method, the raw SQUID data of a sub-run is cut into continuous blocks of equal length. Each block data is fitted to deduce precession phases of both species 3 He and 129 Xe (see Sec. 3.1) and the continuous comagnetometer phase is derived for each block (see Sec. 3.2). For data blinding an additional phase, bound to the measured high voltage signal, can be added to the comagnetometer phase at this point (see Sec. 3.3). The EDM value is acquired by fitting the blinded comagnetometer phases using a polynomial function together with a constructed function comprising the phase evolution introduced by a hypothetical 129 Xe EDM. The unblinded EDM result is obtained by reanalyzing the raw comagnetometer phases, as illustrated in Fig. 3.

The phase of each block
The block length t b is a free parameter with a suitable range from 1 s to 20 s, being short enough to exclude the amplitude decay and frequency drift, and long enough to perform the fit for our data [11]. The SQUID data in each block are fitted to the function where a Xe/He/i , b Xe/He/i , ω Xe/He , c, and d are the fit parameters and ω i=1,2,3,4 = 2π × 50i s −1 represent the power frequency and its harmonics. The constant and linear terms c and d · t describe the background magnetic field and its small drift as seen by the SQUID. The variable projection (VP) method is applied [27], where the nonlinear parameters ω Xe/He are estimated separately from the linear parameters a Xe/He/i , b Xe/He/i , c, and d. To minimize the correlation between the fit terms in Eq. (4), the time of each block is assigned to be symmetrical around zero from −t b /2 to t b /2.  shows the raw SQUID data of a 5 s block from the start of the exemplary sub-run and the residual of the fit to this data. The residual is dominated by the mechanical vibration in the frequency range of 4 Hz -25 Hz as shown in the right plot of Fig. 2. We can assume approximate orthogonality between the precession signal and the vibrational noise of our setup. Therefore, the error on the fit parameter values caused by the latter one is negligible compared to that caused by the white noise, although its integrated power is much larger than the white noise power. This was validated with Monte-Carlo simulations using the recorded vibrational noise (see Appendix A.1). The phase of each species for the block k in the range of [−π, π) can be obtained by where Arg is the function to get the principle argument of a complex number, i is the imaginary unit and m = Xe or He. Note that due to the time centering, the estimated phase φ k is referred to the middle time of each block The time interval of the block k is defined as (t k−1 , t k ). The parameter uncertainties δa k m and δb k m are estimated from the covariance matrix of the fit where r is the residual, ν is the degrees of freedom and J is the Jacobian matrix. The standard deviation of the derived phase δφ k Xe/He is Eq. (6) assumes that the residual r stems from the wideband white noise, which is a conservative approach for our case since the main signal in the residuals is the narrowband vibrational noise, leading to an overestimation of the uncertainty δφ k m . However, the ratio between δφ k m for different blocks reflects the decaying SNR. Therefore, these estimated uncertainties are used as weights in the subsequent GPF routine.

The accumulated comagnetometer phase
The accumulated phase Φ k m in a block k of the continuously precessing spins is the sum of the wrapped phase φ k m and a multiple of 2π where the integer n k m is determined as rounded to the lower integer and n 1 m = 0. Here, the frequencies ω k−1 m are obtained by the fit of the block k using Eq.
is either > π or < −π, n k m is incremented or decremented by one, respectively, to ensure a continuous phase evaluation. The standard deviation of the accumulated phase δΦ k m is equal to δφ k m as Eq. (8) does not introduce any additional uncertainty. According to Eq. (2), the evolved comagnetometer phase Φ k co for each block k is determined by

The fitted EDM value
By integrating Eq. (3), the accumulated phase due to a hypothetical 129 Xe EDM d h at the block k is where E i is the average electric field within the block i. By replacing d h with a computer-generated pseudo-random EDM value d bias , the bias phase Φ k bias is calculated and then used to blind the comagnetometer phase Φ k co,b = Φ k co +Φ k bias in order to avoid operator induced bias during process optimization. The value of d bias has been saved in an independent file in a binary format and Φ k co,b was used for later data analysis. The measured phase Φ k co originates not only from the potential 129 Xe EDM, but also from other sources such as chemical shift [21,26]. These contributions are phenomenologically parametrized by a polynomial of gth order [28]. Hence, the comagnetometer phase is fitted with the function where a, p 0 , p 1 , p 2 , . . . , p g are the global fit parameters. Here the time series t k are normalized to the interval [0,1] and shifted Legendre polynomialsP n (t k ) are applied to decrease the correlation between polynomial coefficients [29]. The fit was conducted by using the iterative least squares estimation method with the built-in function nlinf it in MATLAB. Thereby the inverse values of the phase variances (δΦ k co ) 2 are used as weights. Fig. 5 shows the comagnetometer phase Φ k co , the fit phase Φ k fit , and the EDM function Φ k EDM constructed from the measured E-field pattern of the exemplary sub-run. To determine the order needed for the polynomial function in Eq. (12), we apply an F -test where the significance of adding q terms to the fitting function with g terms was evaluated by the integral probability where P F is the probability density function of the F -distribution and N is the number of data points [30]. The upper bound of the integral is The order of the fit was defined sufficient when P g,g+1 as well as P g,g+2 are both smaller than a chosen threshold of P min . The atomic EDM of 129 Xe is calculated from the fit parameter a as The correlated uncertainties of the parameters are determined as the square root of the reciprocal of the diagonal of the covariance matrix, which inherently includes the uncertainty of the correlations between a and polynomial parameters. The influences of these correlations to the estimation of a are small due to the orthogonality between the constructed function Φ k EDM and the polynomial function of the order up to n − 2 where 2 n is the number of nonzero high voltage segments. The correlation matrix for the exemplary sub-run (see Fig. 2) is given in Table 2. In this case, the correlations between the EDM parameter a and the polynomial coefficients are significantly smaller than 1, but nonzero, since the polynomials of higher than 3rd order are not orthogonal to Φ k EDM . The derived uncertainty is in good agreement with the result using the log profile likelihood method. We also applied the linear regression method with the model in Eq. (12) and obtained consistent results.

The modified Allan deviation
The modified Allan deviation (MAD) is an established tool to evaluate the low-frequency drift of a time series of phases Φ, which is defined as where the integration time τ is n times the block length t b , and the total measurement time T is subdivided into P time intervals of equal length τ , such that P τ ≈ T [31]. As an example, the MAD of the exemplary sub-run is plotted in Fig. 6. σ f of Φ k co reaches the minimum at the integration time τ of 550 s and then increases due to the comagnetometer frequency drift. For the residual phase Φ k co − Φ k fit of this exemplary sub-run, the MAD decreases with increasing integration time according to σ f ∝ τ −3/2 (dashed line in Fig. 6) over the considered range, down to 0.4 nHz. This behavior is an indicator that the comagnetometer phase Φ k co is adequately described by the fit model of Eq. (12), since the residual is dominated by white phase noise.   6. The modified Allan deviation and its error bar of the accumulated comagnetometer phase and the residual phases for the fit with a 7th order polynomial. To fulfill the MAD statistics criteria [31], only data are shown for integration time τ < 4000 s.

The theoretical statistical uncertainty bound
The theoretical limit of the 129 Xe EDM uncertainty can be derived as the CRLB, which also provides insights into optimizing experimental parameters. For the sake of simplicity, only a single species spin-precession signal is considered and its amplitude is assumed to be a constant over the whole sub-run. For the GPF method, d A ( 129 Xe) is estimated with two steps: The VP fitting to obtain the phase of each block and global phase fitting of the sub-runs. Therefore, the overall CRLB is the combination of the results of these two fits.
For the phase φ of a sinusoid embedded in white Gaussian noise (WGN) observed over one block with time being symmetrically around 0 s, the CRLB is where σ 2 w is the variance of the WGN, A the amplitude and N the number of data points in one block [32]. The CRLB for the parameters in the fit model Eq. (12) is the reciprocal of the Fisher information matrix where M is the number of segments in one sub-run, and J is the number of blocks in one segment. For the sake of simplicity, the standard polynomial is used in the fit model Eq. (12). Assuming JM k=1 Φ k EDM t i k = 0 for i going from 0 to g and the phase uncertainty δφ is a constant, the considered CRLB can be simplified to the so called ideal or uncorrelated CRLB as By substituting Eqs. (11), (17) and (18) into Eq. (19), and exploiting the periodic property of the constructed EDM function (see Fig. 5), for our case, where T = M JN ∆t is the total measurement time and ∆t = 1/f s is the sampling interval. Note that the number of segments M should be large enough to ensure the orthogonality between Φ k EDM and the polynomial functions. In case of an exponentially decaying amplitude A of the precession signal, the CRLB has to be calculated with Eq. (18). For the PC method, the CRLB on the 129 Xe EDM for M segments is derived in Ref. [26] as The PC method applies linear fits to the comagnetometer phases within one segment to derive the comagnetometer frequency of each segment, which requires the addition of an interception term as a starting phase, increasing the variance by a factor of four compared to a linear fit without interception term. In the GPF method the accumulated comagnetometer phases within one sub-run are analyzed in a single fit, therefore the uncertainty does not increase as the interception term is orthogonal to the EDM function (see Eq. (18)). Furthermore, the PC method requires the unweighted average of at least four segment frequencies, which increases its statistical uncertainty even further.

Results
A Monte-Carlo study was conducted to confirm that the GPF method can reach the higher sensitivity as shown by the CRLB compared to the PC method. Later, the GPF method was used to obtain the 129 Xe EDM from the data set as taken in Ref. [11] using the same channel and block length for analysis. As there were data sets in the 2017 and 2018 campaigns which were not useable with the PC method but could be analyzed with the GPF method, we gathered all data and optimized the analysis parameters to obtain the minimum uncertainty from the data. Ultimately, an improved upper limit of the 129 Xe EDM was derived using the unblinded data.

Monte-Carlo tests
The accumulated phase of each spin species for the sampling point j was generated as Φ j He,syn = tj 0 γ He B(t) + 2π(f He lin + u He e −t/T He 1 )dt, (24) where the drift of the background field B(t) was parametrized with a 4th order polynomial. f Xe/He lin represent the frequency shifts caused by the chemical shift and Earth's rotation. u Xe/He are the drift amplitudes of the respective precession frequencies. The frequency drift was modeled as exponentially decaying functions with the characteristic time of T 1 [21,22,25]. Thereby it was assumed that T 1 is larger than T 2 and its range is listed in Table 3. f EDM is the frequency shift due to the coupling of a synthetic EDM d syn with the electric field according to Eq. (3). Substituting Eqs. (23) and (24) into Eq. (10) results in the synthetic comagnetometer phase, whose time dependence is designed to mimic the measured data (for details see Appendix A.2). The exponentially decaying spin precession signals of 129 Xe and 3 He atoms can be described by with t j = j∆t , which is the time for the sampling point j.   Table 3. Three different kinds of noise were separately added into the synthetic data, including two WGN with σ = 154 fT, the standard deviation of the white noise in the real data, and σ = 154/5 = 30.8 fT, as well as real SQUID gradiometer noise. The overall EDM values obtained with the GPF method from the 18 synthetic sub-runs for four synthetic values d syn = (1, 2, 5, 10) × 10 −28 e cm are plotted in Fig. 7. The averaged overall EDM uncertainty for WGN data with σ = 154 fT is 1.74 × 10 −28 e cm, which is roughly a factor of 5 larger than that obtained from the data with σ = 30.8 fT and a factor of 1.1 higher than the calculated CRLB for these 18 sub-runs, which is 1.59 × 10 −28 e cm. This mainly results from the correlation between the EDM and the parameters of the polynomials in the phase fit. The uncertainty for the real noise is 1.85 × 10 −28 e cm, being similar to that for the white noise with σ = 154 fT. Most of the 1σ confidence intervals of the derived EDM cover the added EDM values d syn , showing that the GPF method is capable of accurately obtaining d syn ≥ 1 × 10 −28 e cm independent of the realistic noise level.

Statistical uncertainty
Applying the GPF method to the same data set of 41 runs (80 sub-runs) as analyzed by the PC method [11] and using the same channel and analysis parameters, the statistical uncertainty is decreased by a factor of 2.1 from 6.6 ×10 −28 e cm to 3.1 ×10 −28 e cm. Due to fewer constraints in the GPF method, runs with the number of segments M = 4n with n ∈ N or having SQUID jumps could be included in the data analysis, leading to a total of 45 runs (87 sub-runs). Furthermore, the segments with zero high voltage are included into the analysis. For the analysis, the block length is t b = 5 s, the threshold of the F -test is set to P min = 0.6 (see Appendix B) and the minimum order of the polynomial used in the fit is set to 4 in order to adequately describe the comagnetometer phase drift. The average polynomial order used for all sub-runs is 6.4 and the maximum order 13.
The overall result using the full data set is d A ( 129 Xe) = 1.1 ± 3.1 × 10 −28 e cm with χ 2 /dof = 115.5/86. As all sub-run measurements were taken with considerable different background noise a χ 2 /dof ≥ 1 can be expected. According to the PDG guidelines [33] we accounted for these random variations by scaling the statistical uncertainty with the factor χ 2 /dof = 1.16 leading to 3.6 × 10 −28 e cm. Bootstrapping [34] the 87 EDM measurements resulted in an estimate of the statistical uncertainty of 3.14 × 10 −28 e cm. Fig. 8 shows the derived EDM results per sub-run. Sorting all EDM measurements into groups based on the experimental parameters, such as the cell geometry, B 0 field direction, number and duration of segments and the gas pressure, shows no correlation between the deduced EDM value and these parameters, as can be seen in Fig. 9. Furthermore, no correlation between the chosen polynomial order and the derived sub-run EDM values was seen.

Systematic uncertainty
The systematic uncertainties of the two experiment campaigns were extensively studied in Ref. [11]. We applied the same analysis to the full data set used here and the derived systematic uncertainties are summarized in Table 4. The correction to the comagnetometer frequency drift of order higher than 1, as it has been done in Ref. [11], becomes obsolete for the GPF method since the model of Eq. (12) considers the higher order drifts implicitly.
As mentioned above, the GPF uses the full data set, including data during the high voltage rampings. Therefore the charging current can have an impact on the result in two ways. First the charging current could magnetize parts of the experimental equipment, and change the magnetic field seen by the spins. By this mechanism a false EDM may be generated. This effect has been carefully analyzed in Ref. [11] and has been adapted for the data set used for GPF (see Charging current in Table 4). Secondly, the charging currents just as the leakage currents will generate magnetic  fields, which are correlated with the electric field direction. This effect is only present during the ramping lasting for a few blocks per segment, ranging from 20 to 160 blocks. The impact of the charging current acting as a leakage current is calculated and turned out to be negligible, relative to the effect of leakage currents as given in Table 4.
We further looked for the potential effect of the comagnetometer drift and the vibrational noise with Monte-Carlo simulations and did not find observable systematic error, see Appendix A. The overall systematic uncertainty is the weighted average of the systematic uncertainties of the two measurement campaigns 2017 and 2018 using the reciprocal of its statistical variance as weights, yielding 2.0 × 10 −28 e cm. The final result, separating the statistical and systematic uncertainties, is d A ( 129 Xe) = 1.1 ± 3.6 (stat) ± 2.0 (syst) × 10 −28 ecm, from which we set an upper limit |d A ( 129 Xe)| < 8.3 × 10 −28 e cm at the 95% C.L. This reanalysis leads to a limit that is a factor of 1.7 smaller compared to the previous result [11] and a factor of 8.0 compared to the result in 2001 [16].

SUMMARY AND OUTLOOK
We proposed a global phase fitting method to analyze spin precession data. Applying the GPF method to the data set used in Ref. [11] yields a consistent result for d A ( 129 Xe) but a two times smaller statistical uncertainty compared to the PC method, as predicted by the theoretical CRLB analysis. Using additional data which had to be discarded for the PC method due to incomplete electric field patterns and optimizing the analysis parameters, the upper limit of the 129 Xe EDM improves by a factor of 1.7 to |d A ( 129 Xe)| < 8.3 × 10 −28 e cm at the 95% C.L. This enables 129 Xe to be used as a comagnetometer in future neutron EDM experiments [35] with a systematic error contribution down to |d A ( 129 Xe)| × γ n /γ Xe = 2.1 × 10 −27 e cm. Our GPF method relieves the demands on the physical model describing the comagnetometer frequency drift and could be generally used in similar spin precession experiments, such as the Lorentz-invariance test. By optimizing the experimental parameters for the GPF method (see Appendix C), the upper limit for d A ( 129 Xe) could be reduced even further, as planned for an upcoming EDM campaign with optimized high voltage pattern.  To investigate the potential systematic effect caused by the comagnetometer phase drift, we altered the drift amplitude u Xe and u He in Eqs. (23) and (24) in the synthetic phase data. Fig. 12 shows the derived EDM values as a function of the scale ratio of the drift amplitude. No distinct correlation between the obtained EDM value and the drift amplitude could be observed. Therefore, we did not assign a model dependent uncertainty for the comagnetometer drift when applying the GPF method.

B THE F -TEST THRESHOLD
The F -test threshold P min affects the polynomial order used in the GPF method, as listed in Table 5. The EDM values for various P min are overlapped within the 1σ statistical uncertainty and are all consistent with zero. Additionally, the upper limit of the 129 Xe EDM is almost insensitive to the threshold. We have chosen 0.6 as F -test threshold yielding the highest upper bound.

C DESIGN OF EXPERIMENTAL PARAMETERS
The number of segments M in one sub-run has a significant impact on the estimation uncertainty derived by the GPF method. According to the ideal CRLB, a smaller number of segments results in a lower uncertainty, shown as the red line in Fig. 13. To search for the optimum segment number, we used the synthetic comagnetometer phase data with added white Gaussian noise. The phase uncertainty increases with time and starts with 0.1 mrad. The time constants T 2 for 129 Xe atoms and 3 He atoms are a random number ranged from 8000 s to 9000 s. The total measurement time length is fixed to 38400 s, while M is varied from 2 to 64. The averaged EDM value over 100 runs for each M are plotted as the blue crosses. The fit uncertainty is larger than the ideal CRLB due to the correlation between the EDM function and the phase drift. The gap is reduced with the increase of M , since the orthogonality condition is satisfied better. A relatively flat optimum is found around M = 16. Note that this optimum value also depends on the total measurement time. A sub-run with longer measurement time calls for a higher number of segments, hence the optimum number for T = 6400 s and T = 64000 s is 8 and 64, respectively. The improved understanding of the comagnetometer frequency drift behavior may reduce the requirement on the segment number, thus significantly increasing the measurement sensitivity.