Revisiting ¹²⁹Xe electric dipole moment measurements applying a new global phase fitting approach

For ¹²⁹Xe atoms stored in a cell permeated by a uniform magnetic field $\overrightarrow {B}$ and an electric field $\overrightarrow {E}$, that is parallel to $\overrightarrow {B}$, their spin precesses at an angular frequency

$$\begin{equation}{\omega }_{\mathrm{X}\mathrm{e}}=\left\vert {\gamma }_{\mathrm{X}\mathrm{e}}B+\frac{{d}_{\mathrm{A}}\left({}^{129}\mathrm{X}\mathrm{e}\right)E}{\hslash {F}_{\mathrm{X}\mathrm{e}}}\right\vert ,\end{equation} \tag{ 1 }$$

where F_Xe = 1/2 is the total angular momentum number and γ_Xe is the gyromagnetic ratio of ¹²⁹Xe. The magnetic field B in equation (1) becomes an interference term for directly calculating ${d}_{\text{A}}\left({}^{129}\mathrm{X}\mathrm{e}\right)$ from ω_Xe. To overcome the experimental difficulties on controlling and measuring $\overrightarrow {B}$, comagnetometry was introduced with two collated species measured at the same time [15, 18, 24]. ³He is an ideal candidate for comagnetometry due to its potentially high SNR and its negligible EDM compared to ${d}_{\text{A}}\left({}^{129}\mathrm{X}\mathrm{e}\right)$ [25]. The weighted frequency difference between ¹²⁹Xe atoms and ³He atoms is defined as

$$\begin{equation}{\omega }_{\mathrm{c}\mathrm{o}}={\omega }_{\mathrm{X}\mathrm{e}}-\frac{{\gamma }_{\mathrm{X}\mathrm{e}}}{{\gamma }_{\mathrm{H}\mathrm{e}}}{\omega }_{\mathrm{H}\mathrm{e}},\end{equation} \tag{ 2 }$$

and commonly named the comagnetometer frequency. Here ω_He = |γ_He B| is the spin precession frequency of ³He atoms with γ_He being its gyromagnetic ratio. Therefore, ω_co can be written as

$$\begin{equation}{\omega }_{\mathrm{c}\mathrm{o}}=-\frac{2{d}_{\text{A}}\left({}^{129}\mathrm{X}\mathrm{e}\right)}{\hslash } \overrightarrow {E}\cdot \hat{B},\end{equation} \tag{ 3 }$$

showing that ω_co is independent of the magnitude of the background magnetic field but depends on its orientation $\hat{B}$ 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 [22, 23, 26]. 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₁ of the 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 reference [17], these difficulties can be omitted ⁵ .

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 typical values of 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 reference [27]. The spin precession signal of the transverse magnetization of ³He and ¹²⁹Xe was recorded by a dc-SQUID system comprising two channels (Z1, Z2) with Z2 sensitive to the same B-field direction as Z1 but 12 cm further away from the ¹²⁹Xe and ³He sample. The high voltage and leakage current between the two electrodes of the cell were monitored. In the 2017(2018) campaign, a background magnetic field B₀ 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 figure 2). In order to decrease potential crosstalk (bleed-through) from vibrational noise onto the ¹²⁹Xe and ³He frequencies, a software SQUID axial gradiometer (Z1−Z2) was used. The power of vibrational noise was decreased by a factor of about 70, at a cost of 2% lower signal amplitude as compared to the Z1 magnetometer signal.

Table 1. Starting amplitude A₀, transverse relaxation time T₂, background magnetic field B₀, electric field E₀ and segment length t_s of both measurement campaigns.

	2017	2018
A0,Xe/pT	30–40	70–80
A0,He/pT	4–5	20–25
${T}_{2}^{\mathrm{X}\mathrm{e}}$/s	6000–7000	8000–9000
${T}_{2}^{\mathrm{H}\mathrm{e}}$/s	7000–8000	8000–9000
B0/μT	2.6	3
E0/(kV cm−1)	2.75	4.13
ts/s	400 or 800	100–800

Zoom In Zoom Out Reset image size

The left panel of figure 2 shows the raw SQUID gradiometer signal in pT (grey) of one run from the 2018 campaign lasting 35 000 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 figure 2, the segments last 300 s and 600 s, respectively. The first sub-run ranging from 50 s to 12 400 s is used as an example in the data analysis section.

