# Correlative methods for dual-species quantum tests of the weak equivalence principle

B Barrett, L Antoni-Micollier, L Chichet, B Battelier, P-A Gominet, A

Bertoldi, Philippe Bouyer, A Landragin

► To cite this version:

B Barrett, L Antoni-Micollier, L Chichet, B Battelier, P-A Gominet, et al.. Correlative methods for dual-species quantum tests of the weak equivalence principle. New Journal of Physics, 2015, 17 (8), pp.085010. 10.1088/1367-2630/17/8/085010. [hal-01183822v2](https://hal.sorbonne-universite.fr/hal-01183822v2)

HAL Id: hal-01183822

https://hal.sorbonne-universite.fr/hal-01183822v2

Submitted on 10 Mar 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Correlative methods for dual-species quantum tests of the weak equivalence principle

This content has been downloaded from IOPscience. Please scroll down to see the full text.

2015 New J. Phys. 17 085010

(http://iopscience.iop.org/1367-2630/17/8/085010)

View the table of contents for this issue, or go to the journal homepage for more

Download details:

IP Address: 134.157.80.136

This content was downloaded on 10/03/2016 at 14:48

Please note that terms and conditions apply.
<figure><img src="image_2.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

29 March 2015

REVISED

23 June 2015

ACCEPTED FOR PUBLICATION

17 July 2015

PUBLISHED

26 August 2015

Content from this work
may be used under the
terms of the [Creative
Commons Attribution 3.0
licence](http://creativecommons.org/licenses/by/3.0/).

Any further distribution of
this work must maintain
attribution to the
author(s) and the title of
the work, journal citation
and DOI.

<figure><img src="image_3.png" /></figure>

PAPER

# Correlative methods for dual-species quantum tests of the weak equivalence principle

B Barrett¹, L Antoni-Micollier¹, L Chichet¹, B Battelier¹, P-A Gominet¹, A Bertoldi¹, P Bouyer¹ and A Landragin²

¹ LP2N, IOGS, CNRS and Université de Bordeaux, rue François Mitterrand, F-33400 Talence, France

² LNE-SYRTE, Observatoire de Paris, CNRS and UPMC, 61 avenue de l'Observatoire, F-75014 Paris, France

E-mail: brynle.barrett@institutoptique.fr

**Keywords:** atom interferometry, equivalence principle, data analysis

## Abstract

Matter-wave interferometers utilizing different isotopes or chemical elements intrinsically have different sensitivities, and the analysis tools available until now are insufficient for accurately estimating the atomic phase difference under many experimental conditions. In this work, we describe and demonstrate two new methods for extracting the differential phase between dual-species atom interferometers for precise tests of the weak equivalence principle (WEP). The first method is a generalized Bayesian analysis, which uses knowledge of the system noise to estimate the differential phase based on a statistical model. The second method utilizes a mechanical accelerometer to reconstruct single-sensor interference fringes based on measurements of the vibration-induced phase. An improved ellipse-fitting algorithm is also implemented as a third method for comparison. These analysis tools are investigated using both numerical simulations and experimental data from simultaneous ⁸⁷Rb and ³⁹K interferometers, and both new techniques are shown to produce bias-free estimates of the differential phase. We also report observations of phase correlations between atom interferometers composed of different chemical species. This correlation enables us to reject common-mode vibration noise by a factor of 730, and to make preliminary tests of the WEP with a sensitivity of 1.6 × 10⁻⁶ per measurement with an interrogation time of T = 10 ms. We study the level of vibration rejection by varying the temporal overlap between interferometers in a symmetric timing sequence. Finally, we discuss the limitations of the new analysis methods for future applications of differential atom interferometry.

## 1. Introduction

Einstein's equivalence principle (EEP) is a fundamental concept in physics that describes the exact correspondence between the gravitational and inertial mass of any object. It is a central assumption of the theory of general relativity—which interprets gravity as a geometrical feature of space–time, and predicts identical accelerations for different objects in the same gravitational field. Precise tests of the EEP are of great interest in various fields of physics. For instance, some theories that attempt to unify gravity with the other fundamental forces predict a violation of this principle [1, 2]. The detection of such a violation could aid our understanding of dark energy in cosmology, and advance the search for physics beyond the standard model. In contrast, null results are also pivotal for putting bounds on model parameters contained in various extensions to general relativity [3–5]. The equivalence principle is generally divided into three sub-principles that each must be satisfied for the EEP to hold [6, 7]: the local Lorentz invariance, the local position invariance and the weak equivalence principle (WEP). In this article, we will focus on the latter.

The WEP—otherwise known as the universality of free fall—states that a charge-free body will undergo an acceleration in a gravitational field that is independent of its internal structure or composition. Tests of the WEP generally involve measuring the relative acceleration between two different test bodies that are in free fall with the same gravitational field. The WEP is characterized by the Eötvös parameter, $\eta$, given by
$$ \eta = 2 \frac{a_1 - a_2}{a_1 + a_2} = \frac{\Delta a}{a}, \qquad (1) $$

where $a_1$ and $a_2$ are the accelerations of the two bodies, $\Delta a = a_1 - a_2$ is the relative acceleration, and $a = (a_1 + a_2)/2$ is the average acceleration. The WEP is satisfied if and only if $\Delta a = 0$—implying that $\eta = 0$.

The most precise tests of the WEP have been carried out with lunar laser ranging techniques [8], or using a rotating torsion balance [9, 10], which have both measured $\eta$ at the level of a few parts in $10^{13}$. Various Space missions to test the WEP at improved levels ($10^{-15}$ or better) using other classical devices are presently in progress [11–13]. On a separate frontier, a number of groups have carried out tests between cold atoms [14–18] in an effort to probe the WEP at the quantum level. The majority of these tests have been conducted using matter-wave interferometers which, over the past few decades, have been extensively studied both theoretically and experimentally [19–22]. Atom interferometers have been utilized as ultra-precise inertial sensors to measure, for example, the gravitational acceleration $g$ [23–26], the gravitational constant $G$ [27–29], gravity gradients [25, 30–33], gravitational field curvature [34], and rotations [35–38]. A WEP test based on atom interferometry involves measuring the differential phase shift resulting from a relative acceleration between two species with different masses that are in free fall within the same gravitational field. This measurement is based on the same principle as gravity gradiometry, where the quantity of interest is the differential phase between test atoms of the same type but in different spatial locations. The gradient of the gravitational field can be extracted from the differential phase between two sources, while higher derivatives of the field can be accessed if more than two sources are used. This technique was recently demonstrated to measure the curvature of the gravitational field, and has been proposed to detect gravitational waves and to study geophysical effects [39–41]. Presently, the state-of-the-art for WEP tests using matter-wave interferometry corresponds to an uncertainty of $3 \times 10^{-8}$ [18]. A comparison between the gravitational acceleration measured by atoms and a macroscopic object (i.e. a falling corner-cube) have also been carried out, and yield agreement at the level of $\delta\eta \simeq 6.5 \times 10^{-9}$ [42]. A handful of ground-based [43–46] and micro-gravity-based [41, 47–50] cold-atom experiments are currently underway that aim to greatly improve this precision. In addition, there have been a number of proposals for Space-based quantum tests of the WEP [7, 51–54] that target accuracies at the level of $10^{-15}$.

So far, most tests with cold atoms have used two isotopes of the same atomic element, e.g. $^{85}$Rb and $^{87}$Rb [14, 15, 18, 55], or $^{87}$Sr and $^{88}$Sr [16]. Although this class of test bodies has demonstrated a good level of common-mode noise rejection when performing differential phase measurements [15], it is intrinsically less sensitive to possible violations of the equivalence principle because the two atoms are relatively similar in mass and composition. Thus, it is interesting to perform these tests with two entirely different atomic elements. In this article, we will focus on the case of $^{87}$Rb and $^{39}$K. These atoms exhibit a large difference in their number of nuclei —facilitating a mass ratio of $M_{Rb}/M_K \sim 2.2$. Additionally, they have identical hyperfine spin structure, and similar excitation wavelengths (around 780 and 767 nm, respectively), which enables the use of the same laser technology and optics for cooling and interferometry. Dual-species interferometers of this type have the added advantage of being highly independent—that is, atomic sample properties such as the size and temperature, or interferometer parameters such as the interrogation time, Raman phase, and detuning, can be controlled independently. In contrast to dual-isotope setups where many of these parameters are coupled, this feature is ideal for studying a variety of systematic effects that will be important for future precision measurements [55]. For a more complete comparison of alkali atoms as candidates for WEP tests, see for example [56].

One complication that arises with non-common elemental species is a difference in the scale factors, $S_j \simeq k_j^{\text{eff}} T_j^2$, between the interferometers. When the interrogation times $T_j$ are the same, this difference originates from the effective wave vectors $k_j^{\text{eff}}$ of the interferometer beams used for atoms $j = 1$ and 2. Assuming that the WEP is true, the phase shift of the two interferometers due to a common acceleration $a$ is $\Phi_j = S_j a$. Thus, a difference in the scale factors produces a relative phase shift between interferometers for the same acceleration: $\delta\phi_a^{\text{sys}} = (S_1 - S_2)a$. For the case of $^{85}$Rb and $^{87}$Rb, the scale factors can be made the same by a suitable choice of Raman laser detuning that guarantees $k_1^{\text{eff}} = k_2^{\text{eff}}$ [52]. However, this is not generally possible for different chemical elements, and the systematic phase shift must be addressed in other ways.

Another issue related to having different scale factors regards the rejection of common-mode vibration noise between interferometers. From an analysis of the interferometer transfer functions (see appendix C or [52, 57, 58], for instance), one can show that perfect common-mode rejection requires four conditions to be satisfied: (i) the interferometers occur simultaneously with $T_1 = T_2$, such that they experience the same vibration noise, (ii) they have identical wave vectors, $k_1^{\text{eff}} = k_2^{\text{eff}}$, they exhibit (iii) identical effective Rabi frequencies, $\Omega_1^{\text{eff}} = \Omega_2^{\text{eff}}$, and (iv) identical pulse durations, $\tau_1 = \tau_2$. These conditions imply that if $S_1 \neq S_2$, the interferometers do not respond to common-mode noise with the same phase shift.

The scale factors can be made the same by adjusting the interrogation times of the interferometers such that $T_1 = rT_2$, where $r = \sqrt{k_2^{\text{eff}}/k_1^{\text{eff}}}$ [48]. This technique eliminates the systematic phase shift $\delta\phi_a^{\text{sys}}$ resulting from a constant acceleration, and improves the rejection of common vibration noise at frequencies $\lesssim 1/T_1$, but it
degrades the rejection efficiency at frequencies above ~1/T₁ (see appendix C). However, if the ratio r is very close to unity, as it is for some choices of atoms (r ≈ 1.009 for ³⁹K and ⁸⁷Rb), this option represents a good compromise between efficient noise rejection and reducing systematic effects.

In this article, we describe and demonstrate three analysis methods for atom-interferometric WEP tests—including two new techniques that eliminate both aforementioned problems of systematic phase shifts and diminished common-mode rejection between coupled interferometers of different atomic species. The first of these two new methods is a generalized Bayesian analysis of the Lissajous curves formed by plotting the coupled sensor measurements parametrically. The second technique involves restoring the interferometer fringes by correlating with an auxiliary mechanical accelerometer. In this case, the phase shift for each species can be measured directly from the reconstructed fringes regardless of their scale factors or the degree of temporal overlap between the interferometers. Both of these new methods intrinsically account for different scale factors, and return unbiased estimates of the differential phase. Finally, to give a complete picture, we compare these techniques with an improved ellipse-fitting method recently developed by Szpak et al [59]. This numerical procedure yields an estimate of the differential phase shift with reduced bias compared to more commonly implemented algorithms in the presence of significant amounts of uncorrelated noise between sensors.

In this work, we also report correlated phase measurements between simultaneous interferometers of different chemical species (³⁹K and ⁸⁷Rb). When operated in an environment with significant levels of vibration noise, we demonstrate a common-mode vibration rejection factor of γ ≈ 730. These results represent a major step toward precise tests of the WEP with elements exhibiting vastly different masses. We also investigate the accuracy of the three aforementioned methods on experimental data obtained from the K-Rb interferometer.

The article is organized as follows. Section 2 reviews some theoretical background concerning a WEP test with a dual-species interferometer. In section 3, we briefly describe the three methods of extracting the differential phase. We give a brief description of the experimental setup for the K-Rb interferometer in section 4. We present our experimental results in section 5, and we give a discussion of the advantages and limitations of the new methods in section 6. Finally, we conclude in section 7. A detailed description of the three analysis methods, including extensive numerical tests of the generalized Bayesian estimator, can be found in the Appendices.

## 2. Testing the WEP with two atomic species

An atom-interferometric test of the WEP involves measuring the relative acceleration between two atoms of different mass. This can be done in one of two ways: (i) the absolute acceleration of each atom, $a_1$ and $a_2$, can be individually measured and subtracted, or (ii) $\Delta a$ can be measured directly from the differential phase, $\phi_d$. In the ideal case, acceleration measurements are performed simultaneously in order to take advantage of correlated noise between sensors—reducing the total uncertainty in $\Delta a$. Since method (ii) involves a direct measurement of $\phi_d$, it intrinsically requires both simultaneity and phase correlation between atomic sensors to reject common-mode noise. Henceforth, two or more atom interferometers that satisfy these conditions are referred to as ‘coupled sensors’. Method (i) can be carried out regardless of these two constraints. In this section, we outline some theoretical background related to a WEP test with method (ii).

Generally, the output from two coupled atomic sensors is described by the following sinusoids

$$y_1(a) = A_1 \cos(S_1 a + \phi_1) + B_1, \quad (2a)$$

$$y_2(a) = A_2 \cos(S_2 a + \phi_2) + B_2, \quad (2b)$$

where $A_j$ and $B_j$ are, respectively, the amplitude and offset of the interferometer fringes associated with sensor j ($j = 1, 2$). In principle, these two parameters can be measured and equations (2) can be recast in the normalized form $n_j = (y_j - B_j)/A_j$:

$$n_1(a) = \cos(S_1 a + \phi_1), \quad (3a)$$

$$n_2(a) = \cos(S_2 a + \phi_2). \quad (3b)$$

Here, a is an acceleration common to both atoms, $S_j$ is the scale factor for interferometer j, and $\phi_j$ is a phase shift. The scale factors can be computed exactly from the integral of the response function, $f_j(t)$, given by equation (C.5):

$$S_j = k_j^{\text{eff}} \int f_j(t) dt = k_j^{\text{eff}} (T_j + 2\tau_j) \left( T_j + \frac{4\tau_j}{\pi} \right), \quad (4)$$
where $k_j^{\text{eff}}$ is the effective wave-vector for the counter-propagating interferometer beams, $T_j$ is the interrogation time, and $\tau_j$ is the $\pi/2$ Raman pulse duration. A detailed explanation of the response function and its role in WEP tests is outlined in appendix C. For large interrogation times, $T_j \gg \tau_j$, the scale factors reduce to the well-known relation $S_j \approx k_j^{\text{eff}} T_j^2$.

The phases $\phi_j$ each have three contributions, one from the interferometer laser phase $\phi_j^{\text{laser}}$, one from parasitic systematic effects $\phi_j^{\text{sys}}$, and one from the WEP signal $\phi_j^{\text{WEP}}$

$$ \phi_j = \phi_j^{\text{laser}} + \phi_j^{\text{sys}} + \phi_j^{\text{WEP}}. \qquad (5) $$

In an experiment, the laser phase is a control parameter which can be set to zero, and systematic effects are independently nullified as much as possible. The shift due to a WEP violation can be defined as $\phi_j^{\text{WEP}} = S_j(a_j - a)$, which is expected to be very close to zero. In the ideal case, the total interferometer phase $\Phi_j$ contains only the shift due to the mean acceleration, $S_j a$, and a WEP violation. It then follows that

$$ \Phi_j = S_j a + \phi_j = S_j a_j. \qquad (6) $$

In this case, the total phase of each interferometer can be related to the Eötvös parameter in the following way

$$ \eta = \frac{\Phi_1/S_1 - \Phi_2/S_2}{a} = \frac{\Delta a}{a}. \qquad (7) $$

In principle, the sensitivity in this type of WEP test increases as the square of the interrogation time, $T \sim T_1 \sim T_2$, due to the scale factors, $S_j$, that appear inversely in equation (7).

The general form of equations (3) describes a Lissajous curve. For the purposes of this analysis, it is useful to redefine the phases in equations (3) to reduce the number of free parameters. Choosing sensor 2 as a reference to rescale the phase of sensor 1, we define a common phase $\phi_c$ that satisfies

$$ \phi_c \equiv S_2 a + \phi_2, \qquad (8a) $$

$$ \kappa \phi_c + \phi_d \equiv S_1 a + \phi_1. \qquad (8b) $$

Here, we introduce two new parameters—the scale factor ratio $\kappa$ and the differential phase $\phi_d$—which are constrained from equation (8) to be

$$ \phi_d = \phi_1 - \kappa \phi_2, \qquad \kappa = \frac{S_1}{S_2}. \qquad (9) $$

The sensor outputs are now recast with $\phi_c$ as the primarily parameter

$$ n_1(\phi_c) = \cos(\kappa \phi_c + \phi_d), \qquad (10a) $$

$$ n_2(\phi_c) = \cos(\phi_c). \qquad (10b) $$

Comparing equations (6), (7) and (9), it follows that the Eötvös parameter is directly proportional to the differential phase:

$$ \eta = \frac{\phi_d}{S_1 a}. \qquad (11) $$

### 3. Correlative methods of differential phase extraction

In this section, we review three different methods to measure the differential phase from experimental data: ellipse fitting, Bayesian analysis and fringe reconstruction from mirror acceleration measurements.

#### 3.1. Improved ellipse fitting

The ellipse fitting technique was first applied to atom interferometry in [60] for situations in which the phase common to two coupled atomic sensors is sufficiently scrambled to impede individual fringe observation. In this case, when the measurements from each sensor are plotted parametrically, one obtains an ellipse that is free from common phase noise. Using a least-squares ellipse fitting algorithm, the differential phase $\phi_d$ can be extracted. Multiple groups have demonstrated the utility of ellipse fitting for measurements of gravity gradients [32, 33] and the gravitational constant G [27, 28, 61]. However, this technique suffers from a number of drawbacks. First, it is valid only for coupled sensors with the same scale factor ($\kappa = 1$). Second, in the presence of moderate
<figure><img src="image_4.png" /><figcaption>**Figure 1.** The process of Bayesian estimation of the differential phase $\phi_d$ from synthetic data following a Lissajous curve. (a) Randomly chosen points (labelled 1–4) following equations (10) with Gaussian noise added to $n_1$ and $n_2$. For the actual Lissajous curve (shown as the black curve), we chose $\kappa = 0.8$ and $\phi_d^{\text{act}} = 1$ rad for illustrative purposes. (b) The prior probability distribution computed from Bayes’ algorithm after each measurement. The vertical solid line indicates the differential phase used in the simulation. (c) Error in the Bayesian estimate after successive measurements. Points represent the difference between $\phi_d^{\text{est}}$ and $\phi_d^{\text{act}}$ (i.e. the systematic error), where $\phi_d^{\text{est}}$ is the Bayesian estimate based on the maximum likelihood value from the corresponding prior distribution. The error bars indicate the statistical uncertainty, which are computed from the standard deviation of the prior distributions shown in (b).</figcaption></figure>

amounts of noise in the fringe offsets or amplitudes (the parameters $A_j$ and $B_j$ in equations (2)), or in the differential phase, the ellipse fit returns a biased estimate of $\phi_d^{\text{est}}$.

Recently, Szpak *et al* [59] developed an algorithm based on the optimization of the approximate maximum likelihood distance which seeks a balance between costly geometric methods and stable algebraic techniques. This algorithm—termed the FGEF method—exhibits a smaller bias in the differential phase estimate over a relatively large phase range (centered on $\pi/2$) compared to the more commonly used ‘direct ellipse fit’ (DEF) technique [62]. Additionally, [63] includes error estimations for the geometrically meaningful ellipse parameters (center coordinates, axes and orientation). We have extended their work to include an estimate of the statistical uncertainty in the differential phase, $\delta\phi_d$. We provide a more detailed comparison between DEF and FGEF methods of ellipse fitting in appendix A.

## 3.2. Generalized Bayesian analysis

Heuristic approaches to estimating the differential phase, such as ellipse-fitting methods, do not have knowledge of the noise present in experimental data, nor of how various types of noise can affect the outcome of measurements. Bayesian analysis offers an efficient alternative to the problem by constraining the estimate based on a statistical model that describes the distribution of data that results from different noise sources [64]. Bayesian phase estimation was studied in the context of atom interferometry in [65] for two sensors containing the same scale factor ($\kappa = 1$). In that work, a detailed study of each possible noise source (amplitude, offset and differential phase) is presented. Reference [56] also used Bayesian analysis to estimate the differential phase from a hypothetical system with $\kappa < 1$. There, however, only noise in the differential phase is considered, and the range of common phase was constrained to $\phi_c \in [0, \pi]$. To the best of our knowledge, no complete Bayesian estimator exists that (i) is valid for any scale factor ratio, (ii) accounts for noise in all relevant system parameters, and (iii) allows $\phi_c$ to vary over a broad range. Furthermore, this type of analysis has not yet been demonstrated on experimental data from dual-species interferometers.

In this work, we have developed a generalized Bayesian estimator for $\phi_d$—based on the approach of [65]— that satisfies all three of the requirements mentioned above. We demonstrate this technique by measuring $\phi_d$ from both simulated data (see appendix B) and experimental data from our K–Rb interferometer (see section 5). The advantage of using this estimation technique is that the uncertainty in $\phi_d$ converges much faster than other methods (i.e. it scales as $\sim 1/\sqrt{N}$, where N is the number of measurements), so fewer data are required to reach a given level of sensitivity. Furthermore, since $\kappa$ is built directly into the Bayesian estimate of $\phi_d$, it is free from the aforementioned systematic phase shift $\delta\phi_d^{\text{sys}}$ arising between interferometers with different scale factors. However, some of the drawbacks of the Bayesian analysis are that it requires a priori knowledge of the noise in the system, and it is computationally costly due to the large number of integrals that must be evaluated.

Figure 1 illustrates the basic Bayesian estimation procedure. Here, we simulate data that follow the Lissajous equations (10) with added Gaussian noise in the sensor offsets. After each successive measurement from the system, the width of the new ‘prior’ probability distribution decreases and additional peaks are suppressed—

<footnote n="3">³ Rosi *et al* [34] demonstrated that the bias in the estimate of $\phi_d$ can be eliminated under certain conditions when fitting an ellipse in three dimensions from the output of three simultaneous interferometers.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Schematic of the FRAC method for a single interferometer. A mechanical accelerometer is mounted to the back of the retro-reflection mirror for the Raman beam (with wave vectors $k_1$ and $k_2$). Acceleration measurements during the interferometer sequence are weighted by the response function $f(t)$ (see equation (C.5)) and integrated to obtain the phase estimate $\phi_d^{\text{vib}}$. Correlating this phase with the interferometer signal during the same time interval reproduces the interference fringe.</figcaption></figure>

facilitating an improvement in the estimate of $\phi_d$. This is how the Bayesian method builds in information from previous measurements. It is clear from figure 1(c) that after only a few iterations, both the statistical and systematic error in $\phi_d$ have decreased dramatically. A detailed description of the generalized Bayesian analysis can be found in appendix B.

### 3.3. Fringe reconstruction by accelerometer correlation—the differential FRAC method

Differential atom interferometry is often utilized under conditions where each sensor is overwhelmed by external phase noise that is common to both sensors. Typically, one is concerned with only the differential phase and not the common phase $\phi_c$, which is treated as an arbitrary parameter. Both the ellipse-fitting and Bayesian estimation methods for extracting $\phi_d$ take this approach. An alternative technique involves measuring the common phase and correcting for it. For the case of parasitic mirror vibrations, single-sensor interference fringes that are otherwise smeared by phase noise can be restored based on measurements from seismometers [66–68] or mechanical accelerometers [41, 69, 70]. Henceforth, we refer to this as the fringe reconstruction by accelerometer correlation (FRAC) method. In this work, we demonstrate how the FRAC method can be applied to two quasi-simultaneous interferometers of different atomic species to measure the relative phase shift between them. This technique to extract $\phi_d$ is referred to as the *differential FRAC method* throughout the article to differentiate between the (standard) FRAC method, which is generally employed to measure the absolute phase shift of a single atom interferometer.

Figure 2 illustrates the basic schematic of the FRAC method for a single interferometer. A mechanical accelerometer is secured to the back of the reference mirror used to retro-reflect interferometry light, and the time-dependent mirror acceleration, $a^{\text{vib}}(t)$, is recorded during the interferometer sequence. These acceleration measurements are first weighted by the response function of the $j$th interferometer, $f_j(t)$, and are then integrated to find the vibration-induced phase given by<fnref n="4" />

<footnote n="4">⁴ The underlying assumption of this technique is that motion of the reference mirror at frequencies within the response bandwidth of the mechanical accelerometer are responsible for phase changes of the atom interferometer. Although the corresponding acceleration signal is indistinguishable from fluctuations in local gravity (as a consequence of the equivalence principle), we can be confident in our assumption since typical variations in gravity occur on timescales much larger than $2T$.</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 3. (a) Schematic illustrating the source of correlation in the differential FRAC method. The signal from each atom interferometer (AI), exhibiting a known scale factor S<sub>j</sub> and unknown phase shift φ<sub>j</sub>, is directly correlated (indicated by the solid lines) with a common vibration phase measured by a mechanical accelerometer (MA). Since the resulting fringes are derived from a common source, measurements of each φ<sub>j</sub> are highly correlated, as indicated by the dotted lines and characterized by the positive coefficient ϱ<sub>φ<sub>1</sub>,φ<sub>2</sub></sub>. This technique can be extended for multiple coupled interferometers (shown in gray), although we focus on the case of only two. (b) Correlation coefficient ϱ<sub>φ<sub>1</sub>,φ<sub>2</sub></sub> as a function of φ<sub>d</sub> estimated from a large sample of simulated fringes. Here, we assumed κ = 1 and we added non-common Gaussian noise to the phase of the reconstructed fringes with a standard deviation of σ<sub>φ<sub>d</sub></sub> = 0.1 rad.</figcaption></figure>

$$ \phi_j^{\text{vib}} = k_j^{\text{eff}} \int f_j(t) a^{\text{vib}}(t) dt. \qquad (12) $$

For each repetition of the experiment, this random phase is computed and correlated with the interferometer signal. This process allows one to reconstruct the interference fringes point-by-point. Depending on the level of vibrations and the interferometer sensitivity, the range of vibration-induced phases can span multiple fringes—enabling the single-sensor phase shift φ<sub>j</sub> to be measured using, for instance, a sinusoidal least-squares fit to the data. It is straightforward to extend this algorithm for two or more interferometers, which do not need to be overlapped in time. In this case, the only additional requirement is that the time-series of mirror acceleration measurements span the interrogation times for all interferometers. For two coupled sensors, the differential phase is easily computed from the individual sensor phase shifts via φ<sub>d</sub> = φ<sub>1</sub> − κφ<sub>2</sub>. The statistical error in this quantity is governed by

$$ (\delta\phi_d)^2 = (\delta\phi_1)^2 + (\kappa\delta\phi_2)^2 - 2\kappa\rho_{\phi_1,\phi_2}(\delta\phi_1)(\delta\phi_2), \qquad (13) $$

where the δφ<sub>j</sub> represent the statistical uncertainties in the φ<sub>j</sub> obtained from fits to the two fringes, and ϱ<sub>φ<sub>1</sub>,φ<sub>2</sub></sub> is the correlation coefficient for the measurements of φ<sub>1</sub> and φ<sub>2</sub>. In the limit of perfect correlation (ρ<sub>φ<sub>1</sub>,φ<sub>2</sub></sub> = 1), the uncertainty in the differential phase reduces to δφ<sub>d</sub> = |δφ<sub>1</sub> − κδφ<sub>2</sub>|. Figure 3(a) illustrates how the coupled-interferometer correlation is utilized by the differential FRAC method. Since the fringes for each interferometer are recovered using measurements from the same classical device, the phase noise present on each fringe is highly correlated. This induces a correlation between the measurements of φ<sub>1</sub> and φ<sub>2</sub> extracted from the fits, as characterized by ϱ<sub>φ<sub>1</sub>,φ<sub>2</sub></sub>. The key to the differential FRAC method is maximizing this correlation to reduce the uncertainty in φ<sub>d</sub>.

The correlation coefficient for a given set of reconstructed fringes can be estimated numerically from a large sample of simulated data. We find that it is sensitive to experimental parameters such as the level of uncorrelated noise on each sensor, the scale factor ratio and the differential phase. For instance, figure 3(b) shows the dependence of ϱ<sub>φ<sub>1</sub>,φ<sub>2</sub></sub> on φ<sub>d</sub> for synthetic fringes that contain non-common phase noise with a standard deviation of 0.1 rad. The correlation coefficient yields a maximum when the interferometers are perfectly in-phase or π radians out-of-phase. This is an ideal feature for WEP tests, since the maximum sensitivity occurs exactly at the expected signal of φ<sub>d</sub> = 0. This implies that, unlike ellipse-fitting methods where the sensitivity is optimized at φ<sub>d</sub> = π/2, one does not need to engineer an additional phase shift between the atoms to optimize the sensitivity and reduce systematic bias. Furthermore, a recent study of a gradiometer configuration (i.e. κ = 1) has shown that the differential FRAC method can reach sensitivities close to the quantum-projection-noise limit when modest levels of uncorrelated phase noise are present [71].

A number of ideal features make this technique interesting for both absolute and differential atom interferometry experiments.

(1) The differential FRAC estimate of φ<sub>d</sub> is precise and unbiased over the full phase range φ<sub>d</sub> ∈ [0, π], since it relies on least-squares fits to individual fringes.

(2) It is simple, fast, and computationally low in cost—allowing the interferometer phase to be corrected in real-time [70], or by post-processing the data [66, 67, 69].
(3) Unlike the Bayesian analysis, the FRAC method does not require any *a priori* information about the interferometer offsets, contrasts, and noise parameters—which can be challenging to measure accurately without phase stability [72].

(4) Systematic phase shifts in $\phi_d$ due to non-identical pulse durations $\tau_j$ and Rabi frequencies $\Omega_j^{\text{eff}}$ [15] are accounted for in the estimates of $\phi_j^{\text{vib}}$ for each interferometer. Such systematics will be important to consider in future long-baseline differential interferometry experiments [7, 44–46, 51–54].

(5) The relative timing between coupled interferometers can be freely chosen—they need not be overlapped. This is a unique feature to dual-species interferometers that do not share the same Raman beams. Unlike the ellipse-fitting and Bayesian techniques, the FRAC method allows one to extract absolute phase information from each sensor. Varying the temporal overlap between interferometers can be useful for studying a variety of effects, such as the level of correlation between sensors, or systematics related to the interaction between atoms [55].

(6) Single-sensor fringes can be accurately measured in ‘noisy’ environments, which is ideal for mobile sensors such as atomic gravimeters [26, 41, 68, 69].

Although the standard FRAC method is conceptually simple to implement, the drawback is that it is sensitive to errors in the measurements of vibrations. Such errors include the quality of coupling between the mirror and the mechanical device, electronic noise in the signal acquisition, the level of self-noise of the device, drifts in the offset or sensitivity factor, and nonlinearities in both the amplitude and frequency response. The natural low-pass filtering feature of atom interferometers can alleviate some of these effects—particularly those that dominate at frequencies beyond the cut-off frequency $1/2T_j$. Alternatively, the mechanical accelerometer can be used as a course measurement of the phase to identify the correct interferometer fringe. Then one can use the atom interferometer output to refine the phase measurement [69, 70] by inverting the sinusoidal relation (2). On the other hand, measurements of $\phi_d$ using the differential FRAC method are much less sensitive to many of these noise sources since they are common to two simultaneous interferometers. We discuss the limitations of this method in more detail in section 6.

## 4. Description of the Interférométrie Cohérente pour l'Espace (ICE) experiment

ICE is an experiment that aims to measure $\eta$ using a dual-species interferometer of $^{87}$Rb and $^{39}$K. It is designed to be transportable and to operate in the micro-gravity environment provided by the Novespace Zero-g plane [41, 48, 69, 73]. In this section, we give a brief description of the experimental setup.

A detailed description of the telecom-frequency fiber-based laser system used on ICE can be found in [41, 74]. For each atomic species, we utilize a master-slave architecture, where the master laser diode is locked to either a saturated absorption peak (in the case of rubidium), or to a frequency comb (in the case of potassium). The slave lasers are frequency-locked to their corresponding master through an optical beat-note in the 1550 nm telecom band. After second harmonic generation to 780 nm for $^{87}$Rb and 767 nm for $^{39}$K, the frequency of each slave laser can be precisely adjusted over ~1.3 GHz within ~2 ms of settling time. Approximately 1.5 W of total light is available in each slave beam before entering a free-space optical bench. This module is composed of a series of shutters and acousto-optic modulators (AOMs) that are used to split, pulse and frequency shift the light appropriately for cooling, state preparation, interferometry and detection. Finally, the 780 and 767 nm light is coupled into a series of single-mode, polarization-maintaining fibers and sent to the vacuum chamber. The two frequencies required for cooling and repumping, as well as driving Raman transitions in $^{87}$Rb, are generated via a broadband fiber-based electro-optic modulator operating near 6.8 GHz. Similarly, an AOM operating in dual-pass configuration at ~230 MHz is used to generate these frequencies for $^{39}$K.

The sensor head is composed of a non-magnetic titanium vacuum chamber surrounded by a μ-metal shield. The chamber resides within three nested Helmholtz coils used to compensate residual magnetic fields and to generate a bias along the vertical axis. A custom 2-to-6 way fiber splitter is used to combine the 780 and 767 nm light intended for laser cooling without significant power loss via a polarizing cube and a dichroic wave plate. The splitter subsequently divides the light equally into six beams that are re-coupled into independent fibers used for the dual-species vapor-loaded magneto-optical trap (MOT). In a similar way, light for both detection and interferometry is overlapped in a free-space 2-to-1 way fiber combiner for 780 and 767 nm. The ~2 cm diameter beams output from the combiner have the same linear polarization, and are aligned along the vertical direction through the atoms. A quarter-wave plate (fabricated for the intermediate wavelength of 773 nm and mounted in front of the retro-reflection mirror) rotates the polarization of the Raman beams by 90° such that the counter-propagating fields have $\lin$ ⊥ $\lin$ polarization.
<figure><img src="image_7.png" /><figcaption>Figure 4. Symmetric timing sequence for the K–Rb interferometer. The first Raman pulse of the ³⁹K interferometer is delayed by $\alpha_{Rb} + \Delta T_{Rb, K}$ relative to that of ⁸⁷Rb such that the central π-pulse of both interferometers occurs at the same time. The preparation stages (shown in blue) include both the internal state preparation pulses (microwave + push beam for ⁸⁷Rb, depump for ³⁹K), and the external state preparation via time-of-flight. Raman pulses are shown in green, and detection pulses in orange.</figcaption></figure>

A typical experimental sequence for the K–Rb interferometer is shown in figure 4 and is carried out as follows. The MOT beams load approximately 2 × 10⁸ (7 × 10⁷) atoms in 0.5 s, which is followed by a 7 ms (5 ms) molasses cooling stage for the ⁸⁷Rb (³⁹K) sample. In addition to cooling, the rubidium molasses stage also pumps the atoms into the $|F=2\rangle$ ground state. This is followed by a microwave π-pulse that transfers atoms into $|F=1, m_F=0\rangle$, and the remaining atoms are removed with a push beam resonant with the $F=2$ to $F'=3$ transition. During the potassium molasses, the frequency and intensity of the cooling and repump beams are modified in a similar manner to [75, 76]. At the end of the molasses, the atoms are in a superposition of both hyperfine ground states, which is a critical part of the cooling mechanism for potassium [75]. We detune our 767 nm push beam to the red of the $F=2$ to $F'=3$ transition by ~17 MHz (2.9 Γ) to optically pump the atoms into the $F=1$ level with a 3 µs pulse. Following this depumping stage, the ³⁹K atoms are distributed roughly equally amongst the magnetic sub-levels of the lower hyperfine ground state. With this system, we achieve temperatures of ~3 µK for ⁸⁷Rb and ~20 µK for ³⁹K, as confirmed by both time-of-flight imaging and velocity-sensitive Raman spectroscopy. After preparing the internal atomic states, we typically wait ~12 ms for the atoms to fall such that the Doppler resonance of both sets of counter-propagating Raman beams becomes non-degenerate. Additionally, we apply an external magnetic bias field between 1 – 2 Gauss to shift the $|F=1, m_F = \pm 1\rangle$ states of potassium away from the central $m_F=0$ state on which we perform interferometry. The frequency of the Raman beams for both species is detuned by –1.2 GHz (–200 Γ) relative to the $F=2$ to $F'=3$ transition. We then apply the interferometry pulses in a symmetric fashion, such that the central π-pulse for both interferometers occurs at the same time, as shown in figure 4. The delay between the π/2 pulses for either atom, $\Delta T_{Rb, K}$, can be adjusted within the interrogation time of the rubidium interferometer, $T_{Rb}$, in order to study correlations and effects related to the scale factor ratio, $\kappa$. Finally, we measure the atomic state populations for each atom via fluorescence detection on an avalanche photodiode (50 MHz bandwidth) within 100 µs of one another.

## 5. Experimental results

We now describe some experimental results obtained from the K–Rb interferometer. All of the data presented in this work were recorded in a laboratory environment, with the interferometer beams aligned along the vertical direction, and with no anti-vibration platform. To compensate for the Doppler shift due to gravity, the frequency difference between Raman beams for interferometers $j=1 \equiv K$ and $j=2 \equiv Rb$ is chirped at a rate of $\alpha_j \approx k_j^{\text{eff}} g$ to account for the gravity-induced Doppler shift of the falling atoms. This modifies the total phase shift of the interferometers from equations (2), $\Phi_j = S_j a + \phi_j$, to the following

$$ \Phi_j = S_j (a - \alpha_j / k_j^{\text{eff}}) + \phi_j \approx (k_j^{\text{eff}} g - \alpha_j) T_j^2. \quad (14) $$

The last expression represents the case when both interferometers experience the same acceleration, $a_j = a = g$, and the scale factors can be approximated as $S_j \approx k_j^{\text{eff}} T_j^2$. Determining the location of the central fringe, for which $\alpha_j = k_j^{\text{eff}} g$ is fixed for all $T_j$, yields a measurement of $g$. Using the fact that the sensitivity of the interferometer scales as $T_j^2$, absolute measurements of $g$ have been demonstrated at accuracies of a few 10⁻⁹ [24, 68, 77, 78].
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Interferometer fringes from the $^{87}$Rb interferometer operating at $T = 25$ ms at a contrast of ~40% without vibration isolation. The interferometer phase is scanned by both uncontrolled mirror motion, and by varying the chirp rate $\alpha$ between Raman beams about a central value of $\alpha_0 = 25.1355$ MHz s$^{-1}$. The open circles indicate raw measurements of the normalized atomic population in the $|F = 2\rangle$ state, while the closed circles are the same measurements after applying the FRAC phase correction, $\phi^{\text{vib}}$, to each point. A few of these corrections are shown as blue arrows. The solid curve is a least-squares fit to the corrected data, resulting in a signal-to-noise ratio of ~30 and a relative statistical uncertainty of $10^{-7}$ in the determination of $g_{\text{Rb}}$—corresponding to almost an order of magnitude improvement compared to the raw data.</figcaption></figure>

As discussed in the introduction, we are interested in measuring the differential acceleration $\Delta a$ between $^{39}$K and $^{87}$Rb. One way of achieving this is to measure the gravitationally-induced accelerations $g_K$ and $g_{Rb}$ from each interferometer independently by scanning the chirp rates, $\alpha_j$, in a low-noise environment. This is the approach recently employed for WEP tests with $^{39}$K and $^{87}$Rb by Schlippert et al [17]. However, at high levels of sensitivity (i.e. large $T_j$), or in ‘noisy’ environments, mirror vibrations can corrupt the fringes—making individual phase measurements more challenging. We now demonstrate the utility of the FRAC technique for measuring $g$ from a single interferometer under these conditions.

There are typically two approaches in which the FRAC method can be applied to restore the interference fringes of a single interferometer. The first approach is to let the interferometer phase be ‘scanned’ randomly by vibrations while the laser-induced phase is held fixed. The reconstructed fringes in this case are purely a function of $\phi_j^{\text{vib}}$, as shown in figure 2. This mode of operation can be used to precisely calibrate the mechanical accelerometer by rescaling the voltage-to-acceleration sensitivity factor of the device such that the fringe period is$^5 2\pi$. The second approach is to scan the interferometer phase in a controlled manner, for example by varying the phase difference between Raman lasers, and to correct each phase using $\phi^{\text{vib}}$ obtained during the same measurement interval. This procedure is illustrated in figure 5, where the fringes of a $T = 25$ ms $^{87}$Rb interferometer are shown before and after applying the FRAC correction. Here, the interferometer is operated without any vibration isolation in the presence of a root mean squared (rms) dc vibration noise of $a_{\text{rms}}^{\text{vib}} \approx 6 \times 10^{-5}$ m s$^{-2}$ (integrated over the frequency response of the interferometer)—corresponding to an rms phase noise of $\phi_{\text{rms}}^{\text{vib}} = k_{\text{Rb}}^{\text{eff}} a_{\text{rms}}^{\text{vib}} T_{\text{Rb}}^2 \approx 0.6$ rad. Acceleration measurements were performed with a force-balance three-axis accelerometer (Nanometrics Titan, dc to 430 Hz bandwidth, 5 V/g sensitivity). By applying the FRAC correction to these data, we improve the signal-to-noise ratio (SNR) and hence the uncertainty in the central fringe measurement by almost an order of magnitude. We estimate an individual phase correction uncertainty of $\delta\phi^{\text{vib}} = 1/\text{SNR} \approx 33$ mrad based on the improved SNR of ~30. With this method, we emphasize that the interferometer sensitivity is directly linked to the intrinsic noise of the accelerometer + signal acquisition system, and the quality of the coupling between the device and the Raman mirror. Therefore, modest improvements to any of these system components can result in a dramatic increase in the fringe SNR.

## 5.1. K–Rb interferometer correlation

Typically, when mirror motion is the dominant source of phase noise it is advantageous to use differential atom interferometry techniques to measure $\Delta a$ through the differential phase $\phi_d$. This requires a high level of correlation between interferometers in order to reject the common-mode phase noise. We now compare three methods of extracting $\phi_d$ from experimental data recorded in an environment with high vibrational noise, as in the case of onboard applications [58, 69]. These studies are also applicable to future high-sensitivity differential interferometers operated in low-noise environments [44–46].

<footnote>$^5$One advantage of performing this procedure is that the device can be precisely calibrated for the vibration spectrum on site. Depending on the bandwidth and spectral response of the device, the sensitivity can vary significantly with the vibration spectrum.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Measurements of normalized $|F = 2\rangle$ state populations from simultaneous K–Rb interferometers operating at $T \approx 3$ ms. Fringes are scanned randomly by applying vibrational noise to the retro-reflection mirror. Graphs (a) and (b) show a time series of 500 measurements from the $^{87}$Rb and $^{39}$K interferometers, respectively, along with histograms of the populations. (c) Atomic populations from (a) and (b) plotted parametrically—indicating strong correlation between the two species. The solid green line is an ellipse fit to the data using the FGEF method, which yields a differential phase $\phi_d^{\text{ellipse}} = 1.13(2)$ rad. A separate estimate from a Bayesian analysis gives $\phi_d^{\text{Bayes}} = 1.18(2)$ rad. (d) Interferometer fringes reconstructed from measurements of mirror motion using the FRAC method. The red and blue curves correspond to least-squares fits to Rb and K data, respectively. The differential phase estimated from the fits is $\phi_{\text{rms}}^{\text{FRAC}} = 1.17(1)$ rad. Other interferometer parameters: pulse separations: $T_{\text{Rb}} = 3.018$ ms, $T_{\text{K}} = 3$ ms; $\pi/2$-pulse durations: $\tau_{\text{Rb}} = 4$ µs, $\tau_{\text{K}} = 6$ µs; delay between interferometers: $\Delta T_{\text{K, Rb}} = 10$ µs; one-photon Raman detunings: $\Delta_{\text{Rb}} = \Delta_{\text{K}} \approx -1.2$ GHz.</figcaption></figure>

Figure 6 shows data produced by quasi-simultaneous K–Rb interferometers at a total interrogation time of $2T = 6$ ms. Here, we held the chirp rate fixed at $\alpha_f \approx k_j^{\text{eff}} g$ for each species, and we applied strong vibrations to the system ($a_{\text{rms}}^{\text{vib}} \approx 0.05 \text{ m s}^{-2}$) such that the random vibration-induced phase $\phi_j^{\text{vib}}$ spanned multiple fringes ($\phi_{\text{rms}}^{\text{vib}} \approx 7.3$ rad). The vibrations were applied by mounting a heavy industrial fan on top of the support structure surrounding the vacuum system and running it during the experiment. Figure 6(a) shows a histogram of $^{87}$Rb $|F = 2\rangle$ population measurements, $y_{\text{Rb}}$, which clearly indicates the characteristic bimodal probability distribution of a sinusoid. These distributions can be used to estimate the contrast, offset and SNR of the interferometer fringes as described in [69]. We note that the bimodal distribution is less pronounced for $^{39}$K in figure 6(b) owing to a smaller fringe contrast, and thus a lower SNR, compared to $^{87}$Rb. Despite this fact, the two sensors exhibit strong correlations, as confirmed by the ellipse in figure 6(c).

For these experimental parameters the scale factor ratio is $\kappa = S_{\text{K}}/S_{\text{Rb}} = 1.008$, and the Lissajous curve formed by parametrically plotting the atomic state populations, $y_{\text{Rb}}$ and $y_{\text{K}}$, is indistinguishable from an ellipse at the present level of offset noise. We measure a differential phase of $\phi_d^{\text{ellipse}} = 1.13(2)$ rad from a least-squares fit to an ellipse using the FGEF method [59]. We also estimate $\phi_d^{\text{Bayes}} = 1.18(2)$ rad using the Bayesian analysis described in section 3.2 and appendix B. Here, it is worth mentioning that this non-zero differential phase does not originate from a WEP violation, but from systematic phase shifts in the experiment—primarily due to the quadratic Zeeman effect from an external magnetic bias field ($\sim 1$ G) that is used to sufficiently split the ground state magnetic sub-levels in $^{39}$K.

Figure 6(d) shows the output of each interferometer as a function of the vibration-induced phase, $\phi_j^{\text{vib}}$. Here, the single-sensor fringes were reconstructed using the FRAC method using mirror vibration measurements from a broadband micro-electro-mechanical accelerometer (Colibrys SF3600, dc to 1 kHz bandwidth, 1.2 V/g sensitivity). From these data the differential phase shift between interferometers is clearly visible. Sinusoidal least-squares fits to each fringe yield $\phi_d^{\text{FRAC}} = \phi_{\text{K}} - \kappa\phi_{\text{Rb}} = 1.17(1)$ rad. Here, the statistical uncertainty $\delta\phi_d^{\text{FRAC}}$ was computed from the quadrature sum of each interferometer phase error. The value of $\phi_d$ estimated from the Bayesian analysis and the FRAC method are in good agreement. On the other hand, the differential phase from the ellipse fit is underestimated by $\sim 40$ mrad, i.e. $2\sigma$ below $\phi_d^{\text{Bayes}}$ and $\phi_d^{\text{FRAC}}$. We attribute this discrepancy to the inherent bias of ellipse-fitting techniques (see appendix A), which increases with the level of offset noise or differential phase noise in either interferometer.

We emphasize that a crucial input parameter for the Bayesian analysis is the common phase range. We use the accelerometer data to estimate this range once the experiment is complete: $\phi_c \in [\min(\phi_{\text{Rb}}^{\text{vib}}), \max(\phi_{\text{Rb}}^{\text{vib}})]$. However, if an accelerometer is not available, it is also possible to estimate this range using the raw data from a single interferometer. For example, one can reduce the interrogation time until the sensitivity to vibrations reaches a point where interference fringes are clearly visible. By measuring the rms scatter of the phase about a reference sinusoid, one can estimate the level of vibration noise via the relation $a_{\text{rms}}^{\text{vib}} = \phi_{\text{rms}}^{\text{vib}}/S$. Once $a_{\text{rms}}^{\text{vib}}$ is known, this relation can be inverted to determine the range of phase scanned by the same level of vibrations at larger sensitivities/interrogation times.
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 7.** Correlation studies for different levels of temporal overlap. The symmetric, quasi-simultaneous K–Rb interferometer was operated with $T_{Rb} = 10$ ms and the interrogation time for potassium was varied between $T_K = 6 - 10$ ms. The interferometer phase was scanned by externally applied vibrations and individual fringes were restored using the FRAC method. Parametric plots of the atomic populations are shown to the right, along with the expected Lissajous curve (solid green line). These curves result from plotting the fit functions to each reconstructed fringe parametrically. There is a clear disagreement between the predicted Lissajous curves and the data for $T_{Rb} - T_K \gtrsim 2$ ms. Other interferometer parameters: $\tau_{Rb} = \tau_K = 3~\mu s$; $\Delta_{Rb} = \Delta_K \approx -1.2$ GHz.</figcaption></figure>

The data shown in figure 6(d) also indicate that the combined differential-atomic-sensor + mechanical-accelerometer system is capable of efficiently rejecting common vibrational noise. We estimate a rejection factor of $\gamma = k^{\text{eff}} a_{\text{rms}}^{\text{vib}} T^2 / \delta\phi_d^{\text{FRAC}} \approx 730$ for these data.

Figure 7 displays the results of a correlation study between rubidium and potassium interferometers operating at a total interrogation time of $2T = 20$ ms. Similar to figure 6, the interferometer phases are scanned by externally applied vibrations ($a_{\text{rms}}^{\text{vib}} \approx 1.7 \times 10^{-3} \text{ m s}^{-2}$, $\phi_{\text{rms}}^{\text{vib}} \approx 2.7 \text{ rad}$ at $T_{Rb} = T_K = 10$ ms). The vibrations were applied using the same method as previously mentioned, but with the fan set on a slower rotation setting. Here, we vary the interrogation time of potassium, $T_K$, in a symmetric way with respect to rubidium such that the centers of the $\pi$-pulses coincide. This optimizes the degree to which the vibration-induced phase noise remains common-mode, while modifying the degree of temporal overlap between interferometers. It also allows us to control the scale factor ratio since $\kappa$ scales at $(T_K/T_{Rb})^2$.

From figure 7, three features are clearly visible as $T_K$ is decreased. First, the potassium fringes undergo a phase shift that modifies the differential phase relative to the rubidium fringes. This feature, along with the fact that the scale factor ratio is varied, causes the shape of the Lissajous figures to change, as shown by the solid green curves. Second, the phase range scanned by the potassium interferometer reduces, since it scales as $T_K^2$. Finally, the level of correlation between the interferometers degrades as the temporal overlap decreases. This is evident from the lack of agreement between the data and the predicted Lissajous curves, particularly for $T_K \lesssim 8$ ms.

Regardless of this degradation of correlation and temporal overlap between interferometers, the differential FRAC method is able to restore the interference fringes with a good SNR (~30 for $^{87}$Rb, ~10 for $^{39}$K, limited by uncorrelated offset noise). The vibration rejection factor for each of the data sets shown in figure 7 is approximately $\gamma \sim 100$. This permits unbiased estimates of $\phi_d$ with a statistical uncertainty at the level of $\delta\phi_d \sim 25$ mrad with 300 points—corresponding to WEP test with a statistical sensitivity of $\delta\eta \approx \delta\phi_d / k_{\text{Rb}}^{\text{eff}} g T_{\text{Rb}}^{-2} = 1.6 \times 10^{-6}$ per data set. The robustness of the differential FRAC technique under ‘noisy’
<figure><img src="image_22.png" /><figcaption>**Figure 8.** Comparison of differential phase estimates from the Bayesian and FRAC methods for various $\phi_d$ and $\kappa$. (a) Difference between Bayesian and FRAC estimates $\phi_d^{\text{Bayes}} - \phi_d^{\text{FRAC}}$ as a function $\phi_d^{\text{FRAC}} = \phi_K - \kappa\phi_{\text{Rb}}$. The vertical error bars represent the combined statistical uncertainty of both estimates. (b) Difference between Bayesian and FRAC estimates as a function of $\kappa = S_K/S_{\text{Rb}}$ corresponding to each point in (a). The majority of the data points shown in both figures are consistent with zero to within $1\sigma$ of uncertainty, indicating agreement between the two techniques at the level of ~40 mrad. The standard deviation of offset noise for each data set was typically $\sigma_{B_K} \approx 0.20$ and $\sigma_{B_{\text{Rb}}} \approx 0.05$ in the normalized space ($n_1$ and $n_2$). A value of $\sigma_{\phi_d} = 0.05$ rad was used for the differential phase noise of all data sets. The range of common phase noise was estimated from accelerometer measurements.</figcaption></figure>

conditions makes it an ideal candidate for future WEP tests [7, 53, 54], or other differential atom interferometry applications [34, 40].

In contrast, for Bayesian estimation, an increase in uncorrelated phase noise is problematic. When the 'common phase' becomes largely uncorrelated, the Bayesian method can converge on multiple possible $\phi_d$, or may not converge at all. For these data, we find that by $T_K = 9$ ms the Bayesian estimate of $\phi_d$ is not consistent with the FRAC estimate, and for $T_K \lesssim 8$ ms the analysis is not able to converge on a unique value. We note that these particular results are strongly dependent on the level of phase noise, the degree of temporal overlap, the value of $\phi_d$ and the scale factor of each interferometer. In the following section, we study some of these dependencies more quantitatively.

## 5.2. Comparison of Bayesian and FRAC methods as a function of $\kappa$ and $\phi_d$

We have tested the functionality and accuracy of both the Bayesian and FRAC methods for extracting $\phi_d$ from experimental data acquired under various conditions. Specifically, we are interested in the accuracy of these techniques over (i) the full range of differential phase $\phi_d \in [0, \pi]$, and (ii) a broad range of interferometer scale factor ratios $\kappa = S_K/S_{\text{Rb}}$. To investigate these two aspects, we recorded data using the symmetric interferometer configuration shown in figure 4 with different interrogation times, $T_{\text{Rb}}$ and $T_K$. Since $\kappa$ is proportional to $(T_K/T_{\text{Rb}})^2$, each configuration of $T_j$ corresponds to a different scale factor ratio. Additionally, the differential phase is modified with each $T_K$ due to a systematic phase shift of the potassium interferometer from an external magnetic field. Therefore, we are able to study both effects with a single data set.

Figure 8 shows a comparison between Bayesian and FRAC estimates of $\phi_d$, using the FRAC estimate as a reference. We varied $T_{\text{Rb}}$ from 1 to 5 ms, and $T_K$ independently in the vicinity of $T_{\text{Rb}}$ such that the scale factor ratio was modified over a relatively broad range ($\kappa \approx 0.45$ to 1.01). The phase noise due to the externally applied vibrations was kept quasi-common-mode between sensors by ensuring that $T_K$ was within a few $100~\mu\text{s}$ of $T_{\text{Rb}}$. Over this range of $T_{\text{Rb}}$ and $T_K$, we found that the differential phase ranged from roughly $\phi_d = 0$ to 2.8 rad as a result of a systematic shift of the potassium interferometer. It is clear from figure 8 that there is a high degree of correlation between the Bayesian and FRAC estimates, which is consistent with our expectations based on the simulations discussed in section 3.2. The error bars in this figure were computed from the combined statistical uncertainties of both methods, which both typically yield $\delta\phi_d \sim 30$ mrad at the present level of noise.

To summarize, we find that the difference between the two estimates is consistent with zero within a typical total uncertainty of ~40 mrad. These data confirm that the two analysis techniques produce unbiased estimates of $\phi_d$ for dual-species interferometers with vastly different scale factors. We discuss further the advantages and limitations of these two techniques in the following section.

# 6. Advantages and limitations of the methods

As discussed in section 3.2, Bayes' method is optimally efficient and yields a statistical error that scales as $1/\sqrt{N}$, compared to more heuristic fitting techniques which converge more slowly. This improved efficiency is a clear advantage of the Bayesian estimator compared to the FRAC analysis. However, the disadvantage is that it requires *a priori* information about the system, such as noise levels and interferometer contrasts, and it requires
<table><caption>**Table 1.** Comparison between phase noise for a single interferometer, and two coupled interferometers with effective wave vectors $k_{Rb}^{\text{eff}}$ and $k_K^{\text{eff}}$. The rms phase noise (in radians) due to vibrations ($\phi_{rms}^{\text{vib}}$) and the self-noise of the mechanical accelerometer ($\phi_{rms}^{\text{self}}$) are shown for different frequency bands and interrogation times, $T_{Rb}$. The noise from each band is summed in quadrature to obtain the total noise. Contributions less than 1 mrad are not shown. For the simultaneous differential sensor, it is assumed that $k_{Rb}^{\text{eff}} T_{Rb}^2 = k_K^{\text{eff}} T_K^2$. The rms phase noise $\phi_{rms}^{\text{vib}}$ was computed from equation (C.9) using model (C.13) for the power spectral density $S_a(\omega)$ of ground accelerations in a 'quiet' location [66, 67] with integrated rms noise $1.4 \times 10^{-4} \text{ g}$ (see figure C2(b)). The quantity $\phi_{rms}^{\text{self}}$ was computed in a similar manner by replacing $S_a(\omega)$ with the self-noise spectrum of the accelerometer—here assumed to be white noise with $|S_a|^{1/2} \approx 3.2 \times 10^{-8} \text{ g}/\sqrt{\text{Hz}}$. Quantities in the last row indicated by '*' correspond to low-noise conditions that can be achieved with passive vibration isolation (integrated rms noise $1.4 \times 10^{-6} \text{ g}$), and an accelerometer with 10 times smaller self-noise of $|S_a|^{1/2} \approx 3.2 \times 10^{-9} \text{ g}/\sqrt{\text{Hz}}$.</caption><thead><tr><th rowspan="2">T<sub>Rb</sub> (s)</th><th rowspan="2">Noise Source</th><th colspan="4">Single sensor</th><th colspan="4">Differential sensor</th></tr><tr><th>dc − 1</th><th colspan="2">Frequency range (Hz)</th><th>dc − ∞</th><th>dc − 1</th><th colspan="2">Frequency range (Hz)</th><th>dc − ∞</th></tr></thead><tbody><tr><td rowspan="2">0.01</td><td>(&phi;<sub>rms</sub><sup>vib</sup>)</td><td>0.023</td><td>0.779</td><td>1.071</td><td>1.324</td><td>0.000</td><td>0.000</td><td>0.004</td><td>0.004</td></tr><tr><td>(&phi;<sub>rms</sub><sup>self</sup>)</td><td>0.000</td><td>0.002</td><td>0.002</td><td>0.003</td><td>0.000</td><td>0.002</td><td>0.002</td><td>0.003</td></tr><tr><td rowspan="2">0.1</td><td>(&phi;<sub>rms</sub><sup>vib</sup>)</td><td>2.240</td><td>19.62</td><td>2.962</td><td>19.97</td><td>0.001</td><td>0.380</td><td>0.140</td><td>0.405</td></tr><tr><td>(&phi;<sub>rms</sub><sup>self</sup>)</td><td>0.050</td><td>0.078</td><td>0.005</td><td>0.093</td><td>0.050</td><td>0.078</td><td>0.005</td><td>0.093</td></tr><tr><td rowspan="4">1</td><td>(&phi;<sub>rms</sub><sup>vib</sup>)</td><td>58.13</td><td>23.42</td><td>3.779</td><td>62.79</td><td>1.067</td><td>2.556</td><td>1.511</td><td>3.155</td></tr><tr><td>(&phi;<sub>rms</sub><sup>self</sup>)</td><td>2.936</td><td>0.160</td><td>0.006</td><td>2.940</td><td>2.936</td><td>0.160</td><td>0.006</td><td>2.940</td></tr><tr><td>(&phi;<sub>rms</sub><sup>vib</sup>)<sup>*</sup></td><td>28.06</td><td>1.868</td><td>0.021</td><td>28.12</td><td>0.303</td><td>0.098</td><td>0.009</td><td>0.319</td></tr><tr><td>(&phi;<sub>rms</sub><sup>self</sup>)<sup>*</sup></td><td>0.294</td><td>0.016</td><td>0.001</td><td>0.294</td><td>0.294</td><td>0.016</td><td>0.001</td><td>0.294</td></tr></tbody></table>

significant computational resources to evaluate. Furthermore, it is only a viable solution for simultaneous interferometer configurations that exhibit a high-degree of phase correlation.

In contrast to the Bayesian estimator, the FRAC method requires only the interferometer timing parameters and a sensitive accelerometer that is well-coupled to the reference mirror in order to function accurately. It does not assume any particular interferometer configuration or require any additional system information. The FRAC method also has applications in absolute interferometry, as has been previously demonstrated in [26, 66–69]. Additionally, it is fast enough to be used for real-time feedback, which has been shown to improve single-sensor sensitivity [70].

Table 1 contains estimates of the phase noise for a single $^{87}$Rb interferometer, and two coupled interferometers of $^{87}$Rb and $^{39}$K. The rms phase spread due to vibration noise ($\phi_{rms}^{\text{vib}}$) and the self-noise of the accelerometer ($\phi_{rms}^{\text{self}}$) are shown for various frequency bands and interrogation times. For a single sensor analyzed with the FRAC method, the vibration-induced noise $\phi_{rms}^{\text{vib}}$ represents the spread of phase on the uncorrected fringes, while the quantity $\phi_{rms}^{\text{self}}$ indicates the residual phase noise present on the corrected fringes. Since this term is directly linked to the intrinsic noise of the mechanical accelerometer, it represents a fundamental limitation of the standard FRAC method. To give a quantitative example, based in the self-noise of the Titan accelerometer used in our experiments ($3.2 \times 10^{-8} \text{ g}/\sqrt{\text{Hz}}$), the corresponding phase noise reaches ~90 mrad for an interrogation time of 100 ms, and ~3 rad by $T_{Rb} = 1$ s. With this level of self-noise, fringes cannot be reconstructed accurately. However, for a state-of-the-art device with an order of magnitude smaller self-noise ($3.2 \times 10^{-9} \text{ g}/\sqrt{\text{Hz}}$), the phase noise decreases by a further factor of 10—allowing a least-squares fit to accurately converge on the fringe phase. We also point out that the noise contributions from both vibrations and self-noise are smallest at high frequencies—a result of the natural low-pass filtering of atom interferometers. Thus, high-bandwidth accelerometers are generally not required to implement the FRAC method with a single sensor. For instance, for ground-based WEP test facilities targeting a few $10^{-15}$ [43, 44, 46] using interrogation times of $T \sim 1$ s, we estimate that a mechanical accelerometer with a self-noise less than $10^{-11} \text{ g}/\sqrt{\text{Hz}}$ in the DC to 1 Hz frequency band will yield a phase noise contribution below the projected shot-noise limit of ~1 mrad for each interferometer.

When employing the differential FRAC method with two simultaneous interferometers, the self-noise of the accelerometer contributes to the phase of both sensors. Thus, in table 1 (where $\kappa = 1$ and $T_K = T_{Rb} \sqrt{k_{Rb}^{\text{eff}}/k_K^{\text{eff}}}$), we have indicated the same values for $\phi_{rms}^{\text{self}}$ in the corresponding columns for both single and differential sensors. However, we emphasize that one can measure the differential phase significantly more accurately than the self-noise limit for a single sensor. This is because the noise introduced by the accelerometer is correlated between the two interferometers—reducing the uncertainty in the determination of $\phi_d$, as discussed in section 3.3. A recent study [71] has shown that uncertainties close to the quantum-projection-noise limit can be obtained with this method when the interferometers are in phase ($\phi_d = 0$) and the accelerometer exhibits a conservative level of self-noise ($\phi_{rms}^{\text{self}} \lesssim 0.3$ rad). In general, for the differential FRAC method to function well for all values of $\phi_d$
and $\kappa$, the self-noise of the accelerometer should correspond to less than $\pi/2$ in phase noise for each sensor— allowing individual fringes to be accurately fit. However, for the special case of $\phi_d = 0$ and $\kappa = 1$, the requirements on the accelerometer noise are much less stringent. For instance, [71] indicates that reliable fits can be obtained with up to $\phi_{rms}^{\text{self}} \sim 20$ rad. For these reasons, we emphasize that state-of-the-art mechanical accelerometers are not required to make sensitive measurements of $\phi_d$ with long-baseline differential interferometers. We anticipate that competitive levels of accuracy can be achieved with compact devices that feature a moderate level of sensitivity.

For two coupled interferometers exhibiting different wave vectors, the vibration-induced phase noise is not identical and thus cannot be perfectly rejected at all frequencies. The values of $\phi_{rms}^{\text{vib}}$ listed in the last four columns of table 1 contribute directly to $\phi_d$—representing the level of uncorrelated differential phase noise in the system. We estimate that by $T_{\text{Rb}} = 1$ s the differential phase noise reaches a level of $\sim 3$ rad. However, we note that the differential transfer function (equation (C.11)) rejects most efficiently at frequencies below $\sim 1/T$, and this estimate is directly linked to the vibration spectrum used. In a quieter environment, such as that achieved with a vibration isolation platform [66, 67] or in a satellite [52], the phase noise can be reduced by an order of magnitude or more. At this point, the Bayesian method can be employed—which easily handles differential phase noise. Since the sensitivity scales as $\phi_{rms}^{\text{vib}}/\sqrt{N}$, the analysis simply requires more measurements for larger $\phi_{rms}^{\text{vib}}$ to reach a given level of precision.

## 7. Conclusion

We have described and demonstrated experimentally two new analysis techniques for extracting the differential phase from coupled atom interferometers with different scale factors, $S_j$. A non-unity ratio $\kappa = S_1/S_2$ can result from using atoms with different $k_j^{\text{eff}}$, or from interferometers with different interrogation times, $T_j$. We also carried out correlated phase measurements between simultaneous interferometers of two chemical elements exhibiting different scale factors, and we have demonstrated a vibration rejection factor as large as $\gamma \approx 730$. This system was used to validate the Bayesian and FRAC analysis methods, as well as a new ellipse fitting procedure [59], for extracting $\phi_d$. Furthermore, the FRAC method was used to demonstrate a statistical sensitivity for the Eötvös parameter of $\delta\eta = 1.6 \times 10^{-6}$ per measurement with $T = 10$ ms interferometers operating in a ‘noisy’ environment.

Presently, the contrast and SNR of our $^{39}$K interferometer fringes are limited by, respectively, the temperature of the atoms and technical noise present in the slave laser used for cooling, interferometry and detection. The latter issue results in uncorrelated offset noise, which reduces the correlation between interferometers. This, in turn, increases the error on $\phi_d$, which degrades the vibration rejection factor $\gamma$ and the sensitivity $\delta\eta$. We anticipate that a modest reduction of both the technical noise, and further cooling the $^{39}$K sample in a gray molasses [79, 80], will result in a substantial improvement in the correlation between rubidium and potassium. A precise determination of $\eta$ with our apparatus, including a complete evaluation of systematic effects, is beyond the scope of this work, but will be the subject of a future publication.

Both the generalized Bayesian and differential FRAC methods yield unbiased estimates of $\phi_d$ for any scale factor ratio, $\kappa$, and are robust against experimental parameters such as the common phase range scanned by the two interferometers, or the level of uncorrelated offset noise present in the system. These features make both methods ideal for applications of dual-species interferometry where, until now, the available analysis tools could accommodate only systems that exhibit either $\kappa = 1$ or low levels of common phase noise. These new methods are also appealing for gradiometer configurations using the same atoms and the same $T_j$ [71], which have previously been utilized for precisely measuring the gravitational constant $G$ and gravity gradients [27–29, 31–33, 61].

The freedom to vary the scale factor, the interrogation time or phase of either interferometer independently can be advantageous for studying systematic effects, interactions between atomic species [55], or for shifting the differential phase toward a region of higher sensitivity. Examples of such regions include $\phi_d = \pi/2$ in the case of ellipse-fitting methods, and $\phi_d = 0$ or $\pi$ for the FRAC technique [71]. Both the FRAC and Bayesian methods also eliminate the systematic shift introduced on the measurement of $\Delta a$ when using dual-species interferometers with $\kappa \neq 1$—making them well-suited for upcoming WEP tests on ground [43–46], in microgravity [41, 47–50], and in Space [7, 51–54].

## Acknowledgments

This work is supported by the French national agencies CNES (Centre National d'Etudes Spatiales), l'Agence Nationale pour la Recherche, the Délégation Générale de l'Armement, the European Space Agency, IFRAF
(Institut Francilien de Recherche sur les Atomes Froids), action spécifique GRAM (Gravitation, Relativité Astronomie et Métrologie) and RTRA 'Triangle de la Physique'. B Barrett and L Antoni-Micollier thank CNES and IOGS for financial support. P Bouyer thanks Conseil Régional d'Aquitaine for the Excellence Chair. Finally, the ICE team would like to thank the following people: F Pereira Dos Santos of the laboratory SYRTE and V Ménoret of MuQuans for helpful discussions; D Holleville, B Venon, F Cornu of SYRTE and J-P Aoustin of the laboratory GEPI for their technical assistance building vacuum and optical components, and T Rvachov of the Massachusetts Institute of Technology for his assistance with the 'Cicero Word Generator' control software.

## Appendix A. Ellipse fitting methods

In this appendix, we give some background regarding ellipse-fitting techniques and illustrate the problem of parameter bias for two different fitting algorithms.

The general form of an ellipse in a cartesian plane is described by the algebraic equation for a conic

$$F(\lambda, y) = \lambda \cdot y = A y_1^2 + B y_1 y_2 + C y_2^2 + D y_1 + E y_2 + F = 0, \quad (A.1)$$

provided that $B^2 < 4AC$. Here, $\lambda = \{A, B, C, D, E, F\}$ and $y = \{y_1^2, y_1 y_2, y_2^2, y_1, y_2, 1\}$. The center, orientation, major and minor axes of the ellipse are determined by the elements of $\lambda$, and the differential phase can be shown to be

$$\phi_d = \cos^{-1}\left(-\frac{B}{2\sqrt{AC}}\right). \quad (A.2)$$

Generally, two types of ellipse-fitting algorithms exist: those that seek to minimize (i) an algebraic distance or (ii) a geometric/orthogonal distance between the ellipse and the data points. While algebraic methods tend to be simple, efficient and can guarantee an ellipse solution to the conic equation (A.1) (i.e. parabolic and hyperbolic solutions can be eliminated), they tend to suffer highly from bias in the ellipse parameters—resulting in a poor fit under certain circumstances. Geometric methods are usually much more accurate than algebraic algorithms, but at the cost of more complexity, more computation and less stability. Since minimizing the orthogonal distance between a point and an ellipse has no closed-form solution, these routines resort to iterative techniques that are not guaranteed to converge on an ellipse.

A commonly used algebraic method is the simple and robust DEF method developed by Fitzgibbon et al [62] that minimizes the sum of squared algebraic distances between the points and the ellipse, $\sum_{i=1}^N F(\lambda, y_i)^2$, subject to the constraint $B^2 = 4AC - 1$. Recently, Szpak et al [59, 63] developed an algorithm based on the optimization of the approximate maximum likelihood distance which seeks a balance between the costly geometric methods and stable algebraic techniques. This algorithm—termed the FGEF method—also includes error estimation for the geometrically meaningful ellipse parameters (center coordinates, axes and orientation) which we have extended to include an estimate of the differential phase error, $\delta\phi_d$.

Figure A1 illustrates the bias introduced on the differential phase estimated by the DEF and FGEF methods. For moderate amounts of noise in the offset, the DEF method tends to produce fits that are characteristically compressed along the major axis and stretched along the minor axis of the ellipse, as shown by the red curve in figure A1(a). This effect results in a biased estimate of $\phi_d$ that increases monotonically away from $\pi/2$, as shown in figure A1(b). In contrast to the DEF method, the FGEF algorithm predicts an ellipse (shown in blue) that is much more representative of the actual ellipse (shown in green), and also results in less bias in $\phi_d$ in the central region around $\pi/2$. Outside of this region, the bias behaves nonlinearly in a manner that depends on the ellipse parameters and the level of noise. Here, we point out that these bias estimates are dependent on the type of noise (offset, amplitude, or differential phase) and the amount of noise present in the data, but typically the bias is smallest in the vicinity of $\phi_d = \pi/2$, and decreases with the noise level. In general, ellipse-fitting techniques always generate a non-zero systematic on the differential phase estimate, and depending on the level of sensitivity, this bias must be carefully accounted for when performing precise measurements with $\phi_d$ [28, 31, 32, 34, 61].

## Appendix B. Bayesian analysis of Lissajous curves

In this appendix, we describe in detail our generalized Bayesian analysis technique to estimate the differential phase from Lissajous curves. We also demonstrate the effectiveness of this method using numerically simulated data with Gaussian noise in the offset parameters $\{B_1, B_2\}$ and the differential phase, $\phi_d$. Noise in the amplitude parameters $\{A_1, A_2\}$ of the coupled-sensor model (2) can also be included via a trivial modification of the noise model. In what follows, we first provide some relevant theoretical background of the Bayesian estimation technique. For a more comprehensive description of Bayesian analysis in this context, see [65].
<figure><img src="image_23.png" /><figcaption>**Figure A1.** (a) Synthetic data following an ellipse with added offset noise. The solid green curve represents the actual ellipse, and fits to the data using the DEF method (red curve with big dashes) and FGEEF method (blue curve with small dashes). The simulated ellipse contains 500 points with Gaussian-distributed noise on the offset parameters $B_j$ with standard deviations {$\sigma_{B_1}, \sigma_{B_2}$} = {0.01, 0.03} (corresponding to SNR ~{20, 6}). Ellipse parameters: $A_1 = A_2 = 0.2, B_1 = B_2 = 0.5, \kappa = 1, \phi_d = 1$ rad. (b) Measured bias in differential phase estimates, $\phi_d^{est}$, from the DEF (red triangles) and FGEEF (blue points) methods relative to the actual value, $\phi_d^{act}$. The black squares show the estimates from the differential FRAC method for comparison. On all plots, the error bars correspond to the statistical distribution of fits to 100 synthetic data sets.</figcaption></figure>

In a generalized system, where *M* represents a measurement of the system quantities and *V* represents a variable we are interested in measuring, Bayes' rule can be summarized by the following equation

$$P(V|M) = \frac{p(V)L(M|V)}{N(M)}. \qquad (B.1)$$

Here, $P(V|M)$ is called the 'posterior' probability distribution and represents our state of knowledge after a measurement, $M$. $p(V)$ is the 'prior' probability before the measurement, and $L(M|V)$ is called the 'likelihood' to obtain a certain result for *M* given *V*. The key to the entire estimation process is the likelihood distribution, which is computed based on a specific model of the noise present in the system. The quantity $N(M) = \sum_V L(M|V)p(V)$ is the probability of measuring *M* integrated over all possible values of *V*, and is just a normalizing factor for the posterior distribution. Mathematically, $L(M|V)$ can be thought of as a function of *V* with *M* fixed, and vice versa for $P(V|M)$. The essence of Bayes' rule is that knowledge of the variable *V* can be updated on a measurement-by-measurement basis—with each successive measurement contributing additional information that narrows the width of the probability distribution associated with *V*. A well-known example of this type of recursive analysis is a Kalman filter [81], which is used extensively in the fields of guidance, navigation and trajectory optimization.

For the specific case of two coupled atom interferometers, the variable of interest is $\phi_d$ and the ith system measurement is given by the pair of (normalized) atomic state populations $M_i = \{n_1, n_2\}_{i}$. Thus, for a single measurement equation (B.1) becomes

$$P(\phi_d | \{n_1, n_2\}_{i})_i = \frac{p(\phi_d)_i L(\{n_1, n_2\}_{i} | \phi_d)_i}{N(\{n_1, n_2\}_{i})_i}, \qquad (B.2)$$

where $P(\phi_d|\{n_1, n_2\}_{i})$ is referred to as the conditional distribution based on the ith measurement. The basic algorithm for Bayes' estimation can be summarized as follows:

(1) Choose a suitable initial prior distribution, $p(\phi_d)_{i=1}$. In our case, we take this to be a uniform distribution within the range $\phi_d \in [0, \pi]$, and zero elsewhere.

(2) Record a new measurement $\{n_1, n_2\}_{i}$, and calculate the likelihood distribution $L(\{n_1, n_2\}_{i}|\phi_d)_{i}$ from the noise model.

(3) Compute the conditional probability distribution $P(\phi_d|\{n_1, n_2\}_{i})$ from Bayes' rule (B.2).

(4) Set the new prior distribution equal to the previous conditional distribution: $p(\phi_d)_{i+1} = P(\phi_d|\{n_1, n_2\}_{i})$.

(5) Repeat steps (2) through (4) until the width of the conditional distribution reduces to the desired level.
<figure><img src="image_24.png" /><figcaption>**Figure B1.** Example of a Lissajous curve for $\kappa = p/q = 3/7$. (a) The curve is plotted over the phase range $\phi_c \in [0, 2\pi]$ and the solutions $\{n_1, n_2 = s_{2,\ell}\}$ are shown as points. For $n_1 = -0.9$ there are two solutions $s_{2,\ell}$, shown as the blue points. Similarly, for $n_1 = 0$ there is only one possible value of $n_2$ (shown in green), and for $n_1 = 0.8$ no solutions exist. (b) The same Lissajous curve plotted over $\phi_c \in [0, 2\pi q]$ with $q = 7$. In this range, there are always 6 solutions $s_{2,\ell}$ for each value of $n_1$ (although they may not be unique).</figcaption></figure>

(6) Estimate the variable of interest, $\phi_d$, using the maximum likelihood value of the final conditional (i.e. posterior) probability distribution.

This algorithm is illustrated in figure 1 in section 3.2.

### The likelihood distribution

The main challenge in Bayesian analysis is to compute the likelihood distribution $L(\{n_1, n_2\} | \phi_d)$ given a specific model for $n_1$ and $n_2$. For the specific case of coupled interferometers, there are three possible sources of noise: amplitude, offset and differential phase. To illustrate each source, we modify the definitions of the $n_j$ in equation (10) to explicitly include these noise terms

$$n_1(\phi_c) = (1 + \delta A_1) \cos(\kappa \phi_c + \phi_d + \delta \phi_d) + \delta B_1, \qquad (B.3a)$$

$$n_2(\phi_c) = (1 + \delta A_2) \cos(\phi_c) + \delta B_2. \qquad (B.3b)$$

The parameters $\delta A_j$, $\delta B_j$, and $\delta \phi_d$ represent uncorrelated noise in the amplitude, offset and differential phase, respectively, each of which is assumed to follow a Gaussian probability distribution with zero mean and non-zero standard deviation. Using this model, the likelihood distribution can be shown to be [65]

$$L(\{n_1, n_2\} | \phi_d) = \sum_{\ell} \int_{-1}^{1} \frac{P(n_1|s_1) P(n_2|\{s_{2,\ell}; s_1, \phi_d\})}{\sqrt{1-s_1^2}} ds_1. \qquad (B.4)$$

Here, $P(n_1|s_1)$ and $P(n_2|\{s_{2,\ell}; s_1, \phi_d\})$ are the single-sensor conditional probability distributions for $n_1$ and $n_2$, which we discuss in more detail below. The quantities $s_1 \equiv \cos(\kappa\phi_c + \phi_d)$ and $s_{2,\ell} \equiv \cos(\phi_c)$ are the principle variables on which the coupled measurements $n_1$ and $n_2$ depend in the model (B.3). Due to the periodic nature of the Lissajous equations (10), for each value of $n_1$ there are multiple possible solutions for $n_2$ (as shown in figure B1). We assign an integer $\ell$ to each of these solutions. More specifically, $s_{2,\ell}$ is the $\ell$th root of $n_2$ given $n_1 = s_1$. The sum over $\ell$ appearing in equation (B.4) accounts for all possible solutions. In the distribution functions $P(n_1|s_1)$ and $P(n_2|\{s_{2,\ell}; s_1, \phi_d\})$, we denote the implicit dependence on variables $s_1$ and $\phi_d$ by a semicolon. This notation emphasizes that the quantity $s_{2,\ell}$ is coupled to $s_1$ through the common phase$^6 \phi_c$. Finally, we point out that the coupled variables $s_1$ and $s_{2,\ell}$ both depend on $\phi_d$, but we do not write this dependence explicitly.

At this point, we need to know the possible values $n_2 = s_{2,\ell}$ (given a measurement of $n_1 = s_1$) which enter into the likelihood distribution. We devote the remainder of this section to a detailed description of computing the roots of the Lissajous equations (10). As mentioned above, due to the nonlinear nature of Lissajous curves, there are multiple possible solutions for $n_2$ given a single value of $n_1$ within a predefined phase range. We denote these solutions $s_{2,\ell}$ for integer $\ell$. When $\kappa = 1$, the Lissajous curve collapses to an ellipse, and only two values of $n_2$ exist for each $n_1$ over any $2\pi$ range of $\phi_c$. In this case, it is straightforward to compute the two solutions as $s_{2,\pm 1} = \cos[\cos^{-1}(s_1) \pm \phi_d]$. However, when $\kappa \neq 1$, the problem is much more complex. If the scale factor ratio

<footnote>$^6$Since it is assumed that $\phi_c$ is random and unknown, the probability distributions of $s_1$ and $s_{2,\ell}$ are equivalent to that of a sinusoid: $P(s_{2,\ell}|\phi_c) = (1 - s_{2,\ell}^2)^{-1/2}$.</footnote>
can be written in the form $\kappa = p/q$, where p and q are prime numbers, then the period of the Lissajous curve is $2\pi q$—requiring q revolutions to form a closed loop. Within each $2\pi$ interval, there can be either 0, 1 or 2 solutions of $n_2$ for each $n_1$, as illustrated in figure B1.

To calculate these solutions for a given $n_1 = s_1$ and $\phi_d$, it is necessary to know the approximate range of common phase spanned by the data<fnref n="7" />: $\phi_c \in [\phi_c^{\min}, \phi_c^{\max}]$. With this information, we compute the range of phase spanned by sensor 1, $\theta \in \kappa [\phi_c^{\min}, \phi_c^{\max}] + \phi_d$, and we subdivide this range into intervals of $\pi$ such that the $\ell$th interval is defined as the range $\theta_\ell \in [\ell, \ell + 1)\pi$, where $\ell = [\theta/\pi]$. Here, the brackets $[\cdots]$ indicate the floor function. Beginning with the left-most interval, we check for solutions sequentially at each $\pi$ phase bin until the entire range is spanned. Empirically, we find that if a solution exists within the $\ell$th interval given a value $n_1$, then it is unique and can be written explicitly as

$$s_{2,\ell} = \begin{cases} \cos\left[\left(\cos^{-1}s_1 - \phi_d + 2\pi m_{1,\ell}\right)/\kappa\right] & \text{for even } \ell, \\ \cos\left[\left(\cos^{-1}s_1 + \phi_d - 2\pi m_{2,\ell}\right)/\kappa\right] & \text{for odd } \ell, \end{cases} \quad (B.5)$$

where the integers $m_{1,\ell}$ and $m_{2,\ell}$ are defined as

$$m_{1,\ell} = \begin{cases} \lfloor (\ell - 1)/2 \rfloor & \ell < -1, \\ 0 & -1 \le \ell < 2, \\ \lfloor \ell/2 \rfloor & \ell \ge 2, \end{cases} \quad m_{2,\ell} = \begin{cases} \lfloor \ell/2 \rfloor & \ell < -2, \\ 0 & -2 \le \ell < 1, \\ \lfloor (\ell + 1)/2 \rfloor & \ell \ge 1. \end{cases} \quad (B.6)$$

With these solutions in hand, it is possible to compute the likelihood (B.4) given specific noise models for the single-sensor probability distributions $P(n_1|s_1)$ and $P(n_2|\{s_{2,\ell}; s_1, \phi_d\})$. We now investigate the specific cases of offset and differential phase noise on the extraction of $\phi_d$ from simulated data sets. This analysis can also be extended to include noise in the fringe amplitudes through the parameters $\delta A_j$ [65], but we do not consider this case here.

### Offset noise

When the system exhibits noise only in the offset of the atomic state measurements, the parameters $\delta B_j$ are randomly distributed for each repetition of the experiment, and $\delta A_j = 0$ and $\delta \phi_d = 0$ in the model (B.3). Under realistic conditions, these noise parameters follow a Gaussian distribution with zero mean and standard deviations given by $\sigma_{B_j}$, and the single-sensor conditional probabilities can be written as

$$P(n_1|s_1) \propto \exp\left[-(n_1 - s_1)^2 / 2\sigma_{B_1}^2\right], \quad (B.7a)$$

$$P(n_2|\{s_{2,\ell}; s_1, \phi_d\}) \propto \exp\left[-(n_2 - s_{2,\ell})^2 / 2\sigma_{B_2}^2\right]. \quad (B.7b)$$

Figure B2 shows some examples of simulated data in the presence of offset noise, where the differential phase has been extracted using the Bayesian estimation algorithm described above. These simulations show that $\phi_d$ can be precisely estimated over the full range of $0 - \pi$, and for a wide variety of scale factor ratios. Here, we demonstrate the technique for the limited range $\kappa \in [0.6, 1.4]$, but we have also verified that the extraction method works well outside this range. In contrast to ellipse-fitting techniques, no systematic bias in the phase estimates is observed, and fewer points are required to converge to competitive error levels.

### Differential phase noise

Since the noise parameter associated with the differential phase, $\delta\phi_d$, adds directly to the quantity of interest, $\phi_d$, we can account for this type of noise by adding an extra convolution with our noise model at the end of any likelihood calculation. We choose to examine the case of Gaussian noise for the differential phase, such that the conditional probability distribution is

$$P(\phi'_d | \phi_d) \propto \exp\left[-(\phi'_d - \phi_d)^2 / 2\sigma_{\phi_d}^2\right]. \quad (B.8)$$

Here, $\phi'_d$ represents a measured value of the differential phase in the presence of Gaussian noise centered on the most likely value, $\phi_d$, and $\sigma_{\phi_d}$ is the standard deviation of the noise distribution. The modified likelihood function is described by the convolution

$$L(\{n_1, n_2\} | \phi_d) \propto \int_{-\infty}^{\infty} L(\{n_1, n_2\} | \phi'_d) P(\phi'_d | \phi_d) d\phi'_d. \quad (B.9)$$

In a similar fashion to the offset, in the absence of any other noise sources it is necessary to estimate multiple candidate solutions for $\phi_d$ over a given range of $\phi_c$ in order to compute the likelihood function. Before convolving with the conditional probability distribution in equation (B.9), the likelihood function can be

<footnote n="7">⁷ In practice, this range estimate does need not to be very precise—we find that estimating the correct range to within $\pm\pi$ still results in a precise estimate for $\phi_d$. However, overestimating the phase range may result in a slower convergence rate for the estimate. See section 5.1 for a description of how the phase range can be estimated experimentally.</footnote>
<figure><img src="image_25.png" /><figcaption>Figure B2. Simulated Bayesian phase estimates in the presence of offset noise. Top row: simulated data for $\kappa = \{0.6, 0.8, 1.0, 1.2, 1.4\}$. Here, points based on equation (B.3) were generated from common phases chosen randomly over the range $\phi_c \in [-2\pi, 2\pi]$. The differential phase was set to $\phi_d^{\text{act}} = 1$ rad, and offset noise was applied to the points via Gaussian distributions with standard deviations $\{\sigma_{B_1}, \sigma_{B_2}\} = \{0.02, 0.04\}$—corresponding to SNR ~50 and 25, respectively. Bottom row: the systematic error in the Bayesian phase estimate, $\phi_d^{\text{est}} - \phi_d^{\text{act}}$, as a function of $\phi_d^{\text{act}}$ for each value of $\kappa$. Only 10 points were required to bring the statistical uncertainty indicated by the error bars to $\le 50$ mrad.</figcaption></figure>

written as

$$L(\{n_1, n_2\} | \phi_d') = \sum_k \delta(\phi_d' - \phi_{d,k}), \qquad (B.10)$$

where $\delta(x)$ is the Dirac delta function, and the sum over $k$ accounts for all candidate solutions $\phi_{d,k}$ that exist in the common phase range $\phi_c \in [\phi_c^{\text{min}}, \phi_c^{\text{max}}]$. These solutions can be computed by, again, dividing the phase range into intervals of $\pi$, and labeling each of them by an integer $k = |\phi_c/\pi|$. We find that two possible solutions exist for $\phi_d$ within each interval, which we denote as $\phi_{d,k}^{(\pm)}$ for $\phi_c \in [\pm k\pi, \pm(k+1)\pi)$. Explicitly, these phases can be computed from

$$\phi_{d,k}^{(\pm)} = \cos^{-1}(n_1) \pm \kappa[\cos^{-1}(n_2) - 2\pi m_k], \qquad (B.11)$$

where $m_k = (-1)^k [(|k| + 1)/2]$. We transform these phases into the range of $0 - \pi$ using $\phi_{d,k}^{(\pm)} \rightarrow \cos^{-1}(\cos(\phi_{d,k}^{(\pm)}))$. Finally, this result is convolved with the Gaussian noise model to obtain

$$L(\{n_1, n_2\} | \phi_d) = \sum_k \exp\left[-(\phi_d - \phi_{d,k}^{(\pm)})^2 / 2\sigma_{\phi_d}^2\right]. \qquad (B.12)$$

Two subtleties exist with this analysis, however, that warrant discussion. First, when the common phase range exceeds $\phi_c \in [-\pi, \pi]$, the Bayesian analysis may predict multiple equally probable values for $\phi_d$. This is obviously a problem if we are interested in a precise, unique estimate of the differential phase, and we have no pre-existing knowledge of its value. Therefore, we restrict our consideration of the problem to a range of common phase within $-\pi$ to $\pi$. Second, the noise parameter $\delta\phi_d$ can theoretically take any value, i.e. $\delta\phi_d \in (-\infty, \infty)$, although in practice it is limited to a finite range defined by $\sigma_{\phi_d}$. So far, we have considered $\phi_d$ only in the range of $0$ to $\pi$, but for situations where $\sigma_{\phi_d} \ge \pi/4$, the likelihood distribution can have significant contributions from the wings of the adjacent $\pi$ phase intervals. This effect can be taken into account by using the fact that $P(\phi_d) = P(-\phi_d) = P(2\pi - \phi_d)$, and adding mirrored versions of the likelihood to the convolution in equation (B.9). This 'tiling' technique can be extended to account for large noise levels, where more than one $\pi$ phase bin is spanned [65].

Figure B3 shows some examples of simulated data in the presence of differential phase noise. As for the case of offset noise, estimates of $\phi_d$ exhibit no significant bias over the full range of $0 - \pi$, and for a large range of scale factor ratios. Additionally, only a small number of points are required to converge to a level of uncertainty less than that of the noise defined by $\sigma_{\phi_d}$. The convergence of this uncertainty as a function of the number of measurements is the subject of the next section.

### Scaling with measurement number

To test the scaling of the statistical and systematic error of the Bayesian estimator as a function of the number of measurements, we performed the following study. We randomly generated $M = 50$ samples of 'measurements', each containing 100 points following the model (B.3) with noise added to either the differential
<figure><img src="image_26.png" /><figcaption>**Figure B3.** Simulated Bayesian phase estimates in the presence of differential phase noise. Top row: simulated data for $\kappa = \{0.6, 0.8, 1.0, 1.2, 1.4\}$. Similar to figure B2, common phases are randomly chosen over the range $\phi_c \in [-\pi, \pi]$, and the differential phase is set to $\phi_d^{\text{act}} = 1$ rad. Differential phase noise is applied to each point with a standard deviation of $\sigma_{\phi_d} = 0.1$ rad (SNR~10). Bottom row: the systematic error in the Bayesian phase estimate, $\phi_d^{\text{est}} - \phi_d^{\text{act}}$, as a function of $\phi_d^{\text{act}}$ for each value of $\kappa$. Only 10 points were required to bring the statistical uncertainty to $\lesssim 50$ mrad.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure B4.** The statistical and systematic error in the Bayesian phase estimator as a function of the number of measurements, $N$, used in the Bayesian analysis for (a) differential phase noise and (b) offset noise. In both plots, $\kappa = 0.8$, $\phi_d^{\text{est}} = 1$ rad, and $M = 50$ samples were used. Black points represent the statistical uncertainty $\epsilon_{\phi_d}^{\text{stat}}$, and blue triangles indicate the systematic error $\epsilon_{\phi_d}^{\text{sys}}$. The solid red lines indicate the minimum convergence rates based on $\sigma_{\phi_d}/\sqrt{N}$ for differential phase noise and on $1/\sqrt{N I(\phi_d)}$ for offset noise, where $I(\phi_d)$ is the Fisher information given by equation (B.14). The dashed horizontal lines represent the nominal phase resolution used in the simulations. Values of $\sigma_{\phi_d} = 0.1$ rad (SNR ~ 10) and $\sigma_{B_1} = \sigma_{B_2} = 0.2$ (SNR ~ 5) were used as the noise parameters in (a) and (b), respectively.</figcaption></figure>

phase or the offset. As a function of the measurement number, $N$, within each sample, we computed the Bayesian estimate $\phi_d^{\text{est}}(N)$ and the standard deviation of the associated probability distribution $\delta\phi_d^{\text{est}}(N)$. The statistical error for each measurement is taken as the average of $\delta\phi_d^{\text{est}}(N)$ over all $M$ samples, which we denote as $\epsilon_{\phi_d}^{\text{stat}}(N) = \langle\delta\phi_d^{\text{est}}(N)\rangle_M$. Similarly, the systematic error is defined as $\epsilon_{\phi_d}^{\text{sys}}(N) = \langle|\phi_d^{\text{est}}(N) - \phi_d^{\text{act}}|\rangle_M$. The results are shown in figure B4.

For the specific case of noise that contributes directly to the variable of interest (e.g. differential phase noise) the statistical uncertainty of the Bayesian estimator is given by $\epsilon_{\phi_d}^{\text{stat}} = \sigma_{\phi_d}/\sqrt{N}$. As we show in figure B4(a), the measured statistical error closely follows this dependence. Similarly, on average the systematic error drops to a level much less than $\epsilon_{\phi_d}^{\text{stat}}$ after only a few measurements. This level is primarily determined by the grid resolution used when computing the likelihood distribution for $\phi_d$. During the estimation procedure, we initially set the phase grid resolution to ~$\pi/100$, and we refine this grid size on a measurement-by-measurement basis. As the likelihood distribution narrows, grid points are redistributed toward the maximum likelihood value. We find that this grid optimization procedure can improve the resolution by up to an order of magnitude (depending on the level of noise in the system), while keeping the number of integral evaluations per measurement fixed.
For the more general case of noise present in a parameter that is indirectly related to the quantity of interest through some function, the uncertainty is constrained by the Cramer–Rao lower bound

$$ \sigma_{\phi_d}^{\text{est}} \geq \frac{1}{\sqrt{N I(\phi_d)}}. \qquad (B.13) $$

This relationship can be used to compute the minimum convergence rate of the statistical error in the presence of offset or amplitude noise, for example, where the noise affects $\phi_d$ indirectly through the quantities $\{n_1, n_2\}$. The Cramer–Rao lower bound includes the Fisher information, $I(\phi_d)$, of an individual measurement, which can be computed from the likelihood distribution $L(\{n_1, n_2\}|\phi_d)$ as follows

$$ I(\phi_d) = - \left\langle \frac{\partial^2}{\partial \phi_d^2} \ln[L(\{n_1, n_2\}|\phi_d)] \right\rangle_{\{n_1, n_2\}} = \iint \frac{d n_1 d n_2}{L(\{n_1, n_2\}|\phi_d)} \left(\frac{\partial L}{\partial \phi_d}\right)^2. \qquad (B.14) $$

Here, the brackets $\langle \cdots \rangle_{\{n_1, n_2\}}$ denote an average over the random variables $\{n_1, n_2\}$. The Fisher information is a measure of the amount of information that a random variable (or a set of random variables) carries about an unknown parameter. In this case, the unknown parameter of interest is $\phi_d$ and the set of random variables is the set of measurements $\{n_1, n_2\}$, which are governed by the likelihood distribution $L(\{n_1, n_2\}|\phi_d)$—hence its appearance in equation (B.14). This quantity has no closed-form expression for the case of offset or amplitude noise in our system, and must be evaluated numerically. For the parameters used in figure B4(b), we find $I(\phi_d) \approx 3.3$, which gives a minimum convergence rate of $0.55/\sqrt{N}$. This rate is consistent with the measured statistical uncertainties shown in the figure. We note that the Fisher information empirically scales as $I \sim e^{-\beta\sigma_B}$, where $\beta$ is a large factor that depends on the differential phase and the scale factor ratio used (e.g. $\beta \sim 35$ for $\kappa = 0.8$ and $\phi_d = 1$ rad). Thus, with only a moderate reduction to the level of offset noise in the system, one can dramatically improve the convergence rate of the Bayesian estimate.

## Appendix C. Response of a dual-species interferometer to mirror vibrations

Here, we summarize the essential theoretical tools required to evaluate the response of both single- and dual-species interferometers to vibrational noise of the retro-reflection mirror.

First, we provide a review of the sensitivity function for a single atom interferometer, $g(t)$. This function characterizes how the interferometer transition probability behaves in the presence of fluctuations in the Raman laser phase difference, $\phi_L(t)$. Developed previously for use with atomic clocks [82], the sensitivity function is a useful tool that can be applied, for example, to evaluate the response of the interferometer to laser phase noise [57], or to correct for spurious vibrations in the Raman beam optics [41, 58, 69]. We are primarily interested in the latter.

The sensitivity function is a unitless quantity that is defined as follows

$$ g(t) = \lim_{\delta\varphi \to 0} \frac{\delta\Phi(\delta\varphi, t)}{\delta\varphi} = \lim_{\delta\varphi \to 0} \frac{\delta P(\delta\varphi, t)}{\delta\varphi}, \qquad (C.1) $$

where $\delta\varphi$ is a phase jump occurring at time $t$ during the interferometer that modifies the total interferometer phase, $\Phi$, by an amount $\delta\Phi$, and the transition probability $P(\Phi) = (1 - \cos\Phi)/2$ by a corresponding amount $\delta P$. Thus, the interferometer phase due to an arbitrary phase noise function, $\varphi(t)$, can be computed as

$$ \Phi_{\varphi} = \int_{-\infty}^{\infty} g(t) d\varphi(t) = \int_{-\infty}^{\infty} g(t) \frac{d\varphi(t)}{dt} dt. \qquad (C.2) $$

The quantum mechanical nature of the atom plays a crucial role on the sensitivity function—in particular, the evolution of the internal atomic states during each Raman pulse. Using the procedure outlined in [41, 57], the sensitivity function, $g_j(t)$, of an interferometer with timing parameters labeled with subscript 'j' can be shown to be

$$ g_j(t) = \begin{cases} -\sin(\Omega_j^{\text{eff}}(t - \Delta T_j)) & 0 < t - \Delta T_j \le \tau_j, \\ -1 & \tau_j < t - \Delta T_j \le T_j + \tau_j, \\ -\sin(\Omega_j^{\text{eff}}(t - \Delta T_j - T_j)) & T_j + \tau_j < t - \Delta T_j \le T_j + 3\tau_j, \\ 1 & T_j + 3\tau_j < t - \Delta T_j \le 2T_j + 3\tau_j, \\ -\sin(\Omega_j^{\text{eff}}(t - \Delta T_j - 2T_j)) & 2T_j + 3\tau_j < t - \Delta T_j \le 2T_j + 4\tau_j, \\ 0 & \text{otherwise.} \end{cases} \qquad (C.3) $$

Here, $T_j$ is the interrogation time, $\tau_j$ is a pulse duration, $\Omega_j^{\text{eff}}$ is the effective Rabi frequency associated with the two-photon Raman transitions, and $\Delta T_j$ is a delay with respect to $t=0$ that facilitates a difference in the start
<figure><img src="image_28.png" /><figcaption>Figure C1. Weight functions, $w_j(t)$, described by the response function (C.5). These weights determine the phase shift associated with mirror vibrations in equation (C.4). The pulse durations, $\tau_j$, satisfy $\Omega_j^{\text{eff}} \tau_j = \pi/2$. The differential weight function, i.e. the difference between the red and blue curves, is shown in black.</figcaption></figure>

time between interferometers. It is assumed that $\Omega_j^{\text{eff}} \tau_j = \pi/2$, such that the first and third interferometer pulses have pulse areas of $\pi/2$ with duration $\tau_j$, and the second is a $\pi$-pulse of duration $2\tau_j$.

To evaluate the response of an interferometer to Raman mirror motion, the phase noise function is first expressed as $\phi_j(t) = k_j^{\text{eff}} z(t)$, with $z(t)$ representing the time-dependent position of the mirror along the axis of the beams. Then, the phase shift of interferometer j due to movement of the Raman mirror is

$$ \phi_j^{\text{vib}} = \int_{-\infty}^{\infty} w_j(t) a^{\text{vib}}(t) dt = k_j^{\text{eff}} \int_{-\infty}^{\infty} f_j(t) a^{\text{vib}}(t) dt, \quad (C.4) $$

where $a^{\text{vib}}(t) = \ddot{z}(t)$ is the time-dependent acceleration of the mirror due to vibrations, $w_j(t) = k_j^{\text{eff}} f_j(t)$ is a time-dependent weight function for the mirror accelerations, and $f_j(t)$ is called the response function associated with the jth interferometer. This function is given by the integral of the sensitivity function:

$$ f_j(t) = - \int_0^t g_j(t') dt', \quad \text{and can be evaluated as} $$

$$ f_j(t) = \begin{cases} \frac{1}{\Omega_j^{\text{eff}}} \left(1 - \cos \Omega_j^{\text{eff}} (t - \Delta T_j)\right) & 0 < t - \Delta T_j \le \tau_j, \\ (t - \Delta T_j) + \frac{1}{\Omega_j^{\text{eff}}} - \tau_j & \tau_j < t - \Delta T_j \le T + \tau_j, \\ T_j + \frac{1}{\Omega_j^{\text{eff}}} \left(1 - \cos \Omega_j^{\text{eff}} (t - \Delta T_j - T)\right) & T_j + \tau_j < t - \Delta T_j \le T_j + 3\tau_j, \\ 2T_j + 3\tau_j + \frac{1}{\Omega_j^{\text{eff}}} - (t - \Delta T_j) & T_j + 3\tau_j < t - \Delta T_j \le 2T_j + 3\tau_j, \\ \frac{1}{\Omega_j^{\text{eff}}} \left(1 - \cos \Omega_j^{\text{eff}} (t - \Delta T_j - 2T_j)\right) & 2T_j + 3\tau_j < t - \Delta T_j \le 2T_j + 4\tau_j, \\ 0 & \text{otherwise.} \end{cases} \quad (C.5) $$

At its heart, equation (C.4) is a generalization of the well-known interferometer phase shift due to a constant acceleration, $a$

$$ \phi_j = S_j a = k_j^{\text{eff}} (T_j + 2\tau_j) \left( T_j + \frac{4\tau_j}{\pi} \right) a \approx k_j^{\text{eff}} T_j^2 a. \quad (C.6) $$

In this relation, the quantity $S_j \approx k_j^{\text{eff}} T_j^2$ is equivalent to the integral of the weight function, $w_j(t)$, which determines how strongly the mirror vibration at time $t$ contributes to the interferometer phase shift. This function is triangle-shaped, as shown in figure C1, which indicates that the phase contributions are smallest near $t = \Delta T_j$ and $\Delta T_j + 2T_j + 4\tau_j$, where the wavepacket separation is a minimum. Similarly, the weights are largest near the mid-point, $t = \Delta T_j + T_j + 2\tau_j$, where the separation between the interfering states is a maximum.

For the case of two coupled interferometers, the differential phase shift resulting from mirror vibrations can be expressed as

$$ \phi_d^{\text{vib}} = \phi_1^{\text{vib}} - \phi_2^{\text{vib}} = \int_{-\infty}^{\infty} w_d(t) a^{\text{vib}}(t) dt, \quad (C.7) $$

where the differential weight function, $w_d(t)$, is given by the difference between the single-sensor weight functions
<figure><img src="image_29.png" /><figcaption>Figure C2. (a) Normalized transfer functions, $|H_j(\omega)|/k_j^{\text{eff}} T_j^2$, described by equation (C.10) for coupled K-Rb interferometers. These functions determine the response of the single-species (red curve) and dual-species (black curve) interferometers to acceleration noise at different frequencies, $\omega$. Here, $T_1 = 1$ s and $T_2 = (1 - \epsilon)T_1$, with $\epsilon = 1 - \sqrt{k_1^{\text{eff}}/k_2^{\text{eff}}} \approx 0.0087$, and $\tau_1 = \tau_2 = 10$ µs. (b) Model curves for the power spectral density of ground vibrations described by equation (C.13). The solid red curve corresponds to the power spectral density of a 'quiet' location [66, 67] that we simulated with the parameters $\xi = 0.5$, $\chi = 1$, $\omega_0/2\pi = 10$ Hz, $a_{\text{rms}} = 1.4 \times 10^{-4}$ g. The dashed blue curve corresponds to 'low-noise' conditions achievable with passive vibration isolation ($\xi = 2$, $\chi = 0.8$, $\omega_0/2\pi = 1$ Hz, $a_{\text{rms}} = 1.4 \times 10^{-6}$ g). These two curves were used to compute the rms phase noise in table 1.</figcaption></figure>

$$w_d(t) = w_1(t) - w_2(t) = k_1^{\text{eff}} f_1(t) - k_2^{\text{eff}} f_2(t). \qquad (C.8)$$

This function has an intuitive understanding. For the extreme case when $k_1^{\text{eff}} = k_2^{\text{eff}}$ and the two interferometers are perfectly overlapped (i.e. $\Delta T_1 = \Delta T_2$, $T_1 = T_2$, $\tau_1 = \tau_2$), $w_d(t)$ is zero everywhere. This implies that the differential phase shift due to mirror motion is $\phi_d^{\text{vib}} = 0$—corresponding to perfect common-mode phase noise rejection. In the opposite extreme, when either $k_1^{\text{eff}} \neq k_2^{\text{eff}}$ or the interferometers are not well-overlapped, vibration noise induces a differential phase shift $\phi_d^{\text{vib}}$ between the two sensors given by equation (C.7). This non-zero phase shift is directly responsible for uncorrelated contributions to $\phi_d$ in the case of non-overlapped interferometers, and it explains the loss of common-mode rejection in the case of coupled interferometer with different scale factors. For the case of a constant acceleration, equation (C.7) can also be used to derive the systematic shift $\delta\phi_d^{\text{sys}} = (S_1 - S_2)a$ resulting from interferometers exhibiting $S_1 \neq S_2$.

One can characterize how mirror vibrations with a given frequency spectrum affect each interferometer by computing the mean-squared phase noise

$$(\phi_{\text{rms}}^{\text{vib}})^2 = \frac{1}{2\pi} \int_{0}^{\infty} |H_j(\omega)|^2 S_a(\omega) d\omega. \qquad (C.9)$$

Here, $S_a(\omega)$ is the power spectral density of acceleration noise on the mirror, and $H_j(\omega)$ is the transfer function associated with interferometer j given by the Fourier transform of $w_j(t)$. The transfer function describes how acceleration noise at a given frequency affects the phase over the duration of the interferometer. For frequencies $\omega \ll \Omega_j^{\text{eff}}$ and pulse separations $T_j \gg \tau_j$, this function is well-approximated by

$$H_j(\omega) = -i e^{-i\omega(\Delta T_j + T_j + 2\tau_j)} k_j^{\text{eff}} T_j^2 \text{sinc}^2\left(\frac{\omega T_j}{2}\right). \qquad (C.10)$$

For the dual-species interferometer, one uses the differential transfer function in the same fashion:

$$H_d(\omega) = H_1(\omega) - H_2(\omega). \qquad (C.11)$$

These functions are shown in figure C2 (a) for realistic experimental parameters associated with a K-Rb interferometer. Here, there is a clear difference between the transfer functions associated with single-species and dual-species interferometers. For the individual sensors, the transfer function is well-approximated by the square of a sinc function—which exhibits regular zeroes at the fundamental frequency $1/T_j$ and an envelope that decreases as $(2/\omega T_j)^2$. This dependence implies that the interferometer naturally filters the high-frequency components of the vibration spectrum, with a -3 dB cut-off frequency of $\omega_j^{\text{cut}}/2\pi = \sqrt{2}/\pi T_j \approx 1/2T_j$.

The differential transfer function, on the other hand, has a much more complicated frequency dependence. We will focus on the most interesting case for WEP tests, i.e. when the wave vectors satisfy $k_1^{\text{eff}} = (1 - \epsilon)^2 k_2^{\text{eff}}$ where $\epsilon \ll 1$, and the two interferometers are symmetrically overlapped in time as shown in figure C1. Under these conditions, we find that $H_d$ can be approximated by
$$|H_d(\omega)| \approx k_1^{\text{eff}} (T_1^2 - T_2^2) \text{sinc}\left(\frac{\omega(T_1 + T_2)}{2}\right) \text{sinc}\left(\frac{\omega(T_1 - T_2)}{2}\right) - 2\epsilon k_1^{\text{eff}} T_2^2 \text{sinc}^2\left(\frac{\omega T_2}{2}\right). \qquad (C.12)$$

It follows that there is a competition between the two terms in this expression. For the extreme case when $\epsilon = 0$ (i.e. $k_1^{\text{eff}} = k_2^{\text{eff}}$), the differential transfer function is dominated by the first term, which is identically zero for all frequencies only if $T_1 = T_2$. This represents the ideal case for gravity gradiometry applications. On the other hand, when $\epsilon > 0$ and $T_1 = T_2$, the second term in equation (C.12) dominates. Since the two interferometers are assumed to have different wave vectors, it is not possible to make the transfer function zero at all frequencies. However, it is straightforward to show that $H_d = 0$ in dc provided that $k_1^{\text{eff}} T_1^2 = k_2^{\text{eff}} T_2^2$. This criteria optimizes the rejection of common-mode vibration noise at frequencies below the cut-off for a single-sensor, $\omega_j^{\text{cut}}$, and can be achieved by adjusting the interrogation times such that $T_2 = \sqrt{k_1^{\text{eff}}/k_2^{\text{eff}}} T_1 = (1 - \epsilon) T_1$.

Figure C2(a) shows a comparison between single-sensor and differential transfer functions for $T \sim 1$ s interferometers. When operated differentially, the sensitivity to vibrations at frequencies less than $\omega_j^{\text{cut}}$ is typically more than three orders of magnitude below that of the single interferometer, despite the fact that $k_1^{\text{eff}} \approx k_2^{\text{eff}}$.

Figure C2(b) displays the power spectral density function, $S_a(\omega)$, that was used to compute the rms phase noise in table 1. These curves are based on a regression model for ground accelerations [83] described by

$$S_a(\omega) = N(\xi, \chi) \left( \frac{4\xi\omega_0\omega^3}{(\omega^2 - \omega_0^2)^2 + (2\xi\omega_0\omega)^2} \right)^\chi \frac{2\pi a_{\text{rms}}^2}{\omega}. \qquad (C.13)$$

This model has a single peak at $\omega = \omega_0$ and contains two positive shape parameters, $\xi$ and $\chi$, that determine the sharpness of the peak and the scaling of the wings of the distribution, respectively. The quantity $a_{\text{rms}}^2$ is the mean-squared acceleration of the corresponding time-domain acceleration signal, $a(t)$, and $N(\xi, \chi)$ is a normalization factor that depends on the shape parameters. This factor is chosen such that the integral $\int_{-\infty}^{\infty} S_a(\omega) d\omega = 2\pi a_{\text{rms}}^2$.

## References

[1] Damour T, Piazza F and Veneziano G 2002 Phys. Rev. D 66 046007

[2] Will C M 2006 Living Rev. Relativ. 9 3

[3] Damour T and Nordtvedt K 1993 Phys. Rev. Lett. 70 2217

[4] Damour T and Nordtvedt K 1993 Phys. Rev. D 48 3436

[5] Hohensee M a, Müller H and Wiringa R B 2013 Phys. Rev. Lett. 111 151102

[6] Will C M 1993 Theory and Experiment in Gravitational Physics (Cambridge: Cambridge University Press)

[7] Altschul B et al 2015 Adv. Space Res. 55 501–24

[8] Williams J G, Turyshev S G and Boggs D H 2004 Phys. Rev. Lett. 93 261101

[9] Schlamminger S, Choi K Y, Wagner T A, Gundlach J H and Adelberger E G 2008 Phys. Rev. Lett. 100 041101

[10] Wagner T A, Schlamminger S, Gundlach J H and Adelberger E G 2012 Class. Quantum Grav. 29 184002

[11] Nobili A M et al 2012 Class. Quantum Grav. 29 184011

[12] Overduin J, Everitt F, Worden P and Mester J 2012 Class. Quantum Grav. 29 184012

[13] Bergé J, Touboul P and Rodrigues M 2015 J. Phys.: Conf. Ser. 610 012009

[14] Fray S, Diez C A, Hänsch T W and Weitz M 2004 Phys. Rev. Lett. 93 240404

[15] Bonnin A, Zahzam N, Bidel Y and Bresson A 2013 Phys. Rev. A 88 043615

[16] Tarallo M G, Mazzoni T, Poli N, Sutyrin D V, Zhang X and Tino G M 2014 Phys. Rev. Lett. 113 023005

[17] Schlippert D, Hartwig J, Albers H, Richardson L L, Schubert C, Roura A, Schleich W P, Ertmer W and Rasel E M 2014 Phys. Rev. Lett. 112 203002

[18] Zhou L et al 2015 Phys. Rev. Lett. 115 013004

[19] Bordé C J 1989 Phys. Lett. A 140 10

[20] Berman P R (ed) 1997 Atom Interferometry (San Diego: Academic)

[21] Cronin A, Schmiedmayer J and Pritchard D E 2009 Rev. Mod. Phys. 81 1051

[22] Tino G M and Kasevich M A (ed) 2014 Proc. Int. School of Physics Enrico Fermi (Atom Interferometry) (Amsterdam: IOS Press)

[23] Peters A, Chung K Y and Chu S 1999 Nature 400 849

[24] Peters A, Chung K Y and Chu S 2001 Metrologia 38 25

[25] Sorrentino F, Bertoldi A, Bodart Q, Cacciapuoti L, de Angelis M, Lien Y H, Prevedelli M, Rosi G and Tino G M 2012 Appl. Phys. Lett. 101 114106

[26] Gillot P, Francis O, Landragin A, Pereira dos Santos F and Merlet S 2014 Metrologia 15 L15

[27] Fixler J B, Foster G T, McGuirk J M and Kasevich M A 2007 Science 315 74

[28] Rosi G, Sorrentino F, Cacciapuoti L, Prevedelli M and Tino G M 2014 Nature 510 518

[29] Biedermann G W, Wu X, Deslauriers L, Roy S, Mahadeswaraswamy C and Kasevich M A 2015 Phys. Rev. A 91 033629

[30] Snadden M J, McGuirk J M, Bouyer P, Haritos K G and Kasevich M A 1998 Phys. Rev. Lett. 81 971

[31] McGuirk J M, Foster G T, Fixler J B, Snadden M J and Kasevich M A 2002 Phys. Rev. A 65 033608

[32] Wu X 2009 Gravity gradient survey with a mobile atom interferometer PhD Thesis Stanford University
[33] Sorrentino F, Bodart Q, Cacciapuoti L, Lien Y H, Prevedelli M, Rosi G, Salvi L and Tino G M 2014 *Phys. Rev. A* **89** 023607

[34] Rosi G, Cacciapuoti L, Sorrentino F, Menchetti M, Prevedelli M and Tino G M 2015 *Phys. Rev. Lett.* **114** 013001

[35] Gustavson T L, Bouyer P and Kasevich M A 1997 *Phys. Rev. Lett.* **78** 2046–9

[36] Canuel B, Leduc F, Holleville D, Gauguet A, Fils J, Virdis A, Clairon A, Dimarcq N, Bordé C J, Landragin A and Bouyer P 2006 *Phys. Rev. Lett.* **97** 010402

[37] Barrett B, Geiger R, Dutta I, Meunier M, Canuel B, Gauguet A, Bouyer P and Landragin A 2014 *C. R. Phys.* **15** 875–83

[38] Tackmann G, Berg P, Abend S, Schubert C, Ertmer W and Rasel E M 2014 *C. R. Phys.* **15** 884–97

[39] Dimopoulos S, Graham P, Hogan J M, Kasevich M A and Rajendran S 2008 *Phys. Rev. D* **78** 122002

[40] Canuel B, Amand L, Bertoldi A, Chaibi W, Geiger R, Gillot J, Landragin A, Merzougui M, Riou I, Schmid S P and Bouyer P 2014 *E3S Web Conf.* **4** 01004

[41] Barrett B, Gominet P A, Cantin E, Antoni-Micollier L, Bertoldi A, Battelier B, Bouyer P, Lautier J and Landragin A 2014 Mobile and remote inertial sensing with atom interferometers *Proc. Int. School of Physics Enrico Fermi* (Atom Interferometry vol 188) ed G M Tino and M A Kasevich (Amsterdam: IOS) pp 493–555

[42] Merlet S, Bodart Q, Malossi N, Landragin A, Pereira Dos Santos F, Gitlein O and Timmen L 2010 *Metrologia* **47** 9–11

[43] Dimopoulos S, Graham P W, Hogan J M and Kasevich M A 2007 *Phys. Rev. Lett.* **98** 111102

[44] Dickerson S M, Hogan J M, Sugarbaker A, Johnson D M S and Kasevich M A 2013 *Phys. Rev. Lett.* **111** 083001

[45] Sugarbaker A, Dickerson S M, Hogan J M, Johnson D M S and Kasevich M A 2013 *Phys. Rev. Lett.* **111** 113002

[46] Hartwig J, Abend S, Schubert C, Schlippert D, Ahlers H, Posso-Trujillo K, Gaaloul N, Ertmer W and Rasel E M 2015 *New J. Phys.* **17** 035011

[47] Vogel A *et al* 2006 *Appl. Phys. B* **84** 663

[48] Varoquaux G, Nyman R A, Geiger R, Cheinet P, Landragin A and Bouyer P 2009 *New J. Phys.* **11** 113010

[49] Van Zoest T *et al* 2010 *Science* **328** 1540

[50] Müntinga H *et al* 2013 *Phys. Rev. Lett.* **110** 093602

[51] Sorrentino F *et al* 2011 *J. Phys. Conf. Ser.* **327** 012050

[52] Tino G M *et al* 2013 *Nucl. Phys. B* **243–244** 203–17

[53] Schubert C *et al* 2013 Differential atom interferometry with $^{87}$Rb and $^{85}$Rb for testing the UFF in STE-QUEST (arXiv:1312.5963)

[54] Aguilera D *et al* 2014 *Class. Quantum Grav.* **31** 115010

[55] Kuhn C C N, McDonald G D, Hardman K S, Bennetts S, Everitt P J, Altin P A, Debs J E, Close J D and Robins N P 2014 *New J. Phys.* **16** 073035

[56] Chen X, Zhong J, Song H, Zhu L, Wang J and Zhan M 2014 *Phys. Rev. A* **90** 023609

[57] Cheinet P, Canuel B, Pereira Dos Santos F, Gauguet A, Yver-Leduc F and Landragin A 2008 *IEEE Trans. Instrum. Meas.* **57** 1141–8

[58] Geiger R 2011 Senseur inertiel àdes de matière aéroporté *PhD Thesis* Université Paris-Sud XI

[59] Szpak Z L, Chojnacki W and van den Hengel A 2012 Guaranteed ellipse fitting with the sampson distance *Computer Vision—ECCV 2012* (Lecture Notes in Computer Science vol 7576) ed A Fitzgibbon *et al* (Berlin: Springer) pp 87–100

[60] Foster G T, Fixler J B, McGuirk J M and Kasevich M A 2002 *Opt. Lett.* **27** 951

[61] Lamporesi G, Bertoldi A, Cacciapuoti L, Prevedelli M and Tino G M 2008 *Phys. Rev. Lett.* **100** 050801

[62] Fitzgibbon A, Pilu M and Fisher R B 1999 *IEEE Trans. Pattern Anal. Mach. Intell.* **21** 476

[63] Szpak Z L, Chojnacki W and van den Hengel A 2014 *J. Math. Imaging Vis.* **52** 173–99

[64] Gelman A, Carlin J B, Stern H S and Rubin D B 2003 *Bayesian Data Analysis* 2nd edn (London, Boca Raton, FL: Chapman and Hall, CRC Press)

[65] Stockton J, Wu X and Kasevich M A 2007 *Phys. Rev. A* **76** 033613

[66] le Gouët J, Mehlstäubler T E, Kim J, Merlet S, Clairon A, Landragin A and Pereira Dos Santos F 2008 *Appl. Phys. B* **92** 133

[67] Merlet S, le Gouët J, Bodart Q, Clairon A, Landragin A, Pereira Dos Santos F and Rouchon P 2009 *Metrologia* **46** 87–94

[68] Farah T, Guerlin C, Landragin A, Bouyer P, Gaffet S, Pereira Dos Santos F and Merlet S 2014 *Gyroscopy Navig.* **5** 266–74

[69] Geiger R *et al* 2011 *Nat. Commun.* **2** 474

[70] Lautier J, Volodimer L, Hardin T, Merlet S, Lours M, Pereira Dos Santos M and andragin L F 2014 *Appl. Phys. Lett.* **105** 144102

[71] Pereira Dos Santos F 2015 *Phys. Rev. A* **91** 063615

[72] Ménoret V 2012 Accéléromètre à atomes froids aéroporté pour un test du principe déquivalence *PhD Thesis* Université Paris-Sud-Paris XI

[73] Nyman R A *et al* 2006 *Appl. Phys. B* **84** 673

[74] Ménoret V, Geiger R, Stern G, Zahzam N, Battelier B, Bresson A, Landragin A and Bouyer P 2011 *Opt. Lett.* **36** 4128

[75] Landini M, Roy S, Carcagni L, Trypogeorgos D, Fattori M, Inguscio M and Modugno G 2011 *Phys. Rev. A* **84** 043432

[76] Gokhroo V, Rajalakshmi G, Easwaran R K and Unnikrishnan C S 2011 *J. Phys. B: At. Mol. Opt. Phys.* **44** 115307

[77] Schmidt M, Senger A, Hauth M, Freier C, Schkolnik V and Peters A 2011 *Gyroscopy Navig.* **2** 170

[78] Hu Z K, Sun B L, Duan X C, Zhou M K, Chen L L, Zhan S, Zhang Q Z and Luo J 2013 *Phys. Rev. A* **88** 043610

[79] Nath D, Kollengode Easwaran R, Rajalakshmi G and Unnikrishnan C S 2013 *Phys. Rev. A* **88** 053407

[80] Salomon G, Fouché L, Wang P, Aspect A, Bouyer P and Bourdel T 2013 *Euro. Phys. Lett.* **104** 63002

[81] Kalman R E 1960 *J. Basic Eng.* **82** 35–45

[82] Dick G J 1987 Local oscillator induced instabilities in trapped ion frequency standards 19th Annual Precise Time and Time Interval (U.S. Naval Observatory) pp 133–47 (www.dtic.mil/docs/citations/ADA502386)

[83] Matsuda S 2008 A regression model for estimating power spectral density function of ground acceleration *Proc. 14th World Conf. on Earthquake Engineering* (http://iitk.ac.in/nicee/wcee/article/14_02-0113.PDF)