2. Limitations of quantum accelerometers

A quantum accelerometer operates by transforming a cloud of cold atoms into two spatially separated clouds in free fall such that the change of their vertical displacement in time mimics the two arms of an interferometer. The manipulation of the atom cloud is achieved by using one laser pulse to split the cloud and then a second to recombine the cloud. When the sensor has been subject to a specific force (e.g. acceleration and/or gravity), the two clouds of atoms exhibit different phase characteristics that can be measured. Counting the number of atoms in each cloud yields the relative phase, which in turn yields the specific force of the platform relative to the inertial frame defined by the freely falling atoms. Such a sensor has the potential to produce very precise measurements of acceleration or gravity (Freier et al., Reference Freier, Hauth, Schkolnik, Leykauf, Schilling, Wziontek, Scherneck, Müller and Peters2016; Ménoret et al., Reference Ménoret, Vermeulen, Le Moigne, Bonvalot, Bouyer, Landragin and Desruelle2018). However, the sensor suffers from two key limitations that must be overcome.

The first challenge is that, in general, the quantum sensor has a low sampling rate that is governed by two features: the time required to produce a cloud of cold atoms; and the time that the atoms spend in free fall. A larger size of the atom cloud and a longer time these atoms spend in free fall will result in a better precision of the acceleration measurement. A trade-off therefore exists between the performance of the sensor and the time between measurements. Some laboratories are exploring this trade-off to enable sample rates of up to 330 Hz (Butts et al., Reference Butts, Kinast, Timmons and Stoner2011; McGuinness et al., Reference McGuinness, Rakholia and Biedermann2012).

The second challenge the sensor faces is a low dynamic range. Let $T$ be half the total interrogation time of the sensor, which is equivalent to half the time of flight for the atoms. For a constant acceleration $a$, the phase shift between the two atom clouds is (Peters et al., Reference Peters, Chung and Chu2001)

(2.1)\begin{equation} \triangle \phi = k_{\mathrm{eff}} a T^2, \end{equation}

where the effective wavenumber is $k_{\mathrm {eff}} \approx 4\pi /\lambda$ and $\lambda$ is the wavelength of the laser. This phase shift is measured by the cold atom accelerometer as

(2.2)\begin{equation} S = N\sin(k_{\mathrm{eff}}aT^2 + \phi_0), \end{equation}

where $N$ is the number of atoms in the atomic cloud and $\phi _0$ is the initial phase. Given $S$, an acceleration measurement is thus obtained by inverting Equation (2.2). However, although a single acceleration $a$ uniquely determines $S$, the reverse is not true; a given $S$ can be obtained from distinct accelerations $a$. Figure 1 illustrates the relationship in Equation (2.2) with $N = 1{,}000$ atoms, $\lambda = 780$ nm, $T=1$ ms and a laser pulse width of $\tau = 1\,\mu {\rm s}$.Footnote ¹ We see that once the acceleration exceeds ${\pm }0{\cdot }05$ m s$^{-1}$, an ambiguity occurs, and multiple accelerations map to the same $S/N$ value indicated by the black dashed line.

Figure 1. Normalised output signal ($S/N$) of the quantum accelerometer as a function of input $a$. An output signal (black dashed line) corresponds with multiple acceleration values as indicated by blue circles and dark green boxes, which may be the underlying acceleration measured

In this work, we are interested in unwrapping the phase of the accelerometer output $S$ to identify the underlying acceleration $a$. In Section 3, we introduce an algorithm that performs this unwrapping by fusing a classical sensor with a quantum sensor. The resulting algorithm overcomes both the low sampling rate and low dynamic range challenges exhibited by the quantum sensor.

3. Phase unwrapping by data fusion

In this section, we introduce our approach to unwrapping the phase of the quantum accelerometer via fusion with a classical accelerometer.Footnote ² We assume that the classical accelerometer is perfectly aligned with the quantum accelerometer. Based on Equation (2.2), we can write a complete noise-free quantum accelerometer measurement model as follows (Bonnin et al., Reference Bonnin, Diboune, Zahzam, Bidel, Cadoret and Bresson2018):