As mentioned above, one approach to mitigate the effect of the comagnetometer frequency drift is repetitively reversing the direction of the electric field $\overrightarrow {E}$. This allows to separate the impact of ${d}_{\text{A}}\left({}^{129}\mathrm{X}\mathrm{e}\right)$ on ω_co from other interference terms. The E modulation method has been applied in diverse EDM experiments with varied modulation patterns [12, 28]. For the PC method, the common E pattern for one sub-run consisted of 36 segments with an equal time interval t_s, and the sign of $\overrightarrow {E}$ changed according to the following sequence ±[0 + - - + - + + - - + + - + - - + 0, 0 - + + - + - - + + - - + - + + - 0].

The PC method determines the EDM value from averaging the comagnetometer frequencies ω_co from 2ⁿ ($n\in \mathbb{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 reference [27]).

The block length t_b is a free parameter with a suitable range from 1 s to 20 s, to exclude the amplitude decay and frequency drift and to ensure a good performance of the fit for our data [17]. The SQUID data in each block are fitted to the function

$$\begin{equation}y={a}_{\mathrm{X}\mathrm{e}}\enspace \mathrm{sin}\left({\omega }_{\mathrm{X}\mathrm{e}}t\right)+{b}_{\mathrm{X}\mathrm{e}}\enspace \mathrm{cos}\left({\omega }_{\mathrm{X}\mathrm{e}}t\right)+{a}_{\mathrm{H}\mathrm{e}}\enspace \mathrm{sin}\left({\omega }_{\mathrm{H}\mathrm{e}}t\right)+{b}_{\mathrm{H}\mathrm{e}}\enspace \mathrm{cos}\left({\omega }_{\mathrm{H}\mathrm{e}}t\right)+{a}_{i}\enspace \mathrm{sin}\left({\omega }_{i}t\right)+{b}_{i}\enspace \mathrm{cos}\left({\omega }_{i}t\right)+c+d\cdot t,\end{equation} \tag{ 4 }$$

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⁻¹ 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 [29], 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 equation (4), the time of each block is assigned to be symmetrical around zero from −t_b/2 to t_b/2.

Figure 4 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–25 Hz as shown in the right plot of figure 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) as well as a synthetic noise with a sinusoidal function (see reference [30]). The phase of each species for the block k in the range of [−π, π) can be obtained by

$$\begin{equation}{\phi }_{\text{m}}^{k}=\mathrm{A}\mathrm{r}\mathrm{g}\left({a}_{\text{m}}^{k}+\text{i}\cdot {b}_{\text{m}}^{k}\right),\end{equation} \tag{ 5 }$$

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 t_k = (k − 1/2)t_b. The time interval of the block k is defined as (t_k−1, t_k ). The parameter uncertainties $\delta {a}_{\text{m}}^{k}$ and $\delta {b}_{\text{m}}^{k}$ are estimated from the covariance matrix of the fit

$$\begin{equation}\mathrm{C}\mathrm{o}\mathrm{v}=\frac{{ \overrightarrow {r}}^{\prime }\cdot \overrightarrow {r}}{\nu {\left({J}^{\prime }\cdot J\right)}^{-1}},\end{equation} \tag{ 6 }$$

where $\overrightarrow {r}$ is the residual, ν is the number of degrees of freedom and J is the Jacobian matrix. The standard deviation of the derived phase $\delta {\phi }_{\mathrm{X}\mathrm{e}/\mathrm{H}\mathrm{e}}^{k}$ is

$$\begin{equation}\delta {\phi }_{\text{m}}^{k}=\frac{\sqrt{{\left({a}_{\text{m}}^{k}\cdot \delta {b}_{\text{m}}^{k}\right)}^{2}+{\left({b}_{\text{m}}^{k}\cdot \delta {a}_{\text{m}}^{k}\right)}^{2}-2{a}_{\text{m}}^{k}{b}_{\text{m}}^{k}\mathrm{C}\mathrm{o}\mathrm{v}\left({a}_{\text{m}}^{k},{b}_{\text{m}}^{k}\right)}}{{\left({a}_{\text{m}}^{k}\right)}^{2}+{\left({b}_{\text{m}}^{k}\right)}^{2}}.\end{equation} \tag{ 7 }$$

Equation (6) assumes that the residual $\overrightarrow {r}$ stems from the wideband white noise, which is a conservative approach for our case since the main signal in the residuals was the narrowband vibrational noise, leading to an overestimation of the uncertainty $\delta {\phi }_{\text{m}}^{k}$. However, the ratio between $\delta {\phi }_{\text{m}}^{k}$ for different blocks reflects the decaying SNR. Therefore, these estimated uncertainties are used as weights in the subsequent GPF routine.

Zoom In Zoom Out Reset image size

The accumulated phase ${{\Phi}}_{\text{m}}^{k}$ in a block k of the continuously precessing spins is the sum of the wrapped phase ${\phi }_{\text{m}}^{k}$ and a multiple of 2π

$$\begin{equation}{{\Phi}}_{\text{m}}^{k}={\phi }_{\text{m}}^{k}+2\pi {n}_{\text{m}}^{k},\end{equation} \tag{ 8 }$$

where the integer ${n}_{\text{m}}^{k}$ is determined as

$$\begin{equation}{n}_{\text{m}}^{k}=\frac{{{\Phi}}_{\text{m}}^{k-1}+{\omega }_{\text{m}}^{k-1}{t}_{\mathrm{b}}}{2\pi },\end{equation} \tag{ 9 }$$

rounded to the lower integer and ${n}_{\text{m}}^{1}=0$. Here, the frequencies ${\omega }_{\text{m}}^{k-1}$ are obtained by the fit of the block k using equation (4). If ${{\Phi}}_{\text{m}}^{k}-\left({{\Phi}}_{\text{m}}^{k-1}+{\omega }_{\text{m}}^{k-1}{t}_{\mathrm{b}}\right)$ is either > π or < −π, ${n}_{\text{m}}^{k}$ is incremented or decremented by one, respectively, to ensure a continuous phase evaluation. The standard deviation of the accumulated phase $\delta {{\Phi}}_{\text{m}}^{k}$ is equal to $\delta {\phi }_{\text{m}}^{k}$ as equation (8) does not introduce any additional uncertainty. According to equation (2), the evolved comagnetometer phase ${{\Phi}}_{\mathrm{c}\mathrm{o}}^{k}$ for each block k is determined by

$$\begin{equation}{{\Phi}}_{\mathrm{c}\mathrm{o}}^{k}={{\Phi}}_{\mathrm{X}\mathrm{e}}^{k}-\frac{{\gamma }_{\mathrm{X}\mathrm{e}}}{{\gamma }_{\mathrm{H}\mathrm{e}}}{{\Phi}}_{\mathrm{H}\mathrm{e}}^{k}.\end{equation} \tag{ 10 }$$

By integrating equation (3), the accumulated phase due to a hypothetical ¹²⁹Xe EDM d_h at the block k is

$$\begin{equation}{{\Phi}}_{\text{EDM}}^{k}=-{t}_{\text{b}}\frac{2{d}_{\mathrm{h}}}{\hslash }\sum\limits _{i=1}^{k}\left({ \overrightarrow {E}}_{i}\cdot \hat{B}\right),\end{equation} \tag{ 11 }$$