(3.1)\begin{align} a & = f(S, s, n) \nonumber\\ & =\left\{\begin{array}{ll} \dfrac{s}{k_{\mathrm{eff}}T^2}\left[\arcsin\left(\dfrac{S}{N}\right) + 2n\pi + \phi_0 \right], & s = 1;\\ \dfrac{s}{k_{\mathrm{eff}}T^2}\left[\arcsin\left(\dfrac{S}{N}\right) - (2n+1)\pi-\phi_0 \right], & s ={-}1, \end{array}\right. \end{align}

where $s = \pm 1$ is the sign function and $n\in \mathbb {Z}$ is an integer. In the presence of shot noise, to a good approximation, this gives a signal model of the form (Close et al., Reference Close, Kealy, Moran, Williams, Haine, Szigeti, White, Robins, Freier, Hardman and Legge2019)

(3.2)\begin{equation} a_{q} = f(S,s,n) + \nu_q, \end{equation}

where $\nu _q \sim \mathcal {N}(0, \sigma _q^2)$, $\mathcal {N}(a,b)$ stands for a Gaussian distribution with mean $a$, variance $b$ and $\sigma _q = 1/ (k_{\mathrm {eff}}T^2\sqrt {N})$ signifies the measurement noise dominated by shot noise which is consistent with the experimental result reported by Cheinet et al. (Reference Cheinet, Canuel, Pereira Dos Santos, Gauguet, Yver-Leduc and Landragin2008). Note that the two solution sets to Equation (3.1) are indicated by the green and blue circles in Figure 1. Accordingly, to obtain the unwrapped acceleration, we need to determine which equation to use and estimate the integer $n$.

We solve this dual estimation problem through fusion with a classical accelerometer. Our fusion algorithm is based on the maximum likelihood estimation and comprises two steps. First, using the classical specific force measurement $a_c$ and the output of the cold atom sensor $S$, a rough estimate of $n$ is obtained by inverting Equation (3.1) as follows:

(3.3a)\begin{align} \hat{n}_1 & = \frac{k_{\mathrm{eff}}T^2}{2\pi}a_c - \frac{1}{2\pi}\left(\arcsin\left(\frac{S}{N}\right)-\phi_0\right), \end{align}

(3.3b)\begin{align} \hat{n}_2 & = \frac{k_{\mathrm{eff}}T^2}{2\pi}a_c + \frac{1}{2\pi}\left(\arcsin\left(\frac{S}{N}\right)-\phi_0\right) - \frac{1}{2}. \end{align}

Note that $\hat {n}_1$ and $\hat {n}_2$ are rounded to the nearest integer. Second, the ambiguity-corrected quantum acceleration $a_q$ is obtained by evaluating Equation (3.1) for a finite set of integers centred around $\hat {n}_1$ and $\hat {n}_2$, and choosing the specific force value closest to $a_c$ in the maximum likelihood sense. The initial phase $\phi _0$ is a known constant determined by the hardware configuration. For simplicity, unless stated otherwise, we shall henceforth assume that $\phi _0 = 0$. It is worth mentioning that a poor signal sensitivity with respect to the underlying acceleration occurs near the peaks of the sinusoid (see Figure 1), which may lead large estimation error in Equation (3.3), and we will discuss this more toward the end of Section 4.

While the phase unwrapping using data from a classical accelerometer provides the fusion output of quantum grade accuracy, the data are only available at the sampling rate of quantum accelerometer, which is much lower than that of a high-end classical accelerometer. One way to solve this problem and obtain quantum grade data at the rate of classical data is to estimate the slow varying bias term by filtering, which is observable from the fusion output data. In this work, a simple error state Kalman filter is used for estimating the bias of classical accelerometer data. The fusion process at time $t_k$, based on a maximum likelihood estimation approach, is specifically given by the following steps.

1. 1. Input: $a_c(k)$ from classical accelerometer; $S_k$ from quantum accelerometer.

2. 2. Maximum likelihood method for parameter estimation:

   (3.4)\begin{equation} (n_k^o, s_k^o) = \arg\max_{n_k \in \mathbb{Z}, s_k \pm 1} p(a_{q}(k)\,|\,S_k,a_c(k)), \end{equation}
   where, following Equation (3.3),
   \[ n^o_k = \left\{\begin{array}{ll} \dfrac{k_{\mathrm{eff}}T^2}{2\pi}a_c(k) -\dfrac{1}{2\pi}\arcsin\left(\dfrac{S_k}{N}\right), & s^o_k = 1; \\ \dfrac{k_{\mathrm{eff}}T^2}{2\pi}a_c(k) + \dfrac{1}{2\pi}\arcsin\left(\dfrac{S_k}{N}\right) - \dfrac{1}{2}, & s_k^o={-}1. \end{array}\right. \]

3. 3. Estimate $a_{f}(k)$ using Equation (3.1), the ambiguity-corrected quantum acceleration output as shown in Figure 2, given $S_k$ and $n_k^o$ and $s_k^o$. Its conditional distribution is

   \[ a_q(k)\,|\,S_k,n^o_k,s^o_k \sim \mathcal{N}(a_{q}(k), \sigma_q^2), \]
   where $\sigma _q = 1/(k_{\mathrm {eff}}T^2\sqrt {N})$.

4. 4. Estimate measurement bias $\hat {b}(k)$ from $a_c(k)$. As shown in Figure 2, a simple error state Kalman filter (Groves, Reference Groves2008) is adopted to estimate the bias of $a_c(k)$ based on the ambiguity-corrected quantum acceleration output $a_q$, which serves as the filtering measurement, at the sampling rate of quantum accelerometer.

5. 5. The final fusion output $a_{\mathrm {out}}(k)$ is given by

   (3.5)\begin{equation} a_{\mathrm{out}}(k) = a_c(k) - \hat{b}(k). \end{equation}

Figure 2. Fusion of classical accelerometer and quantum accelerometer. Errors of the classical accelerometer are corrected using low rate quantum accelerometer measurements after ambiguity correction. The estimated (slow varying) bias term $\hat {b}$ of the classical accelerometer reading is obtained using an error state Kalman filter whose measurement is the signal difference $a_c-a_q$

Figure 2 illustrates the above procedure based on the simulation scenario where the classical accelerometer sampling frequency is 200 Hz and quantum accelerometer sampling frequency is 1 Hz.

We now evaluate the performance of the proposed algorithm via Monte Carlo simulations under the condition that the dynamic range of ground truth acceleration is well beyond that of the quantum accelerometer but within that of the classical accelerometer (though this is impossible in practice). At each run, the data of ground truth acceleration are drawn from a uniform distribution. The following configuration for the cold atom sensor is used: $N = 10{,}000$ atoms, $T = 1$ ms half interrogation time, $\tau =1\,{\rm \mu} {\rm s}$ beam splitter pulse width of the laser and a laser wavelength of $\lambda = 780$ nm.

We assume that the output $a_c$ of the classical accelerometer at $t$ can be expressed as (Titterton and Weston, Reference Titterton and Weston2004):

(3.6)\begin{equation} a_{c}(t) = a(t) + b(t) + w(t), \end{equation}

where $a(t)$ signifies the true specific force, $b(t)$ represents measurement bias which is the sum of a constant bias term and a random bias term, the latter is modelled as a first-order Gauss–Markov process with time constant $\tau _a$, and $w(t)$ represents the random errors associated with the sensor. Note that for this simulation, $b(t)$ may also include the accelerometer scale-factor error which has the effect similar to that of bias error, except that the error term increases as the input acceleration increases. In the simulation, both the constant bias and the standard deviation of random bias are set to $1\times 10^{-3}$ m s$^{-2}$, and the time constant is $\tau _a = 2{,}000$ s, which is equivalent to the bias of a high-end tactical grade sensor (Chow, Reference Chow2011). The random error of the accelerometer is assumed to be a zero-mean Gaussian random variable (Quinchia et al., Reference Quinchia, Falco, Falletti, Dovis and Ferrer2013), i.e.

(3.7)\begin{equation} w(t) \sim \mathcal{N}(0, \sigma_c^2), \end{equation}

where $\sigma _c =1{\cdot }4\times 10^{-3}$ m s$^{-2}$ is the standard deviation.

The error statistics for this evaluation are shown in Figures 3 and 4, based on 10,000 Monte Carlo simulations. These histograms show the error between the ambiguity-corrected quantum acceleration estimate $a_{q}$ and the original classical measurement $a_c$ for two situations: the cold atom sensor with shot noise corruption (blue) and without (orange). In Figure 3, the classical acceleration is drawn from $\mathcal {U}(-10\,{\rm m\,s}^{-2},10\,{\rm m\,s}^{-2})$ whereas in Figure 4, the measurement is drawn from $\mathcal {U}(-1{,}000\,{\rm m\,s}^{-2},1{,}000\,{\rm m\,s}^{-2})$. We observe that the main contribution to the fusion error spread without shot noise (orange) is the nonlinear sensitivity of the mapping between the sensor output signal and the underlying acceleration in the low input dynamical acceleration range. We also observe that the spread of the fusion errors increases as the input acceleration dynamical range increases, because the maximum likelihood estimator error increases since the likelihood from Equation (3.4) becomes flat. This error is discussed further in Section 4.

Figure 3. Fusion error statistics ($a_{q}-a_c$) from 10,000 Monte Carlo runs for $a$ drawn from $\mathcal {U}(-10\,{\rm m\,s}^{-2},10\,{\rm m\,s}^{-2})$ in the presence (blue) and absence (orange) of shot noise

Figure 4. Fusion error statistics ($a_{q}-a_c$) via 10,000 Monte Carlo runs for $a$ drawn from $\mathcal {U}(-1{,}000\,{\rm m\,s}^{-2},1{,}000\,{\rm ms}^{-2})$ in the presence (blue) and absence (orange) of shot noise. This graph shows that for high input acceleration magnitude from the classical accelerometer, the fusion error distribution is dominated by the algorithm nonlinear sensitivity. Compared with Figure 3, the fusion error spread increases with the input acceleration dynamical range

In the next section, we demonstrate the performance of our fusion approach using simulations of a one-dimensional inertial navigation scenario.

4. Navigation simulation results

In this section, we evaluate the performance of the proposed fusion algorithm in the context of inertial navigation via Monte Carlo simulations. Practically, there are many error sources including those due to gravity, rotation of the earth and imperfect inertial sensors which will contribute to the acceleration, velocity and position errors of an inertial navigation system on a moving vehicle (Braasch, Reference Braasch2015). To concentrate on the accelerometer performance comparison, we will consider a one-dimensional inertial navigation case where the vehicle is on a perfectly flat, perfectly straight track and only the accelerometer is needed. Real examples of this are the rocket sleds which are used for testing equipment and various sensors. We idealise such a sled that moves in one dimension by assuming that the accelerometer is mounted in parallel to the track and ignoring sway and vibration, so that velocity is the result of integrating the specific force once and the distance along the track is obtained from an additional integration. Note that this is not a realistic scenario as even a nominally horizontal accelerometer will exhibit some sensitivity to the gravity reaction which we ignored here.

As shown in Figure 2, the system consists of a classical accelerometer and a quantum accelerometer. The initial position and velocity of the vehicle are set to zero. At each run, the underlying body specific force is randomly generated (a single realisation is shown in Figure 5), and its spectrum follows a zero-mean Gaussian distribution $a \sim \mathcal {N}(0, \sigma _a^2)$ with standard deviation $\sigma _a = 1$ m s$^{-2}$. For the classical accelerometer, the standard deviation of measurement noise is $\sigma _q = 1{\cdot }4\times 10^{-3}$ m s$^{-2}$ and both the constant bias and the standard deviation of random bias in $b(t)$ are $1\times 10^{-3}$ m s$^{-2}$. For the quantum accelerometer, we assume that its measurement noise is dominated by shot noise with distribution $v_f \sim \mathcal {N}(0, \sigma _s^2)$, where $\sigma _s = 1/(k_{\mathrm {eff}}T_{pi}^2\sqrt {N})$, where the effective wavenumber is $k_{\mathrm {eff}} = 4\pi /\lambda$, $\lambda = 780$ nm, the half interrogation time $T = 1$ ms, the duration of the laser beam splitter is assumed to be $\tau = 1\times 10^{-6}$ s and the average number of atoms per shot is $N = 10{,}000$. Based on these parameters, we have $\sigma _s = 3{\cdot }1\times 10^{-4}$ m s$^{-2}$ (note that the noise level is larger than the achievable value of $3\times 10^{-8}$ m s$^{-2}$ mentioned by Jekeli, Reference Jekeli2005; Cheiney et al., Reference Cheiney, Barrett, Templier, Jolly, Battelier, Bouyer, Porte and Napolitano2019). The sampling rate of the classical accelerometer is 200 Hz, while the quantum accelerometer is set at 1 Hz.

Figure 5. Comparison of specific force values between classical, quantum, fusion and ground truth in a single run. Apart from the quantum accelerometer output which is indicated in green, the other three are overlapped in this figure. We highlight the RMS error difference between the classical and fusion in Figure 6

In this special inertial navigation scenario, we examine two cases. In the first, the computed inertial navigation is driven by the output of the classical accelerometer alone, and in the second, navigation is driven by the fusion of the classical with the quantum following the proposed fusion procedure, as illustrated in Figure 2. We compare the inertial navigation performance of the two cases in terms of root-mean-squared (RMS) errors.

Figure 5 shows a comparison of the estimated and ground truth acceleration values over 1,000 s. Although they are largely overlapped, the curves inside the blue and green boxes are enlarged to highlight their differences. We see from the figure that the measurement of the quantum accelerometer (green curve) cannot follow the ground truth acceleration caused by the signal with wrapped phases. However, our proposed fusion process extends the dynamic range of the quantum accelerometer, and the output fusion signal ($a_{\mathrm {out}}$) yields a smaller error than that of the classical accelerometer. The RMS error difference between classical ($a_c$) and fusion ($a_{\mathrm {out}}$) specific force measurements with respect to ground truth are illustrated in Figure 6, which is computed from 1,000 Monte Carlo runs driven by a random acceleration profile. It shows that the proposed fusion process shown in Figure 2 is able to eliminate the bias and drift of the classical accelerometer by fusing ambiguity-corrected quantum accelerometer output.

Figure 6. RMS errors of specific forces versus time averaged over 1,000 runs. The equivalent standard deviation of white noise on the quantum accelerometer used in this simulation is $3{\cdot }1\times 10^{-4}$ m s$^{-2}$, and on the classical accelerometer is $2\times 10^{-3}$ m s$^{-2}$

The statistical results on the RMS velocity and position errors are shown in Figures 7 and 8, respectively. These results demonstrate a substantial improvement in velocity and position error performances in the inertial navigation scenario through the proposed fusion process over using the classical accelerometer alone.

Figure 7. RMS velocity errors versus time in the inertial navigation experiment averaged over 1,000 runs

Figure 8. RMS position errors versus time in the inertial navigation experiment averaged over 1,000 runs

Nevertheless, our experiments show that the fused acceleration error exists, and it increases with the scale of the underlying body specific force (see Figures 3 and 4). One major contribution to the fusion errors is from the mapping between the output signal $S$ of the quantum accelerometer and the ambiguity-corrected quantum acceleration $a_q$, which presents different sensitivities due to the nonlinear nature of the sine function. Reading $a_q$ from the linear part of the sine function in $S$ can address the problem. This was shown by Bonnin et al. (Reference Bonnin, Diboune, Zahzam, Bidel, Cadoret and Bresson2018), who used two orthogonal phased quantum accelerometers to remove the nonlinear sensitivity problem of the quantum accelerometer.

As shown in Figure 9 – a normalised plot for the expected output signal versus the value of acceleration – we assume that two ‘exactly orthogonal phased’ quantum accelerometers are available for the fusion process. At this stage, we assume that the output switching between the two quantum sensors is determined by selecting the normalised signal output which satisfies $S < N\sqrt {2}/2$, and that there is no switching error in the simulation. Enhanced error performances are observed from the simulation results similar to those shown in Figures 6 and 8 for the two orthogonal phased quantum accelerometer configurations. In addition, it is observed that the fusion error spread shown in Figure 3 is reduced by half with an orthogonal phased quantum accelerometer configuration.

Figure 9. Illustration of the normalised output signals of two orthogonal phased quantum accelerometers versus acceleration. The bold curve highlights the parts of linear sensitivity from the two sensors across the acceleration range, with $\phi _1 = 0$ (red) and $\phi _2 = \pi /2$ (blue)

It is worth mentioning that the errors of an inertial navigation system depend on many factors, and improving the performance of the accelerometers simply leads to other error sources dominating the inertial navigation performance, including gyro errors (Groves, Reference Groves2008), gravity-modelling errors (Jekeli, Reference Jekeli2012), sensor misalignment, bandwidth limitations and positive feedback of height errors through the gravity model (Titterton and Weston, Reference Titterton and Weston2004).

As ongoing research, we will continue our investigation along these lines for handling phase noise between the two orthogonal-phased quantum accelerometers, and simulating an adaptive phase-locked loop implementation for inertial navigation systems of predictable accelerations.