where ${ \overrightarrow {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 ${{\Phi}}_{\text{bias}}^{k}$ is calculated and then used to blind the comagnetometer phase ${{\Phi}}_{\mathrm{c}\mathrm{o},\mathrm{b}}^{k}={{\Phi}}_{\mathrm{c}\mathrm{o}}^{k}+{{\Phi}}_{\text{bias}}^{k}$ in order to avoid operator induced bias during process optimization.

The measured phase ${{\Phi}}_{\mathrm{c}\mathrm{o}}^{k}$ originates not only from the potential ¹²⁹Xe EDM, but also from other sources such as chemical shift [22, 27]. These contributions are phenomenologically parametrized by a polynomial of gth order. Hence, the comagnetometer phase is fitted with the function

$$\begin{equation}{{\Phi}}_{\text{fit}}^{k}=a{{\Phi}}_{\text{EDM}}^{k}+{p}_{0}+{p}_{1}{\tilde {P}}_{1}\left({t}_{k}\right)+{p}_{2}{\tilde {P}}_{2}\left({t}_{k}\right)+\cdots +{p}_{g}{\tilde {P}}_{g}\left({t}_{k}\right),\end{equation} \tag{ 12 }$$

where a, p₀, p₁, p₂, ..., p_g are the global fit parameters, and the atomic EDM of ¹²⁹Xe is ${d}_{\text{A}}\left({}^{129}\mathrm{X}\mathrm{e}\right)=a\cdot {d}_{\mathrm{h}}$. Here the time series t_k are normalized to the interval [0, 1] and shifted Legendre polynomials ${\tilde {P}}_{n}\left({t}_{k}\right)$ are applied to decrease the correlation between polynomial coefficients [31]. The fit was performed by using the iterative weighted least squares estimation method ⁶ with the built-in function nlinfit in MATLAB. Thereby the inverse values of the phase variances ${\left(\delta {{\Phi}}_{\mathrm{c}\mathrm{o}}^{k}\right)}^{2}$ are used as weights. Figure 5 shows the comagnetometer phase ${{\Phi}}_{\mathrm{c}\mathrm{o}}^{k}$, the fit phase ${{\Phi}}_{\text{fit}}^{k}$, and the EDM function ${{\Phi}}_{\text{EDM}}^{k}$ constructed from the measured E-field pattern of the exemplary sub-run.

Zoom In Zoom Out Reset image size

To determine the order needed for the polynomial function in equation (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

$$\begin{equation}{P}_{g,g+q}={\int }_{0}^{{F}_{g,g+q}}{P}_{\mathrm{F}}\left(F;q,N-g-q\right)\mathrm{d}F,\end{equation} \tag{ 13 }$$

where P_F is the probability density function of the F-distribution and N is the number of data points [33]. The upper bound of the integral is

$$\begin{equation}{F}_{g,g+q}=\frac{\left(N-g-q\right)\left({\chi }_{g}^{2}-{\chi }_{g+q}^{2}\right)}{q\cdot {\chi }_{g+q}^{2}}.\end{equation} \tag{ 14 }$$

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 correlation matrix for the exemplary sub-run (see figure 2) using 7th order (fulfilling the criteria for P_min = 0.6) is given in table 2. The correlated uncertainties of the parameters are determined as the square root of the diagonal of the covariance matrix, which inherently includes the uncertainty of the correlations between a and all polynomial parameters. The influences of these correlations to the estimation of a are small due to the orthogonality between the constructed function ${{\Phi}}_{\text{EDM}}^{k}$ and the polynomial function of the order up to n − 2 where 2ⁿ is the number of nonzero high-voltage segments. For the exemplary run, 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 ${{\Phi}}_{\text{EDM}}^{k}$ (with 2⁵ = 32 high-voltage segments). The derived uncertainty is in good agreement with the result using the log profile likelihood method (see reference [30]). We also applied the linear regression method with the model in equation (12) and obtained consistent results.

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

$$\begin{equation}{\sigma }_{f}\left(\tau \right)=\frac{1}{2\pi }\sqrt{\frac{\sum\limits _{j=1}^{P-3n-1}{\left(\sum\limits _{k=j}^{j+n-1}{{\Phi}}^{k+2n}-2{{\Phi}}^{k+n}+{{\Phi}}^{n}\right)}^{2}}{2{n}^{2}{\tau }^{2}\left(P-3n+1\right)}},\end{equation} \tag{ 15 }$$

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 [34]. The residual phase ${{\Phi}}_{\mathrm{c}\mathrm{o}}^{k}-{{\Phi}}_{\text{fit}}^{k}$ to the sub-run B881 is shown in figure 6 (left) with the adjacent histogram showing a Gaussian-like distribution. The residual generally increases with time due to the signal decay. The MAD of this sub-run is plotted on the right in figure 6. σ_f of ${{\Phi}}_{\mathrm{c}\mathrm{o}}^{k}$ reaches the minimum at the integration time τ of 550 s and then increases due to the comagnetometer frequency drift. For the residual phase of this sub-run, the MAD decreases with increasing integration time according to σ_f ∝ τ^−3/2 (dashed line in figure 6) over the full range, down to 0.4 nHz. This behaviour that is also observed in all other sub-runs is an indicator that the comagnetometer phase ${{\Phi}}_{\mathrm{c}\mathrm{o}}^{k}$ is adequately described by the fit model of equation (12), since the residual is white phase noise dominated.

Zoom In Zoom Out Reset image size

The theoretical limit of the ¹²⁹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}_{\text{A}}\left({}^{129}\mathrm{X}\mathrm{e}\right)$ is estimated with two steps: the VP fitting to obtain the phase of each block and GPF 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 with white Gaussian noise (WGN) observed over one block with time being symmetrically around 0 s, the CRLB is

$$\begin{equation}\mathrm{v}\mathrm{a}\mathrm{r}\left(\hat{\phi }\right){\geqslant}\frac{2{\sigma }_{\mathrm{w}}^{2}}{{A}^{2}N},\end{equation} \tag{ 16 }$$

where ${\sigma }_{\mathrm{w}}^{2}$ is the variance of the WGN, A the amplitude and N the number of data points in one block [35]. The CRLB for the parameters in the GPF model equation (12) is the reciprocal of the Fisher information matrix

$$\begin{equation}I=\left(\begin{matrix}\hfill \sum\limits _{k=1}^{JM}\frac{{\left({{\Phi}}_{\text{EDM}}^{k}\right)}^{2}}{\delta {\phi }_{k}^{2}}\hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \sum\limits _{k=1}^{JM}\frac{{{\Phi}}_{\text{EDM}}^{k}{t}_{k}^{g}}{\delta {\phi }_{k}^{2}}\hfill \\ \hfill \sum\limits _{k=1}^{JM}\frac{{{\Phi}}_{\text{EDM}}^{k}}{\delta {\phi }_{k}^{2}}\hfill & \hfill \sum\limits _{k=1}^{JM}\frac{1}{\delta {\phi }_{k}^{2}}\hfill & \hfill \dots \hfill & \hfill {\vdots}\hfill \\ \hfill {\vdots}\hfill & \hfill {\vdots}\hfill & \hfill \ddots \hfill & \hfill {\vdots}\hfill \\ \hfill \sum\limits _{k=1}^{JM}\frac{{{\Phi}}_{\text{EDM}}^{k}{t}_{k}^{g}}{\delta {\phi }_{k}^{2}}\hfill & \hfill \sum\limits _{k=1}^{JM}\frac{{t}_{k}^{g}}{\delta {\phi }_{k}^{2}}\hfill & \hfill \dots \hfill & \hfill \sum\limits _{k=1}^{JM}\frac{{t}_{k}^{g}{t}_{k}^{g}}{\delta {\phi }_{k}^{2}}\hfill \end{matrix}\right),\end{equation} \tag{ 17 }$$

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 here, the standard polynomial is used in equation (12). Assuming $\sum _{k=1}^{JM}{{\Phi}}_{\text{EDM}}^{k}{t}_{k}^{i}=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

$$\begin{equation}\mathrm{v}\mathrm{a}\mathrm{r}{\left({\hat{d}}_{\text{A}}\left({}^{129}\mathrm{X}\mathrm{e}\right)\right)}_{\text{GPF}}{\geqslant}\frac{\delta {\phi }^{2}}{\sum\limits _{k=1}^{JM}{\left({{\Phi}}_{\text{EDM}}^{k}\right)}^{2}}.\end{equation} \tag{ 18 }$$

By substituting equations (11), (16) and (17) into equation (18), and exploiting the periodic property of the constructed EDM function (see figure 5), for our case, $\sum _{k=1}^{JM}{\left({{\Phi}}_{\text{EDM}}^{k}\right)}^{2}=M\sum _{k=1}^{J}{\left({{\Phi}}_{\text{EDM}}^{k}\right)}^{2}$, the overall CRLB for ${d}_{\text{A}}\left({}^{129}\mathrm{X}\mathrm{e}\right)$ becomes

$$\begin{equation}\mathrm{v}\mathrm{a}\mathrm{r}{({\hat{d}}_{\text{A}}({}^{129}\mathrm{X}\mathrm{e}))}_{\text{GPF}}{\geqslant}\frac{2{\sigma }_{\mathrm{w}}^{2}}{{A}^{2}N}/({(\frac{2\vert E\vert {t}_{\text{b}}})}^{2}M\sum\limits _{k=1}^{J}{k}^{2})\end{equation} \tag{ 19 }$$

$$\begin{equation}{\geqslant}\frac{{\sigma }_{\mathrm{w}}^{2}}{{A}^{2}}{\left(\frac{\hslash }{2\vert E\vert }\right)}^{2}\frac{6{M}^{2}\delta t}{{T}^{3}},\end{equation} \tag{ 20 }$$

where T = MJNΔ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 ${{\Phi}}_{\text{EDM}}^{k}$ and the polynomial functions. In case of an exponentially decaying amplitude A of the precession signal, the CRLB has to be calculated with equation (17). For the PC method, the CRLB on the ¹²⁹Xe EDM for M segments is derived in reference [27] as

$$\begin{equation}\mathrm{v}\mathrm{a}\mathrm{r}{\left({\hat{d}}_{\text{A}}\left({}^{129}\mathrm{X}\mathrm{e}\right)\right)}_{\mathrm{P}\mathrm{C}}{\geqslant}\frac{{\sigma }_{\mathrm{w}}^{2}}{{A}^{2}}{\left(\frac{\hslash }{2\vert E\vert }\right)}^{2}\frac{24{M}^{2}\delta t}{{T}^{3}}.\end{equation} \tag{ 21 }$$

The PC method applies linear fits to the comagnetometer phases within each segment to derive the comagnetometer frequency, 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 single interception term is orthogonal to the EDM function (see equation (17)), showing up as the difference of equation (20) to equation (21). Furthermore, the PC method requires the unweighted average of at least four segment frequencies, which increases its statistical uncertainty even further.

The accumulated phase of each spin species for the sampling point j was generated as

$$\begin{equation}{{\Phi}}_{\mathrm{X}\mathrm{e},\mathrm{s}\mathrm{y}\mathrm{n}}^{j}={\int }_{0}^{{t}_{j}}{\gamma }_{\mathrm{X}\mathrm{e}}B\left(t\right)+2\pi \left({f}_{\text{lin}}^{\mathrm{X}\mathrm{e}}+{u}^{\mathrm{X}\mathrm{e}}\enspace {\text{e}}^{-t/{T}_{1}^{\mathrm{X}\mathrm{e}}}+{f}_{\text{EDM}}\right)\mathrm{d}t,\end{equation} \tag{ 22 }$$

$$\begin{equation}{{\Phi}}_{\mathrm{H}\mathrm{e},\mathrm{s}\mathrm{y}\mathrm{n}}^{j}={\int }_{0}^{{t}_{j}}{\gamma }_{\mathrm{H}\mathrm{e}}B\left(t\right)+2\pi \left({f}_{\text{lin}}^{\mathrm{H}\mathrm{e}}+{u}^{\mathrm{H}\mathrm{e}}\enspace {\text{e}}^{-t/{T}_{1}^{\mathrm{H}\mathrm{e}}}\right)\mathrm{d}t,\end{equation} \tag{ 23 }$$

where the drift of the background field B(t) was parametrized with a 4th order polynomial. ${f}_{\text{lin}}^{\mathrm{X}\mathrm{e}/\mathrm{H}\mathrm{e}}$ 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₁ [22, 23, 26]. Thereby it was assumed that T₁ is larger than T₂ with range as 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 equation (3). Substituting equations (22) and (23) into equation (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 ¹²⁹Xe and ³He atoms can be described by

$$\begin{equation}{V}_{\mathrm{X}\mathrm{e}/\mathrm{H}\mathrm{e}}^{j}={A}_{0}^{\mathrm{X}\mathrm{e}/\mathrm{H}\mathrm{e}}\enspace {\text{e}}^{-{t}_{j}/{T}_{2}^{\mathrm{X}\mathrm{e}/\mathrm{H}\mathrm{e}}}\enspace \mathrm{sin}\enspace {{\Phi}}_{\mathrm{X}\mathrm{e}/\mathrm{H}\mathrm{e},\mathrm{s}\mathrm{y}\mathrm{n}}^{j}\end{equation} \tag{ 24 }$$

with t_j = jΔt, which is the time for the sampling point j.

Table 3. The range of the parameter values used for generating synthetic spin precession data for Monte-Carlo simulations.

Para.	Range
uHe	3.5–4.5 μHz
uXe	9–11 μHz
${f}_{\text{lin}}^{\mathrm{H}\mathrm{e}}$	4–10 μHz
${T}_{1}^{\mathrm{H}\mathrm{e}}$	9000–14 000 s
${T}_{1}^{\mathrm{X}\mathrm{e}}$	9000–14 000 s
${f}_{\text{lin}}^{\mathrm{X}\mathrm{e}}$	4–10 μHz

The parameters used to generate synthetic data were taken from 18 sub-runs of high sensitivity from the 2018 campaign. The starting amplitude of ¹²⁹Xe and ³He were set to ${A}_{0}^{\mathrm{X}\mathrm{e}}=70$ pT, ${A}_{0}^{\mathrm{H}\mathrm{e}}=25$ pT and ${T}_{2}^{\mathrm{X}\mathrm{e}/\mathrm{H}\mathrm{e}}=8000$ s. The electric field contained 36 segments of 200 s up to 800 s length, as used in the measurement campaign. The values of other parameters in equations (22) and (23) were random and uniformly distributed in the ranges listed in table 3. Three different kinds of noise were separately added into the synthetic data of equation (24), including two WGN with σ = 154 fT, the standard deviation of the white noise in the real data, and a 5 times lower value of σ = 30.8 fT as well as real SQUID gradiometer noise derived as shown in appendix A.1. The overall EDM values obtained with the GPF method from the 18 synthetic sub-runs for four synthetic values d_syn = (1, 2, 5 and 10) × 10⁻²⁸ e cm are plotted in figure 7. The averaged overall EDM uncertainty for WGN data with σ = 154 fT is 1.74 × 10⁻²⁸ 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⁻²⁸ 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⁻²⁸ e cm, being similar to that for the white noise with σ = 154 fT, confirming the negligible influence of vibrational noise crosstalk. 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 obtaining d_syn ⩾ 1 × 10⁻²⁸ e cm in the presence of the chosen noise.

Zoom In Zoom Out Reset image size

4.2.1. Statistical uncertainty

Applying the GPF method to the same data set of 41 runs (80 sub-runs) as analyzed by the PC method [17] and using the same channel and analysis parameters, the statistical uncertainty is decreased by a factor of 2.1 from 6.6 × 10⁻²⁸ e cm to 3.1 × 10⁻²⁸ e cm.

Due to fewer constraints in the GPF method, runs with the number of segments M ≠ 4n with $n\in \mathbb{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 a zero voltage were included in the analysis. For the analysis, the block length t_b = 5 s was chosen and the threshold for the F-test set to P_min = 0.6 (see appendix B) and the minimum order of the polynomial was set to 4 in order to adequately describe the comagnetometer phase drift. The average polynomial order used for all sub-runs was 6.4 and the maximum order 13.

The overall averaged result using the full data set is ${d}_{\text{A}}\left({}^{129}\mathrm{X}\mathrm{e}\right)=1.1{\pm}3.1{\times}1{0}^{-28}\enspace \text{e}\enspace \mathrm{c}\mathrm{m}$ with χ²/dof = 115.5/86. The detailed results and parameters for each sub-run are given in reference [30]. According to the PDG guidelines [36] we accounted for these random variations by scaling the statistical uncertainty with the factor $\sqrt{{\chi }^{2}/\mathrm{dof}}=1.16$ leading to 3.6 × 10⁻²⁸ e cm. To cross-check the derived confidence interval of the weighted mean on ${d}_{\text{A}}\left({}^{129}\mathrm{X}\mathrm{e}\right)$, we applied the bootstrap method with 5000 data samples [37] to the 87 EDM values and obtained a similar estimate of the statistical uncertainty of 3.14 × 10⁻²⁸ e cm. Figure 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₀ 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 figure 9. Furthermore, no correlation between the chosen polynomial order and the derived sub-run EDM values was seen.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2.2. Systematic uncertainty

The systematic uncertainties of the two experiment campaigns were extensively studied in reference [17]. 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 one (due to the unweighted average of four segments in the PC method), as it has been done in reference [17], becomes obsolete for the GPF method since the model of equation (12) considers the higher order drifts implicitly. Furthermore, when applying the GPF method one has to check the additional impact of the charging current acting as leakage currents since it uses data during the high-voltage ramps. This effect is calculated in appendix C and turned out to be negligible.

Table 4. The systematic uncertainties determined as done in reference [17] based on the data set used for the GPF method.

	2017 (e cm)	2018(e cm)
Leakage current (incl. impact of ICharging during ramping)	1.20 × 10−28	4.40 × 10−31
Magnetization due to charging current	1.70 × 10−29	1.20 × 10−29
Cell motion (rotation)	4.20 × 10−29	4.0 × 10−29
Cell motion (translation)	2.60 × 10−28	1.90 × 10−28
|E|2 effect	1.20 × 10−29	2.20 × 10−30
|E| uncertainty	9.90 × 10−29	5.70 × 10−30
Geometric phase	⩽2 × 10−31	⩽1 × 10−29
Total systematic uncertainty	3.07 × 10−28	1.95 × 10−28
Scaled statistical uncertainty	15.57 × 10−28	3.67 × 10−28

We further looked for the potential systematic effect caused by the comagnetometer drift and the vibrational noise with Monte-Carlo simulations and did not find observable systematic errors, 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⁻²⁸ e cm. The final result, separating the statistical and systematic uncertainties, is

$$\begin{equation}{d}_{\text{A}}\left({}^{129}\mathrm{X}\mathrm{e}\right)=\left[1.1{\pm}3.{6}_{\left(\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}\right)}{\pm}2.{0}_{\left(\mathrm{s}\mathrm{y}\mathrm{s}\mathrm{t}\right)}\right]{\times}1{0}^{-28}\enspace \text{e}\enspace \mathrm{c}\mathrm{m},\end{equation} \tag{ 25 }$$

from which we can set an upper limit $\vert {d}_{\text{A}}\left({}^{129}\mathrm{X}\mathrm{e}\right)\vert {< }8.3{\times}1{0}^{-28}\enspace \text{e}\enspace \mathrm{c}\mathrm{m}$ at the 95% C.L. This reanalysis leads to a limit that is a factor of 1.7 smaller compared to the previous result [17] and a factor of 8.0 compared to the result in 2001 [15].

The effect of the real measurement noise (e.g. vibrational noise) on the estimated phase is quantitatively analyzed. Here, a synthetic precession signal was generated, using a single sinusoidal function with a constant amplitude A = 30 pT and a length of 10 000 s yielding 2000 blocks. Furthermore, white noise with σ = 154 fT (the standard deviation of the white noise in real gradiometer data) generated with MATLAB, or real noise (from the exemplary sub-run with a total noise power of (0.4 pT)² and the precession signals filtered out) were added separately to the simulated precession signal. The error of the fitted phase for block i is defined as _i = ϕ_fit,i − ϕ_real,i . Here ϕ_real,i was known and ϕ_fit,i was obtained from the fit to block i.

The histograms of _i for these two synthetic data sets are plotted in figure A1. The phase error for the white noise data is in good agreement with the normal distribution with σ = 1.11 × 10⁻⁴ rad, which is close to 1.08 × 10⁻⁴ rad, the CRLB on the phase estimator in equation (16). The phase error for the real noise data also satisfies the Gaussian distribution with σ = 1.24 × 10⁻⁴ rad, which is similar to the result of the white noise data. This implies that the vibrational noise did not cause evident additional phase deviation, although the standard deviation of the vibrational noise is around 7 times bigger than the white noise. As evident in figure A1, the vibrational noise does not cause an observable systematic error on the derived phase. Modelling the vibrational noise with a constant sinusoidal signal provides another way to study the phase error caused by the crosstalk effect. By doing so, we found that the error due to the modelled vibrational signal is more than ten-fold smaller than the uncertainty caused by the white noise. Detailed analysis results can be found in reference [30].

Zoom In Zoom Out Reset image size

The analyzed comagnetometer phase drift ${{\Phi}}_{\mathrm{c}\mathrm{o}}^{k}$ for 9 sub-runs with high sensitivities on ¹²⁹Xe EDM are plotted in the left panel of figure A2. Even though the subtraction of chemical shift and Earth's rotation is unnecessary for the EDM analysis using the GPF method, we subtracted both effects from the measured phase using the deterministic equations [18] to show only the comagnetometer phase drift. The synthetic comagnetometer phase drifts generated with two exponential functions for 18 random sub-runs are plotted in the right panel of figure A2, showing a similar behaviour as the experimentally obtained ones. The parameter ranges are listed in table 3.

Zoom In Zoom Out Reset image size

To investigate the potential systematic effect caused by the comagnetometer phase drift, we changed the ranges of the drift amplitudes u^Xe and u^He as given in table 3 by factors of 0.05 to 20 in the synthetic phase data. White phase noise was added into the phase data and the noise deviation σ increased with time in an exponential way with T₂ = 8000 s. Figure A3 shows the derived EDM values as a function of the scale ratio of the drift amplitudes for two starting noise uncertainties σ₀ = 100 μrad, similar to the real data, and σ₀ = 1 μrad. No significant dependence between the obtained EDM value and the drift amplitudes could be observed. For the low noise case, the maximum deviation of the central EDM value among 7 results is less than 1 × 10⁻²⁹ e cm. So, we did not assign a model dependent uncertainty for the comagnetometer drift when applying the GPF method, also supported by the analysis shown in figure 7 and the F-test results as shown in appendix B.

Zoom In Zoom Out Reset image size