# Optimized squeezing for accurate differential sensing under large phase noise

Robin Corgier, Marco Malitesta, Leonid A Sidorenkov, Franck Pereira dos Santos, Gabriele Rosi, Guglielmo Tino, Augusto Smerzi, Leonardo Salvi, Luca Pezzè

► To cite this version:

Robin Corgier, Marco Malitesta, Leonid A Sidorenkov, Franck Pereira dos Santos, Gabriele Rosi, et al.. Optimized squeezing for accurate differential sensing under large phase noise. Quantum Science and Technology, 2025, 10 (4), pp.045016. 10.1088/2058-9565/adf2d8. hal-05363835

HAL Id: hal-05363835

https://hal.science/hal-05363835v1

Submitted on 13 Nov 2025

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 CC BY 4.0 - Attribution - International License
# Squeezing-enhanced accurate differential sensing under large phase noise

Robin Corgier¹,†, Marco Malitesta²,³,†, Leonid A. Sidorenkov¹, Franck Pereira Dos Santos ¹, Gabriele Rosi³,⁴, Guglielmo M. Tino²,³,⁴,⁵, Augusto Smerzi²,⁵,⁶, Leonardo Salvi³,⁴,⁵, Luca Pezzè²,⁵,⁶

¹LTE, Observatoire de Paris, Université PSL, Sorbonne Université, Université Lille, LNE, CNRS, 61 avenue de l'Observatoire, 75014 Paris, France.

²Istituto Nazionale di Ottica, Consiglio Nazionale delle Ricerche (INO-CNR), Largo Enrico Fermi 6, 50125 Firenze, Italy.

³Department of Physics and Astronomy, Universita di Firenze, Via Sansone 1 50019 Sesto Fiorentino, Italy.

⁴Istituto Nazionale di Fisica Nucleare (INFN), Via Sansone 1, 50019 Sesto Fiorentino, Italy.

⁵European Laboratory for Nonlinear Spectroscopy (LENS), Via Nello Carrara 1, 50019 Sesto Fiorentino, Italy.

⁶QSTAR, Largo Enrico Fermi 2, 50125 Firenze, Italy.

E-mail: robin.corgier@obspm.fr, leonardo.salvi@unifi.it and
luca.pezze@ino.cnr.it

## Abstract.

Atom interferometers are reaching sensitivities fundamentally constrained by quantum fluctuations. A main challenge is to integrate entanglement into quantum sensing protocols to enhance precision while ensuring robustness against noise and systematics. Here, we theoretically investigate differential phase measurements with two atom interferometers using spin-squeezed states, accounting for common-mode phase noise spanning the full $2\pi$ range. We estimate the differential signal using model-free ellipse fitting, a robust method requiring no device calibration and resilient to additional noise sources. Our results show that spin-squeezing enables sensitivities below the standard quantum limit (SQL). Specifically, we identify optimal squeezed states that minimize the differential phase uncertainty, scaling as $N^{-2/3}$, thus overcoming the SQL, while eliminating the bias inherent in ellipse fitting methods. We benchmark our protocol against the Cramér-Rao bound and compare it with hybrid methods that incorporate auxiliary classical sensors. Our findings provide a pathway to robust and high-precision atom interferometry, in realistic noisy environments and using readily available states and estimation methods.
# 1. Introduction

Precision measurements with atomic devices are often achieved through interferometry techniques, where a physical quantity of interest (e.g. an acceleration, a force, or a frequency) is encoded into a detectable phase shift [1-5]. Atom interferometers are now reaching sensitivities fundamentally bounded by quantum fluctuations, in particular given by the standard quantum limit (SQL) [6-8], which applies when uncorrelated atoms are used [9]. To surpass the SQL, several proof-of-principle experiments have explored the engineering of quantum correlations among atoms [10], e.g. generated through a spin-squeezing process [10-13]. This consists in reducing fluctuations in a collective pseudo-spin observable, leading to smaller measurement uncertainty and thus higher precision in sensing.

The creation of spin-squeezed states has been reported in atomic clocks [14-20]. In particular, an optical clock frequency measurement below the SQL has been realized [21]. Entangled sources have also been proposed for inertial sensors [22-26]. Delocalization of entanglement in a superposition of well-separated atomic momentum states has been demonstrated with ultracold atoms [27, 28] and Bose-Einstein condensates [29], leading to the realization of an atomic gravimeter with sensitivity below the SQL [30]. However, these pioneering results were achieved in controlled environments with low phase noise. In practice, stochastic phase fluctuations can displace squeezed states from their optimal operating point: this causes the antisqueezed quadrature to leak in the measurement and deplete the advantage gained over coherent states [31, 32]. Current efforts focus on mitigating this issue, including optimized feedback loops with multiple squeezed ensembles [33, 34], non-demolition measurements [35] or variational optimization of the probe state and the measurement [36, 37].

This work explores squeezing-enhanced interferometry in a regime characterized by large phase noise, specifically with fluctuations uniformly distributed in the full $2\pi$ interval [20]. We consider a differential configuration, where the simultaneous interrogation of multiple atomic sources allows the rejection of arbitrary common-mode noise [38-40]. Differential schemes are ubiquitous in atom interferometry using uncorrelated atoms [41-48]. They have been exploited for precise measurements of fundamental constants [49-51], to demonstrate that gravity induces Aharonov-Bohm phase shifts [52] and to propose innovative tests of general relativity [53, 54], as well as gravitational wave detection in the milli-Hertz range [55, 56]. Tests of the weak equivalence principle of general relativity have been conducted by measuring the differential free-fall acceleration of distinct atomic species, isotopes, or internal states [57-62]. Differential schemes have also been used for precise measurements of rotations [63], gravity cartography [64] and synchronous clock comparison [65, 66]. Finally, it is worth emphasizing that the high-sensitivity interferometers that reached SQL sensitivities [6-8], as mentioned above, operate in a differential configuration.

Here, we consider a differential interferometer scheme that utilizes combinations of spin-coherent and spin-squeezed states. In the presence of significant large common-mode phase noise, the correlated data from the two interferometers, on average, trace an ellipse. Quantum
noise, which we compute from first principles, introduces fluctuations around this ellipse. To estimate the phase shift, we employ a robust method based on ellipse fitting of correlated data [39], a standard technique in differential interferometry that is inherently resistant to additional noise sources. This approach extracts the phase shift of interest from the parameters of the ellipse by fitting the data obtained from the two coupled interferometers. Besides cold atoms, ellipse fitting has been applied to optical Michelson-Morley interferometry [67, 68] and is implemented across various fields. Unlike other methods, such as maximum likelihood or Bayesian estimation [69], this model-free technique does not require a prior calibration of the experimental output probabilities. Our approach is thus readily experimentally accessible and extends beyond previous theoretical studies, which mainly focused on spin-squeezed state advantages under vanishing [70] or sufficiently small [26] phase noise, or addressing the computation of ultimate sensitivity bounds [71, 72].

Incorporating quantum projection noise from entangled states is a central feature of our analysis. This distinguishes our approach from existing methods [67, 68, 73–76], which primarily consider Gaussian fluctuations around a mean ellipse. Notably, quantum fluctuations contribute to the bias in ellipse fitting methods. We identify a class of spin-squeezed states that simultaneously reduce this bias and enhance the sensitivity in estimating the differential phase shift, outperforming schemes with uncorrelated particles even under significant phase noise. Specifically, we determine an optimal value of the squeezing strength that leads to nearly unbiased differential phase measurements in a broad range of signal phase shifts $0 \le \delta\phi \le \pi/2$ and provides a sensitivity improvement by a factor $N^{1/6}$ relative to the SQL, where $N$ is the total number of particles. To emphasize the robustness of the ellipse fitting method, we compare our results with a simplified analytical ellipse model, developed in this study, as well as with a fringe fitting method relying on noise correlation with auxiliary sensors, and finally with the optimal phase sensitivity given by the Cramér-Rao bound.

## 2. Model and Methods

### 2.1. The differential interferometer scheme

The differential scheme under consideration is shown in Fig. 1. It consists of two Ramsey interferometers operating in parallel and using a common laser to generate beam splitters and mirrors. Our approach applies to Mach-Zehnder interferometers as well. Each interferometer, denoted as A and B, is modeled as a two-mode device, $a_1, a_2$ being the modes of interferometer A and $b_1, b_2$ those of interferometer B. The associated bosonic annihilation and creation operators are respectively $\hat{a}_1, \hat{a}_2, \hat{b}_1, \hat{b}_2$ and $\hat{a}_1^\dagger, \hat{a}_2^\dagger, \hat{b}_1^\dagger, \hat{b}_2^\dagger$. We also introduce collective pseudo-spin operators $\hat{J}_{M,x} = (\hat{m}_1^\dagger\hat{m}_2 + \hat{m}_2^\dagger\hat{m}_1)/2$, $\hat{J}_{M,y} = (\hat{m}_1^\dagger\hat{m}_2 - \hat{m}_2^\dagger\hat{m}_1)/2i$ and $\hat{J}_{M,z} = (\hat{m}_1^\dagger\hat{m}_1 - \hat{m}_2^\dagger\hat{m}_2)/2$, where $\{m, M\} = \{a, A\}$ or $\{b, B\}$. These operators satisfy the usual commutation relations $[\hat{J}_{M,i}, \hat{J}_{M,j}] = i\epsilon_{ijk}\hat{J}_{M,k}$ with $\epsilon_{ijk}$ being the Levi-Civita symbol. For simplicity, we assume that both interferometers use a probe state with the same number of particles, $N$. The two-mode output state of the interferometer M is [77]

$$|\psi_M^{\text{out}}(\phi_M)\rangle = \exp(-i\phi_M \hat{J}_{M,y}) |\psi_M^{\text{in}}\rangle, \quad (1)$$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Differential interferometer scheme. (a) Wigner distributions of the two-mode input states of each interferometer considered in this manuscript: (top) a coherent state given by Eq. (2) and corresponding to $\tau = 0$; (middle) a spin-squeezed state given by Eq. (3) with strength $\tau_{\text{ref}}$; (bottom) a spin-squeezed state with strength $\tau^*$. (b) The differential scheme consists of two interferometers interrogated simultaneously by a common laser. The accumulated phase differences $\phi_A$ and $\phi_B$ are read-out by detecting the particle numbers at the output ports of each interferometer. (c-e) Mean value of the normalized population imbalance, $\bar{z}_M$, as a function of the single interferometer phase shift $\phi_M$, for $\tau = \{0, \tau_{\text{ref}}, \tau^*\}$. The shaded area is the corresponding mean square fluctuation. The different panels correspond to the different input states in (a). (f-h) Probability distribution Eq. (5) for a phase noise spanning the full $\phi_{\text{cn}} \in [0, 2\pi]$ interval, and different probe states: (f) two coherent states, (g-h) two spin-squeezed states with strength $\tau_{\text{ref}}$ (g) and $\tau^*$ (h). Here, $\delta\phi = \pi/16$ and $N = 100$.</figcaption></figure>

where $\phi_M$ is the accumulated phase shift and $|\psi_M^{\text{in}}\rangle$ is the generic input state. The transformation $\exp(-i\phi_M \hat{J}_{M,y}) = \exp(i\frac{\pi}{2}\hat{J}_{M,x})\exp(-i\phi_M \hat{J}_{M,z})\exp(-i\frac{\pi}{2}\hat{J}_{M,x})$ is obtained as a sequence of two beam splitters (collective x-rotation) and a phase precession (z-rotation). In this article, we consider the probe state of each individual interferometer being either a spin-coherent state, described by the binomial particle-number distribution [78]

$$|\psi_M^{\text{Coh}}\rangle = \frac{1}{2^{N/2}} \sum_{n=0}^{N} \binom{N}{n}^{1/2} |N-n\rangle_{m_1} |n\rangle_{m_2}, \qquad (2)$$

or a state generated through one-axis twisting (OAT) dynamics [11],

$$|\psi_M^{\text{Squ}}\rangle = \exp(-i\nu\hat{J}_{M,x})\exp(-i\tau\hat{J}_{M,z}^2)|\psi_M^{\text{Coh}}\rangle. \qquad (3)$$

The state Eq. (3) is spin-squeezed [12] for $0 < \tau \le 1/\sqrt{N}$, where $\tau$ is the parameter controlling the OAT evolution. In Eq. (3), OAT is followed by an appropriate rotation of angle $\nu$ around the x axis to minimize the fluctuations of $\hat{J}_z$. In the inset of Fig. 1(a) we show the Wigner distribution for different input states considered in this manuscript: the top distribution is for
the coherent state of Eq. (2), while the middle and bottom distributions are for Eq. (3) and specific values of $\tau$ (indicated as $\tau_{\text{ref}}$ and $\tau^*$, respectively, see below).

In a differential interferometer configuration,

$$ \phi_A = \phi_{\text{cn}} + \delta\phi/2, \qquad (4a) $$

$$ \phi_B = \phi_{\text{cn}} - \delta\phi/2, \qquad (4b) $$

where $\phi_{\text{cn}}$ is an uncontrolled phase that is common to both interferometers and fluctuates from shot to shot with a probability distribution $\mathcal{P}(\phi_{\text{cn}})$. The differential phase shift $\delta\phi = \phi_A - \phi_B$ is estimated from measurement results of the observable $\hat{z}_M = 2\hat{J}_{M,z}/N$, corresponding to the relative number of particles between the two output interferometer modes, normalized to the total atom number. A collection of $n$ measurements $\{z_{A,j}, z_{B,j}\}_{j=1,...,n}$ can be regarded as a set of points in the $z_A - z_B$ plane. The measurement data are distributed according to

$$ \mathcal{P}(z_A, z_B | \delta\phi) = \int_0^{2\pi} d\phi_{\text{cn}} \mathcal{P}_0(z_A | \phi_A) \mathcal{P}_0(z_B | \phi_B) \mathcal{P}(\phi_{\text{cn}}), \qquad (5) $$

where $\mathcal{P}_0(z_M | \phi_M) = |\langle z_M | \psi_M^{\text{out}}(\phi_M) \rangle|^2$ is the probability to measure the normalized population imbalance $-1 \le z_M \le 1$ in the output state $|\psi_M^{\text{out}}(\phi_M)\rangle$ and $|z_M\rangle$ is the eigenstate of the $\hat{z}_M$ operator with eigenvalue $z_M$. Using Eqs. (4a) and (4b), we find that Eq. (5) only depends on $\delta\phi$.

The ultimate sensitivity bound of an unbiased estimate of $\delta\phi$ is given by the Cramèr-Rao bound

$$ \sigma_{\delta\phi}^{\text{CRB}} = (nF)^{-1/2}, \qquad (6) $$

where $n$ denotes the number of repeated independent measurements and

$$ F = \sum_{z_A, z_B} \frac{1}{\mathcal{P}(z_A, z_B | \delta\phi)} \left( \frac{\partial \mathcal{P}(z_A, z_B | \delta\phi)}{\partial \delta\phi} \right)^2 \qquad (7) $$

is the Fisher information computed from the probability distribution Eq (5), the sums running over $z_M = -1, -1 + 1/N, ..., 1 - 1/N, 1$.

## 2.2. Noiseless environment

In a noiseless environment, $\mathcal{P}(\phi_{\text{cn}})$ can be approximated as a Kroeneker delta function, such that $\mathcal{P}(z_A, z_B | \delta\phi) \approx \mathcal{P}(z_A | \phi_A)\mathcal{P}(z_B | \phi_B)$. The phases $\phi_A$ and $\phi_B$ are estimated independently, giving

$$ \delta\phi_{\text{est}} = \phi_{A,\text{est}} - \phi_{B,\text{est}}, \qquad (8) $$

where $\phi_{M,\text{est}}$ is an estimate of $\phi_M$. The corresponding uncertainty, $\sigma_{\delta\phi}^2 = \sigma_{\phi_A}^2 + \sigma_{\phi_B}^2$, is given by the sum of estimation variances in each interferometer. A practical strategy to estimate $\phi = \phi_M$ (we drop the pendix $M$ where unnecessary) is the method of moments, based on inverting the relation

$$ \bar{z}(\phi) \equiv \langle \psi^{\text{out}}(\phi) | \hat{z} | \psi^{\text{out}}(\phi) \rangle = -C_{\tau} \sin(\phi), \qquad (9) $$
where $C_{\tau} = \cos^{N-1}(\tau)$ is the amplitude of sinusoidal phase oscillations and $\bar{z}$ is estimated by averaging over $n$ measurements. The phase is estimated as

$$ \phi_{\text{est}} = \arcsin(-\bar{z}/C_{\tau}). \qquad (10) $$

The corresponding sensitivity is computed by error propagation,

$$ \sigma_{\phi}(\phi) = \frac{\sigma_{z}(\phi)}{\sqrt{n|\partial\bar{z}/\partial\phi|}}, \qquad (11) $$

where $\sigma_{z}^{2}(\phi) = \langle \psi^{\text{out}}(\phi) | [\hat{z} - \bar{z}(\phi)]^{2} | \psi^{\text{out}}(\phi) \rangle$. We find

$$ \sigma_{z}^{2}(\phi) = \cos^{2}(\phi) \sigma_{z}^{2}|_{\phi=0} + \sin^{2}(\phi) \sigma_{z}^{2}|_{\phi=\frac{\pi}{2}}, \qquad (12) $$

with

$$ \sigma_{z}^{2}|_{\phi=0} = \frac{1}{N} + \frac{N-1}{4N} \left( \mathcal{K}_{1} - \sqrt{\mathcal{K}_{1}^{2} + \mathcal{K}_{2}^{2}} \right), \qquad (13a) $$

$$ \sigma_{z}^{2}|_{\phi=\frac{\pi}{2}} = (1 - C_{\tau}^{2}) - \frac{(N-1)\mathcal{K}_{1}}{2N}, \qquad (13b) $$

$\mathcal{K}_1 = 1 - \cos^{N-2}(2\tau)$ and $\mathcal{K}_2 = 4 \sin \tau \cos^{N-2}(\tau)$ [11]. For two coherent states Eq. (2), we have

$$ \sigma_{\delta\phi}^{\text{SQL}} = \sqrt{2} \sigma_{\phi}^{\text{SQL}} = \sqrt{2} n^{-1/2} N^{-1/2}, \qquad (14) $$

and is indicated as standard quantum limit (SQL) for the differential scheme, where $\sigma_{\phi}^{\text{SQL}}$ is the SQL for the single interferometer. With squeezed states, it is possible to minimize Eq. (11) with respect to the squeezing strength $\tau$. The minimum is achieved for $\tau_{\text{ref}} \approx 3^{1/6}/N^{2/3}$ in the limit $N \gg 1$ [9]. From Eq. (11), we have that

$$ \sigma_{\delta\phi}|_{\phi=0, \tau_{\text{ref}}} \approx 3^{1/3} n^{-1/2} N^{-5/6}, \qquad (15) $$

where the optimal phase sensitivity of a single interferometer is obtained at mid-fringe, $\phi = 0$ (corresponding to $\bar{z} = 0$). Figures 1(c) and (d) show $\bar{z}(\phi)$ as a function of $\phi$, for a coherent state, $\tau = 0$, and a spin-squeezed state of strength $\tau = \tau_{\text{ref}}$ respectively. The mean squared error of the population imbalance $\sigma_z(\phi)$ is added vertically to $\bar{z}(\phi)$ and provides the shaded area. Notice that, according to the error propagation formula, Eq. (11), the phase sensitivity, $\sigma_{\phi}$ is obtained as the horizontal width of the shaded area.

## 2.3. Large noise environment

We refer to large noise environment when $\phi_{\text{cn}}$ has a uniform distribution in $[0, 2\pi]$, namely $\mathcal{P}(\phi_{\text{cn}}) = 1/(2\pi)$. In this case, the standard error propagation analysis outlined above cannot be applied since the large phase noise averages out the mean signal of Eq. (9). By exploiting correlations between the output measurements of two interferometers, it is nevertheless possible to estimate the differential phase from the shape of a Lissajous curve. In fact, when
$\phi_{cn}$ spans the full interval $[0, 2\pi]$, and $\bar{z}_M$ has a sinusoidal dependence on the phase as in Eq. (9), the average moments $\bar{z}_A$ and $\bar{z}_B$ satisfy the ellipse equation

$$ \frac{\bar{z}_A^2}{C_{\tau_A}^2} - 2 \cos(\delta\phi) \frac{\bar{z}_A \bar{z}_B}{C_{\tau_A} C_{\tau_B}} + \frac{\bar{z}_B^2}{C_{\tau_B}^2} - \sin^2(\delta\phi) = 0, \quad (16) $$

where $C_{\tau_M}$ denotes the fringe contrast for a squeezed state of strength $\tau_M$. Equation (16) depends strongly on $\delta\phi$: for $\delta\phi = 0$ or $\pi$, Eq. (16) collapses on a straight line, while for $\delta\phi = \pi/2$ Eq. (16) is a circle. Quantum projection noise determines a spreading of measurement data $\{z_{A,j}, z_{B,j}\}_{j=1,...,n}$ around the average ellipse. In particular, the width of the ellipse is dominated by the largest term between $\sigma_{z_M}^2|_{\phi=0}$ and $\sigma_{z_M}^2|_{\phi=\pi/2}$ in Eq. (12). Panels (f) and (g) of Fig. 1 show color plots of the probability distribution, $\mathcal{P}(z_A, z_B|\delta\phi)$ of Eq. (5), for two coherent states, and two spin-squeezed states with $\tau = \tau_{ref}$, respectively. In particular, when considering $\tau_{ref}$, we find that $\sigma_{z_M}^2|_{\phi=\pi/2}$ is much larger than $\sigma_{z_M}^2|_{\phi=0}$, see Eqs. (13a) and (13b). For $\tau \ge \tau_{ref}$, the Lissajous curve deforms into a “butterfly-like shape”.

The dispersion of the data point around the average ellipse is minimized when $\sigma_{z_M}^2|_{\phi=0} = \sigma_{z_M}^2|_{\phi=\pi/2}$, such that the variance $\sigma_{z_M}^2$ becomes phase-independent. This condition is satisfied for a specific squeezing strength $\tau$ given by $\sigma_{z_M}^2|_{\phi=0} = \sigma_{z_M}^2|_{\phi=\pi/2}$ and leads to [79]

$$ \tau = \tau^* \approx \left(\frac{2}{N^5}\right)^{1/6} \approx \left(\frac{2}{3N}\right)^{1/6} \tau_{\text{ref}}. \quad (17) $$

The corresponding variance is $\sigma_{z_M}^2 \approx 2^{-1/3}N^{-4/3}$. Figure. 1(e) shows $z_M$ (line) with mean square error (shaded region) for the case $\tau = \tau^*$. Using Eq. (11) as in the small-noise case would give the following differential phase sensitivity at mid-fringe:

$$ \sigma_{\delta\phi|_{\phi=0,\tau^*}} = \sqrt{2}\sigma_{\phi|_{\phi=0,\tau^*}} \approx 2^{1/3} N^{-1/2} N^{-2/3}. \quad (18) $$

Interestingly, when the fluctuations of $z_M$ are minimized for all $\phi_M$, a gain over the shot noise limit is possible: the price to pay is a less favorable scaling with atom number compared to the optimal squeezing strength, $N^{-2/3}$ rather than $N^{-5/6}$, Eq. (15). As shown in Fig. 1(h), $\tau^*$ provides a narrower distribution of data compared to i) the case of two coherent states ($\tau = 0$) and ii) the case of two optimized spin-squeezed states in a low noise environment ($\tau = \tau_{ref}$). While the probability density concentrates around the major axis of the ellipse for $\tau = 0$ and around the major and minor axes for $\tau = \tau_{ref}$, a trade-off with nearly uniform distribution is obtained for $\tau = \tau^*$.

### 2.3.1. Ellipse fitting and statistical analysis
A fit of the measurement data aims at recovering the “true” ellipse in Eq. (16) and provides an estimate $\delta\phi_{est}$ of the differential phase. The figures of merit we consider are both the variance

$$ \sigma_{\delta\phi_{est}}^2 = \langle \delta\phi_{est}^2 \rangle - \langle \delta\phi_{est} \rangle^2, \quad (19) $$

which quantifies the precision of the estimate, and the deviation of the mean estimate from the true value of the differential phase,

$$ \mathcal{B}(\delta\phi_{est}) = \langle \delta\phi_{est} \rangle - \delta\phi, \quad (20) $$
indicated as bias and quantifying the accuracy. In the following, we outline two different approaches depending on the a priori knowledge about the system.

If the squeezing strength is unknown, we cannot make use of Eq. (16). We thus fit the data {$z_{A,j}, z_{B,j}$}$_{j=1,...,n}$ to the conic function

$$a z_A^2 + b z_A z_B + c z_B^2 + d z_A + e z_B + f = 0, \quad (21)$$

to extract the estimate $\delta\phi_{est}$ by means of the relation [39],

$$\delta\phi_{est} = \arccos\left(\frac{-b}{2\sqrt{ac}}\right). \quad (22)$$

The conic parameters $v = \{a, b, c, d, e, f\}^\top$ are obtained from a least-square approach that minimizes the sum of squared distances,

$$\sum_{j=1}^{n} d^2(z_{A,j}, z_{B,j}; c_v), \quad (23)$$

between data points and the general conic $c_v$, Eq. (21), with parameter vector $v$. The most popular choices for $d(z_{A,j}, z_{B,j}; c_v)$ are the algebraic and the geometric distances, respectively given by, $d_{AL}(z_{A,j}, z_{B,j}; C) = k_j^T v$ with $k_j = (z_{A,j}^2, z_{A,j}z_{B,j}, z_{B,j}^2, z_{A,j}, z_{B,j}, 1)^T$ and $d_G(z_{A,j}, z_{B,j}; C) = \min_{(z_A, z_B) \in C} \sqrt{(z_A - z_{A,j})^2 + (z_B - z_{B,j})^2}$. These distances, respectively, lead to the algebraic and the geometric fits [80]. The geometric distance is a non-linear function and the geometric fit has to be optimized over a large parameter space, that contains the coordinates of the $n$ points $(\tilde{z}_{A,j}, \tilde{z}_{B,j})_{j=1,...,n}$ and the conic parameter vector $v$, with an iterative procedure. The algebraic distance is a linear function but the algebraic fit procedure requires to impose constraints on the conic parameters. In this work, we consider the linear constraints $a + c = 1$, known as the trace constraint, referred to as the “trace method” in the following [81], and the quadratic constraint $b^2 - 4ac = -1$, referred to as ellipse-specific constraint [82–84]. Other fitting constraints can be found in the literature [85–87]. The estimation depends on the expression of $d(z_{A,j}, z_{B,j}; c_v)$ and the choice of the constraint: indeed least-square fit on curved lines is generally biased [73, 88]. In Fig. 2, we compare the performance of a geometric and an algebraic fit on data points sampled according to Eq. (5). In panel (a), with $\tau = 0$, we see that the geometric fit (green dotted line) gets closer to the true ellipse (solid line) than the algebraic one (blue dashed line), at the cost of a non-linear fitting procedure of higher computational complexity. In contrast, for squeezing strength $\tau = \tau^*$, as shown in panel (b), both fits are close to Eq. (16), leading to an accurate estimation of $\delta\phi$, as discussed in the following Section.

Instead, if the squeezing strength is known, then we can fit the data to Eq. (16), where $\delta\phi$ is the only unknown parameter. We indicate this method as one-parameter fit. This provides the estimate:

$$\delta\phi_{est} = \arccos(h), \quad (24)$$

where $h$ is the real solution to the cubic equation

$$G_0 + G_1 h + G_2 h^2 + G_3 h^3 = 0. \quad (25)$$
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Illustration of different ellipse fitting methods for (a) two coherent states and (b) two squeezed states of strength $\tau^*$. The black line is the “true” ellipse given by Eq. (16). The empty gray circles are $n = 500$ randomly data points $\{z_{A,j}, z_{B,j}\}_{j=1,...,n}$ sampled according to Eq. (5). In each panel, the dashed blue and dotted green line show, respectively, the linear algebraic fit (with trace method) and the non-linear geometric fit. Here, $\delta\phi = \pi/4$ and $N = 100$. The insets are zooms.</figcaption></figure>

The cubic coefficients $G_0, G_1, G_2$ and $G_3$ are functions of $\{z_{A,j}, z_{B,j}\}_{j=1,...,n}$ and require the knowledge of the contrast $C_{\tau_M}$ in Eq. (16), see Sec. Appendix A. In the limit of large $n$, we have

$$ \langle \delta\phi_{est} \rangle \approx f(\langle G_0 \rangle, \langle G_1 \rangle, \langle G_2 \rangle, \langle G_3 \rangle), \qquad (26) $$

$$ \sigma_{\delta\phi_{est}}^2 \approx \sum_{j,l=0}^{3} \left( \frac{\partial f}{\partial G_j} \right) \mathrm{Cov}(G_j, G_l) \left( \frac{\partial f}{\partial G_l} \right). \qquad (27) $$

In the above equations $f(\{G_i\}_{i\in[0,3]}) = \arccos[p_{cube}\{G_i\}_{i\in[0,3]}]$ where $p_{cube}(\{G_i\}_{i\in[0,3]})$ corresponds to the Cardano formula expressing the real solution of a cubic equation as a function of its four coefficients $G_i$. $\mathrm{Cov}(G_j, G_l) = \langle G_j G_l \rangle - \langle G_j \rangle \langle G_l \rangle$ and the partial derivatives on the right-hand side of Eq. (27) are evaluated for the arguments $\langle G_i \rangle$ as for the right-hand side of Eq. (26).

### 3. Results

In the following, we show how spin-squeezing enables quantum-enhanced and accurate differential phase measurements in presence of large common phase noise. We first present our results focusing on the trace method (namely, algebraic ellipse-fitting method with linear constraints), leading to large bias with coherent states, but numerically very efficient. We then compare with other methods which are potentially less biased but numerically more demanding, as it is the case of the non-linear geometric fitting method. We focus on different aspects of quantum-enhanced metrology and analyze the bias, the effective single-point phase sensitivity defined as $\sigma_{\delta\phi_{est}}^{\mathrm{eff}} = \sqrt{n} \sigma_{\delta\phi_{est}}$, and the quantum gain with respect to the SQL, defined as $\mathcal{G} = \sigma_{\delta\phi}^{\mathrm{SQL}}/\sigma_{\delta\phi}$. We emphasize the scaling of these quantities with atom number.
3.1. Bias and phase sensitivity

<figure><img src="image_4.png" /><figcaption>**Figure 3.** The bias $\mathcal{B}(\delta\phi_{\text{est}})$ (top panels) and the effective single-point differential phase sensitivity, $\sigma_{\delta\phi_{\text{est}}}^{\text{eff}} = \sqrt{n} \sigma_{\delta\phi_{\text{est}}}$ (bottom panels) are evaluated as a function of the true differential phase $\delta\phi$ (a, c) and as a function of the reduced squeezing strength $\tilde{\tau} = \tau/\tau_{\text{ref}}$ (b, d). Different symbols represent calculations for different values of $\tilde{\tau}$ in (a, c) and $\delta\phi$ in (b, d), as indicated by proper legends. In all panels, filled (open) symbols are for the case of two squeezed states (one squeezed and one coherent state) as input, and the dotted black lines in (a,c) hold for two coherent input states, i.e. $\tilde{\tau} = 0$. The filled symbols are joint with the line for clarity purpose. In (c, d), the purple stars show the Cramèr-Rao bound for the unbiased estimators, Eq. (6). In (b, d), the vertical dashed line denotes the optimal squeezing strength $\tau^*$. For all data, we employed the algebraic fit method with trace constraint on $n_{\text{ell}} = 1000$ ellipses of $n = 1000$ points with $N = 500$ atoms per interferometer. Statistical errors are within the size of the symbols. Note that $n \gtrsim 10^2$ ensures convergence (see Sec.C of the Appendix).</figcaption></figure>

Previous works [39–48] have shown that the bias $\mathcal{B}(\delta\phi_{\text{est}})$ in an ellipse fitting procedure in the presence of noise is reduced for $\delta\phi \approx \pi/2$, while being large around $\delta\phi = 0$. This result is observed for various combinations of probe states, in Fig. 3(a). Using one spin-squeezed state, shown by empty symbols, or two spin-squeezed states, shown by filled symbols, leads to bias reduction for any differential phase when compared to the case of two coherent-states (black dotted line). In Fig. 3(b) we further analyze the bias for two squeezed states of strength $\tau$. For all values of the differential phase we analyzed [the cases $\delta\phi = \pi/16$ and $\delta\phi = \pi/8$ are reported in Fig. 3(b)], we always found a squeezing strength at which the bias curve crosses the zero-bias line, implying a complete cancellation of the bias. Such bias-free configuration depends on $N$ and is realized in the neighborhood of $\tau^*$. A slight shift of the optimal $\tau$ with respect to $\tau^*$ is a finite-size effect due to the relatively small value of $N$ considered: we have
checked that the discrepancy reduces with increasing $N$. Below, in Sec. 3.3, we discuss how the bias at $\tau = \tau^*$ scales with $N$. We anticipate that the bias is much smaller than that found for the coherent state and has a faster decrease with $N$. The configuration $\tau = \tau^*$ is shown by the filled green diamonds in Fig. 3a and highlighted by the vertical dashed green line in Fig. 3b. An unbiased estimation can be also obtained when only a single spin-squeezed is used but at the expense of larger squeezing strength, $\tau > \tau^*$, as shown Fig. 3b.

Figure 3(c,d) shows the effective single-point differential phase sensitivity, $\sigma_{\delta\phi_{est}}^{\text{eff}}$, following the same color code and symbols as in Fig. 3(a,b). First, the effective single-point differential phase sensitivity improvement with two spin-squeezed states (filled symbols) is significantly larger than that obtained with a single squeezed state (empty symbols) where the gain with respect to the SQL is at best $\sqrt{2}$ [26]. Furthermore, the best ellipse fitting result is obtained for two spin-squeezed states of strength $\tau = \tau^*$ where the differential phase sensitivity is essentially independent of the differential phase, as shown by the filled green diamonds in Fig. 3c. Interestingly, in the case $\tau = \tau^*$, high-sensitivity ellipse fitting is possible in the full range $0 \le \delta\phi \le \pi/2$. The sensitivity of the ellipse fitting for the optimal differential spin-squeezed configuration ($\tau_A = \tau_B = \tau^*$) is found only slightly above the Cramèr-Rao bound for unbiased estimators, Eq. (6), shown as the filled purple stars. Saturating the Cramèr-Rao bound would require a more involved parameter estimation analysis, for instance, using a maximum likelihood or a Bayesian method. However, these approaches require the knowledge of the conditional probability distribution Eq. (5) for all possible measurement events $z_A$ and $z_B$. This may be accessed by preliminary calibration of the differential interferometer, which seems rather impractical as it requires collecting a large number of data.

## 3.2. Comparison of different ellipse-fitting methods

To strengthen the above results, we compare four different ellipse fitting methods in Fig. 4. The algebraic fit with trace constraint used in Fig. 3 is shown as the filled blue circles, the algebraic fit with ellipse-specific constraint is shown as the empty red squares, the geometric fit is shown as the filled green diamonds and the one-parameter analytical fit, is shown by the solid gray line. For simplicity, we restrict to the case of two spin-squeezed states with identical squeezing strength.

Independently of the fitting procedure, we recover that $\tau \approx \tau^*$ minimizes the bias, as shown in Fig. 4(a-b) for two different phase differences, and is the best choice of $\tau$ leading to the optimum phase sensitivity (apart from finite size effects, as discussed above), as shown Fig. 4(c-d). In particular, only the “trace” and “geometric” ellipse fitting allow to cancel the bias for any differential phase. In addition, a direct comparison with the Fisher information, leads to the conclusion that geometric ellipse fitting is optimal in the regime of small differential phase when $\tau \approx \tau^*$, as shown in Fig. 4(c).
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Bias, $\mathcal{B}(\delta\phi_{est})$, and effective single-point differential phase sensitivity, $\sigma_{\delta\phi_{est}}^{\text{eff}} = \sqrt{n}\sigma_{\delta\phi_{est}}$, obtained for different ellipse fitting methods: algebraic with trace constraint used in Fig. 3 in filled blue circles, algebraic one-parameter (solid gray line), algebraic with ellipse-specific constraint in red squares and geometric in green diamonds. All methods are discussed in Sec. 2.3. The vertical dashed line highlights the value of $\tau^*$, see Eq. (17). The calculations are performed using $n_{ell} = 2000$ ellipses each containing $n = 1000$ points. Bias and sensitivity for the one-parameter fit are obtained by evaluation of Eqs. (26) and (27). Here we used $N = 500$ atoms per interferometer.</figcaption></figure>

### 3.3. Bias and phase sensitivity scaling with atom number

Figure 5 shows the scaling of the bias absolute value, $|\mathcal{B}(\delta\phi_{est})|$, and of the effective differential phase sensitivity, $\sigma_{\delta\phi_{est}}^{\text{eff}}$, as functions of the atom number $N$ in each interferometer. Here, ellipse fitting is performed using the algebraic fit with trace constraint (empty symbols) and the one-parameter fit (solid lines). In each panel, empty symbols or solid lines refer to the combinations of: (i) two coherent states (black circles), (ii) one coherent and one squeezed state (blue squares) and (iii) two squeezed states (green diamonds and solid gray lines). The squeezing strength is fixed to $\tau^*$. To account for the statistical errors, the one-sigma error bars on the bias, $\sigma_{\delta\phi_{est}} / \sqrt{n}_{ell}$, are shown in Fig. 5(a-b) for each configuration. Here, $n_{ell}$ denotes the total number of sampled ellipses while the number of points per ellipse, $n$, is already accounted for in the differential phase sensitivity, see Eq. (18).

In order to derive the scaling laws for the bias and the sensitivity of the fit with trace constraint, we have performed power-law fitting on the relative data. For each series of symbols, the dashed line corresponds to a fit of the form,

$$ \log_{10}(y) = \alpha - \beta \log_{10}(N), \qquad (28) $$

where the values of $\alpha$ and $\beta$ are reported in table. 1. Interestingly, the scaling of the bias with
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Scaling of the bias, $\mathcal{B}(\delta\phi_{est})$ (a-b), and effective single-point differential phase sensitivity, $\sigma_{\delta\phi_{est}}^{\text{eff}} = \sqrt{n} \sigma_{\delta\phi_{est}}$ (c-d). We consider coherent-coherent (black circles), coherent-squeezed (blue squares) and squeezed-squeezed (green diamonds and solid gray lines) input state configurations. We indicate as $N$ the number of atoms in each interferometer. Vertical rows correspond to different values of $\delta\phi$. In panels (a) and (b) the error bars correspond to the statistical error given by $\sigma_{\delta\phi_{est}} / \sqrt{n}_{\text{ell}}$ where $n_{\text{ell}} = 1000$ is the total number of ellipses containing each $n = 1000$ points. The purple stars are Eq. (6) for two spin-squeezed states. In all panels, the squeezing strength is fixed to $\tau^*$ (see Eq. 17). The vertical dotted lines highlight the configuration with $N = 500$ atoms shown Figs. 3 and 4. The dashed lines corresponds to a linear fit restricted to the range $300 < N < 1000$. The fit coefficients are given in Table. 1.</figcaption></figure>

atom number is more advantageous for two squeezed states where $\mathcal{B}(\delta\phi_{est}) \propto 1/N^{1.3}$, than for only one or zero squeezed states where $\mathcal{B}(\delta\phi_{est}) \propto 1/N$. In this analysis, the numerical fluctuations caused by the finite number $n_{\text{ell}}$ of sampled ellipses, mostly affecting the points that correspond to the smallest bias values, can hinder a careful evaluation of the scaling exponents. Thus, it is worth noting that the same scaling laws are recovered with the one-parameter analytical fit. Using a first-order Taylor expansion of Eq. (26), the bias can be expressed through the compact formula (see Sec. Appendix A of the Appendix for details),

$$ \mathcal{B}(\delta\phi_{est}) \approx -4 \cot(\delta\phi) \frac{H_0 + H_2 h^2}{1 + 2h^2}, \qquad (29) $$

where $h = \cos(\delta\phi)$ and $H_0, H_2$ are functions of $\tau$ and $N$. For $\tau = 0$, we have $H_0 = -7/(4N)$ and $H_2 = 1/N$, implying $\mathcal{B}(\delta\phi_{est}) \sim 1/N$ as found above. On the other hand, for $\tau = \tau^* \sim 1/N^{5/6}$, we have $H_0 \sim \sigma_z^2 \sim 1/N^{4/3}$ (see also Sec. 2.3) and $H_2 = 0$, thus recovering $\mathcal{B}(\delta\phi_{est}) \sim 1/N^{4/3}$.
<table><caption>**Table 1.** Coefficient of Eq. (28) for the data shown in Fig. 5. Notice that scaling coefficients for $δφ = π/16$ are identical to those obtained at $δφ = π/4$.</caption><thead><tr><th>y</th><th>δφ</th><th>{τ<sub>A</sub>, τ<sub>B</sub>}</th><th>α</th><th>β</th></tr></thead><tbody><tr><td rowspan="3">|δφ<sub>est</sub> − δφ|</td><td rowspan="6">π/16</td><td>{0, 0}</td><td>3.67</td><td rowspan="3">≈ 1</td></tr><tr><td>{0, τ<sup>*</sup>}</td><td>3.39</td></tr><tr><td>{τ<sup>*</sup>, τ<sup>*</sup>}</td><td>3.84</td></tr><tr><td rowspan="2">σ<sub>δφ<sub>est</sub></sub></td><td>{0, 0}</td><td>3.37</td><td rowspan="2">≈ 1/2</td></tr><tr><td>{0, τ<sup>*</sup>}</td><td>3.31</td></tr><tr><td rowspan="2">σ<sub>F</sub></td><td>{τ<sup>*</sup>, τ<sup>*</sup>}</td><td>3.38</td><td rowspan="2">≈ 2/3</td></tr><tr><td>{τ<sup>*</sup>, τ<sup>*</sup>}</td><td>3.23</td></tr></tbody></table>

In panel (c) and (d), we recover the scaling of the differential phase sensitivity with atom number discussed above in Sec. 2 in the specific case $φ_M = 0$, i.e. at mid-fringe (see Eqs. 14 and 18). In the case where one or zero squeezed states are used, as respectively shown by the blue squares and black circles, $σ_{δφ_{est}} ∝ N^{-1/2}$, giving a gain $G = σ_{δφ}^{SQL}/σ_{δφ}$ that does not scale with $N$. On the contrary, for two squeezed states we find $σ_{δφ_{est}} ∝ N^{-2/3}$, and thus $G ∝ N^{1/6}$, as shown by the green diamonds. For all configurations, the corresponding phase sensitivity inferred from the Fisher information has the same scaling, as shown by the purple stars for the specific case of two squeezed states of strength $τ^*$. For the latter case, the difference between ellipse and Fisher information is only 1.5 dB. We notice that the gain factor $N^{1/6}$ has been also observed previously when studying the performance of a single squeezed state in an atomic clock [31, 33].

### 3.4. Ellipse vs Hybridization with classical sensors

So far, we have studied a differential atom interferometer scheme shown in Fig. 1 subject to an unknown random common phase noise $φ_{cn}$. Any information about the actual value of $φ_{cn}$ at each measurement shot was considered unnecessary, as this noise serves to distribute the data points around the ellipse and thus facilitates the conic fitting to estimate the differential phase of interest. Now we extend the differential scheme of Fig. 1 by associating it with an additional classical sensor that measures the common phase noise $φ_{cn}$. For instance, the classical sensor can be a seismometer recording the common-mode vibration of the experimental platform that shifts the laser phase with respect to free-falling atoms. This hybrid quantum-classical configuration was proposed in Ref. [40] and demonstrated experimentally in Ref. [43]. This method allows to correlate the contribution of the phase noise $φ_{cn}$ at the j-th iteration with both interferometer outputs $z_{M,j}$. As a consequence, the output data of each single interferometer becomes distributed along sinusoids, similar to the one shown in Fig. 1(c-e), with a phase noise limited by the accuracy of the correlation with the classical sensor. The inertial phase of interest is then directly extracted by the fringe fit [43], simultaneously providing $φ_{A}^{est}$ and $φ_{B}^{est}$, and thus $δφ^{est}$ according to Eq. (8).

This hybrid method has demonstrated a virtually unbiased extraction of differential phase $δφ$, i.e. $B(δφ_{est}) = 0$, with sensitivities close to the SQL when $δφ ≈ 0$ with coherent spin
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Ellipse fitting vs fringe fitting. In all panels, the configuration with hybridization is shown by the empty orange diamonds and obtained by fringe fitting for $\delta\phi = 0$. The configuration without hybridization is shown by the filled green diamonds and obtained by ellipse fitting for the specific case $\delta\phi = 1$ rad. (a-b) Evolution of the effective single shot differential phase sensitivity, $\sigma_{\delta\phi est}^{\text{eff}} = \sqrt{n} \sigma_{\delta\phi est}$, for varying squeezing strengths. The vertical dashed line denotes the value of $\tau^*$. The Cramèr-Rao bound, given by Eq. 6, is shown as the filled purple stars in (a) for $N = 10^3$ atoms. (c) Evolution of the quantum gain, $\mathcal{G} = \sigma_{\delta\phi}^{\text{SQL}}/\sigma_{\delta\phi est}$, as a function of the atom number, $N$, for $\tau = \tau^*$. The lines correspond to a fit given in the caption. The total number of ellipses is $n_{\text{ell}} = 1000$ with each $n = 1000$ points.</figcaption></figure>

states [43]. At $\delta\phi = 0$, the noise correlation error (stemming, for example, from the finite bandwidth of the auxiliary classical sensor) responsible for the residual differential phase noise $\sigma_{\delta\phi}^{\text{corr}} \lesssim 0.5$ rad barely affects the sensitivity of the fringe fitting method. Away from this optimal point, however, the sensitivity rapidly degrades with increasing $\sigma_{\delta\phi}^{\text{corr}}$. Given the in-applicability of the ellipse fitting approach in the vicinity of $\delta\phi = 0$, we emphasize that both methods are complementary since they optimize the sensitivity in different ranges.

Figure 6 compares the performance of the ellipse-fitting method (without hybridization, filled green diamonds) with that achievable using a fringe-fitting method in the hybrid approach (empty orange diamonds). The data employed in the fringe fitting method are sampled from the two input squeezed states of $\tau = \tau^*$ (optimal for ellipse fitting method) given by Eq. (3), with differential phase $\delta\phi = 0$. Here we also restrict the study to the $\sigma_{\delta\phi}^{\text{corr}} = 0$ case, i. e. assuming a perfect knowledge of the $\phi_{cn}$ values. Note that the differential phase sensitivity with ellipse fitting is independent of the differential phase at $\tau^*$, as shown in Fig. 3c and discussed in sec. 3.1, and we choose $\delta\phi = 1$ rad by default. The two methods manifest the same scaling with atom number, while the hybridization method provides about 20% better sensitivity at its optimal point: $\delta\phi = 0$. For a variable squeezing strength $\tau$ in Fig. 6 (a), the sensitivity values provided by the two methods remain very close up to the optimal squeezing point $\tau^*$, where the fringe fitting sensitivity reaches the Cramèr-Rao bound (purple star). A somewhat less smooth behavior of the ellipse as compared to fringe fitting results (i.e. larger point-to-point variation) is likely attributed to stronger impact of data dispersion due to small atom number on the ellipse fitting. The results for a large atom number ($N \gtrsim 10^4$) in Fig. 6b are obtained by approximating the distributions $P_0(z_M|\delta\phi_M)$ in Eq. (5) by Gaussian functions centered in $\bar{z}_M(\phi_M)$, Eq. (9), and with variance given by Eq. (12). Overall, the results confirm
the $N^{1/6}$ scaling of the metrological gain, as shown in Fig. 6 (c).

## 4. Discussion and Conclusion

To summarize, we have explored the possible advantages of using squeezed states for estimating differential phase shifts through standard ellipse fitting techniques. Ellipse fitting offers several key benefits: it requires no calibration of the experimental setup and no prior knowledge of the phase noise model. It is inherently robust against large common-phase noise and operates under the minimal assumptions that the mean signal depends sinusoidally on the phase shift and the interferometric setup is affected by large common-mode noise. The nonlinear nature of ellipse fitting allows for various techniques, with differing complexities, which we have compared in this study. Our analysis shows the possibility to reach sensitivities below the SQL across a broad signal range $0 \le \delta\phi \le \pi/2$. Upon optimizing the squeezing strength, we reach a gain in instrument sensitivity over the SQL with an optimal scaling $N^{1/6}$, and also show that spin-squeezing can remarkably suppress the bias in differential phase estimation with respect to using spin-coherent states. This result is also obtained when considering a hybrid approach where the differential interferometer is assisted by a classical sensor measuring the phase noise value $\phi_{cn}$. Our protocol is practical and of immediate applicability in experiments: it uses quantum states that have been demonstrated in several experiments [10]. Notably, squeezing dynamic in presence of large common phase noise has been realized [20] with Rydberg interactions in a programmable optical atomic clock and used to perform a synchronous frequency comparison between independent squeezed states at times beyond the atom-laser coherence. In this experiment [20], the differential phase shift was extracted by using a maximum-likelihood estimation method given an appropriate noise model.

The possibility to reach the Heisenberg limit of sensitivity with ellipse fitting may require engineering specific quantum probe states [36, 89] that, differently from the squeezed states considered in this manuscript, may not be straightforward to generate in atom interferometry experiments. As atom interferometers are now reaching control levels dominated by quantum noise [6–8], our work represents a significant step toward ultra-precise differential estimation, paving the way for advanced applications in realistic noisy environments.

**Acknowledgments.** This research has been carried out in the frame of the QuantERA project SQUEIS (Squeezing enhanced inertial sensing), funded by the European Union's Horizon Europe Programme, the Agence Nationale de la Recherche (ANR-22-QUA2-0006). We also acknowledge financial support from the European Union's Horizon 2020 research and innovation programme—Qombs Project, FET Flagship on Quantum Technologies Grant No. 820419. RC, LAS and FP acknowledge the support from a government grant managed by the Agence Nationale de la Recherche under the Plan France 2030 with the reference "ANR-22-PETQ-0005". LS and GR acknowledge financial support from the PRIN 2022 project "Quantum Sensing and Precision Measurements with Nonclassical States". MM and GMT acknowledge acknowledge funding from the European Union's Next Generation
EU Programme IR0000016 I-PHOQS "Integrated Infrastructure Initiative in Photonic and Quantum Sciences. LS and GMT acknowledge financial support from PNRR MUR Project No. PE0000023-NQSTI

References

[1] A. D. Cronin, J. Schmiedmayer and D. E. Pritchard, Optics and interferometry with atoms and molecules, *Rev. Mod. Phys.* **81**, 1051 (2009).

[2] G. M. Tino and M. A. Kasevich (editors) *Atom Interferometry* (SIF and IOS Press, Bologna, Amsterdam, 2014).

[3] M. S. Safronova, D. Budker, D. DeMille, D. F. J. Kimball, A. Derevianko, and C. W. Clark, Search for new physics with atoms and molecules, *Rev. Mod. Phys.* **90**, 025008 (2018).

[4] R. Geiger, A. Landragin, S. Merlet, and F. Pereira Dos Santos, High-accuracy inertial measurements with cold-atom sensors, *AVS Quantum Sci.* **2**, 024702 (2020).

[5] K. Bongs, M. Holynski, J. Vovrosh, P. Bouyer, G. Condon, E. Rasel, C. Schubert, W. P. Schleich, and A. Roura, Taking atom interferometric quantum sensors from the laboratory to real-world applications, *Nat. Rev. Phys.* **1**, 731 (2019).

[6] A. Gauguet, B. Canuel, T. Lévèque, W. Chaibi, and A. Landragin, Characterization and limits of a cold-atom Sagnac interferometer, *Phys. Rev. A* **80**, 063604 (2009).

[7] F. Sorrentino, Q. Bodart, L. Cacciapuoti, Y.-H. Lien, M. Prevedelli, G. Rosi, L. Salvi, and G. M. Tino, Sensitivity limits of a Raman atom interferometer as a gravity gradiometer, *Phys. Rev. A* **89**, 023607 (2014).

[8] C. Janvier et al., Compact differential gravimeter at the quantum projection-noise limit, *Phys. Rev. A* **105**, 022801 (2022).

[9] L. Pezzè and A. Smerzi, Entanglement, nonlinear dynamics, and the Heisenberg limit, *Phys. Rev. Lett.* **102**, 100401 (2009).

[10] L. Pezzè, A. Smerzi, M. K. Oberthaler, R. Schmied, and P. Treutlein, Quantum metrology with nonclassical states of atomic ensembles, *Rev. Mod. Phys.* **90**, 035005 (2018).

[11] M. Kitagawa and M. Ueda, Squeezed spin states, *Phys. Rev. A* **47**, 5138 (1993).

[12] D. J. Wineland, J. J. Bollinger, W. M. Itano, and D. J. Heinzen, Squeezed atomic states and projection noise in spectroscopy, *Phys. Rev. A* **50**, 67 (1994).

[13] J. Ma, X. Wang, C.P. Sun and F. Nori, Quantum spin squeezing, *Phys. Rep.* **509**, 89 (2011).

[14] A. Louchet-Chauvet, J. Appel, J. J. Renema, D. Oblak, N. Kjaergaard, and E. S. Polzik, Entanglement-assisted atomic clock beyond the projection noise limit, *New J. Phys.* **12** 065032 (2010).
[15] I. Kruse, K. Lange, J. Peise, B. Lücke, L. Pezzè, J. Arlt, W. Ertmer, C. Lisdat, L. Santos, A. Smerzi, and C. Klempt, Improvement of an atomic clock using squeezed vacuum, *Phys. Rev. Lett.* **117**, 143004 (2016).

[16] O. Hosten, et al., Measurement noise 100 times lower than the quantum-projection limit using entangled atoms, *Nature* **529**, 505 (2016).

[17] T. Laudat, et al., Spontaneous spin squeezing in a rubidium BEC, *New J. Phys.* **20**, 073018 (2018).

[18] E. Pedrozo-Penafiel, et al., Entanglement on an optical atomic-clock transition, *Nature* **588**, 414 (2020).

[19] M. Huang et al., Observing Spin-Squeezed States under Spin-Exchange Collisions for a Second, *PRX Quantum* **4**, 020322 (2023).

[20] W. J. Eckner, N. D. Oppong, A. Cao, A. W. Young, W. R. Milner, J. M. Robinson, J. Ye and A. M. Kaufman, Realizing spin squeezing with Rydberg interactions in an optical clock, *Nature* **621**, 734 (2023).

[21] J. M. Robinson, M. Miklos, Y. M. Tso, C. J. Kennedy, T. Bothwell, D. Kedar, J. K. Thompson and J. Ye, Direct comparison of two spin-squeezed optical clock ensembles at the 10⁻¹⁷ level, *Nat. Phys.* **20**, 208 (2024).

[22] L. Salvi, N. Poli, V. Vuletić, and G. M. Tino, Squeezing on Momentum States for Atom Interferometry, *Phys. Rev. Lett.* **120**, 033601 (2018).

[23] S. S. Szigeti et al., High-Precision Quantum-Enhanced Gravimetry with a Bose-Einstein Conden-sate, *Phys. Rev. Lett.* **125**, 100402 (2020).

[24] R. Corgier, L. Pezzè and A. Smerzi, Nonlinear Bragg interferometer with a trapped Bose-Einstein condensate, *Phys. Rev. A* **103**, L061301 (2021).

[25] R. Corgier, N. Gaaloul, A. Smerzi and L. Pezzè, Delta-kick Squeezing, *Phys. Rev. Lett.* **127**, 183401 (2021).

[26] R. Corgier, M. Malitesta, A. Smerzi and L. Pezzè, Quantum-enhanced differential atom interferometers and clocks with spin-squeezing swapping, *Quantum* **7**, 965 (2023).

[27] G. P. Greve, C. Luo, B. Wu, and J. K. Thompson, Entanglement-Enhanced Matter-Wave Interferometry in a High-Finesse Cavity, *Nature* **610**, 472 (2022).

[28] B. K. Malia, Y. Wu, J. Martínez-Rincón and M. A. Kasevich, Distributed quantum sensing with mode-entanglement spin-squeezed atomic states, *Nature* **612**, 661 (2022).

[29] F. Anders, A. Idel, P. Feldmann, D. Bondarenko, S. Loriani, K. Lange, J. Peise, M. Gersemann, B. Meyer-Hoppe, S. Abend, N. Gaaloul, C. Schubert, D. Schlippert, L. Santos, E. Rasel, and C. Klempt, Momentum Entanglement for Atom Interferometry, *Phys. Rev. Lett.* **127**, 140402 (2021).

[30] C. Cassens, B. Meyer-Hoppe, E.M Rasel, and C. Klempt, Entanglement-enhanced atomic gravimeter, arXiv:2404.18668.

[31] A. Andrè, A. S. Sørensen, and M. D. Lukin, Stability of Atomic Clocks Based on Entangled Atoms, *Phys. Rev. Lett.* **92**, 230801 (2004).
[32] B. Braverman, A. Kawasaki, and V. Vuletić, Impact of non-unitary spin squeezing on atomic clock performance, *New. J. Phys.* **20**, 103019 (2018)

[33] L. Pezzè ans A. Smerzi, Heisenberg-Limited Noisy Atomic Clock Using a Hybrid Coherent and Squeezed State Protocol, *Phys. Rev. Lett.* **125**, 210503 (2020).

[34] L. Pezzè ans A. Smerzi, Quantum Phase Estimation Algorithm with Gaussian Spin States, *Phys. Rev. X Quantum* **2**, 040301 (2013).

[35] J. Borregaard and A. S. Sørensen, Near-Heisenberg- Limited Atomic Clocks in the Presence of Decoherence, *Phys. Rev. Lett.* **111**, 090801 (2013).

[36] R. Kaubruegger, D. V. Vasilyev, M. Schulte, K. Hammerer, and P. Zoller, Quantum Variational Optimization of Ramsey Interferometry and Atomic Clocks, *Phys. Rev. X* **11**, 041045 (2021).

[37] C. D. Marciniak, T. Feldker, I. Pogorelov, R. Kaubruegger, D. V. Vasilyev, R. van Bijnen, P. Schindler, P. Zoller, R. Blatt, and T. Monz, Optimal metrology with programmable quantum sensors, *Nature* **603**, 604 (2022).

[38] M. J. Snadden, J. M. McGuirk, P. Bouyer, K. G. Haritos, M. A. Kasevich, Measurement of the Earth’s gravity gradient with an atom interferometer-based gravity gradiometer, *Phys. Rev. Lett.* **81** 971 (1998).

[39] G. T. Foster, J. B. Fixler, J. M. McGuirk, and M. A. Kasevich, Method of phase extraction between coupled atom interferometers using ellipse-specific fitting, *Opt. Lett.* **27**, 951 (2002).

[40] F. Pereira dos Santos, Differential phase extraction in an atom gradiometer, *Phys. Rev. A* **91**, 063615 (2015).

[41] G. Rosi, L. Cacciapuoti, F. Sorrentino, M. Menchetti, M. Prevedelli, and G. M. Tino, Measurement of the Gravity-Field Curvature by Atom Interferometry, *Phys. Rev. Lett.* **114**, 013001 (2015).

[42] B. Barrett, L. Antoni-Micollier, L. Chichet, B. Battelier, T. Lévèque, A. Landragin, and P. Bouyer, Dual Matter-Wave Inertial Sensors in Weightlessness, *Nature Communications*, **7** 2041-1723 (2016).

[43] M. Langlois, R. Caldani, A. Trimeche, S. Merlet, F. Pereira Dos Santos, Differential phase extraction in dual interferometers exploiting the correlation between classical and quantum sensors, *Phys. Rev. A* **96**, 053624 (2017).

[44] R. H. Parker, C. Yu, W. Zhong, B. Estey, and H. Müller, Measurement of the fine-structure constant as a test of the Standard Model, *Science* **360**, 191 (2018).

[45] G. Rosi, G. D'Amico, L. Cacciapuoti, F. Sorrentino, M. Prevedelli, M. Zych, Č. Brukner, and G. M. Tino, Quantum test of the equivalence principle for atoms in coherent superposition of internal energy states, *Nat. Commun.* **8**, 15529 (2019).

[46] L. Hu, N. Poli, L. Salvi, and G. M. Tino, Atom Interferometry with the Sr Optical Clock Transition, *Phys. Rev. Lett.* **119**, 263601 (2017).

[47] L. Salvi, L. Cacciapuoti, G. M. Tino, and G. Rosi, Atom Interferometry with Rb Blue Transitions, *Phys. Rev. Lett.* **131**, 103401 (2023).
[48] E. R. Elliott, et al. Quantum gas mixtures and dual-species atom interferometry in space, *Nature* **623**, 502–508 (2023)

[49] J. B. Fixler, G. T. Foster, J. M. McGuirk, and M. A. Kasevich, Atom Interferometer Measurement of the Newtonian Constant of Gravity, *Science* **5** 74 (2007)

[50] G. Lamporesi, A. Bertoldi, L. Cacciapuoti, M. Prevedelli, and G. M. Tino, Determination of the Newtonian Gravitational Constant Using Atom Interferometry, *Phys. Rev. Lett.* **100** 050801 (2008).

[51] G. Rosi, F. Sorrentino, L. Cacciapuoti, M. Prevedelli, and G. M. Tino, Precision measurement of the Newtonian gravitational constant using cold atoms, *Nature* **510**, 518 (2014).

[52] C. Overstreet, P. Asenbaum, J. Curti, M. Kim, M. A. Kasevich, Observation of a gravitational Aharonov-Bohm effect *Science* **375**, 226 (2022).

[53] S. Dimopoulos, P. W. Graham, J. M. Hogan, M. A. Kasevich, Testing general relativity with atom interferometry, *Phys. Rev. Lett.* **98** 111102 (2007).

[54] G. Varoquaux, R. A. Nyman, R. Geiger, P. Cheinet, A. Landragin, and P. Bouyer, How to estimate the differential acceleration in a two-species atom interferometer to test the equivalence principle, *New J. Phys.* **11** 113010 (2009).

[55] B. Canuel et al., ELGAR-a European Laboratory for Gravitation and Atom-interferometric Research, *Class. Quantum Grav.* **37** 225017 (2020).

[56] S. Dimopoulos, P. W. Graham, J. M. Hogan, M. A. Kasevich, S. Rajendran, Atomic gravitational wave interferometric sensor, *Phys. Rev. D* **78** 122002 (2008)

[57] P. Asenbaum, C. Overstreet, M. Kim, J. Curti, M. A. Kasevich, Atom-Interferometric Test of the Equivalence Principle at the $10^{-12}$ Level, *Phys. Rev. Lett.* **125**, 191101 (2020).

[58] C. Struckmann, R. Corgier, S. Loriani, G. Kleinsteinberg, N. Gox, E. Giese, G. Métris, N. Gaaloul, P. Wolf, Platform and environment requirements of a satellite quantum test of the Weak Equivalence Principle at the $10^{-17}$ level, Phys. Rev. D **109**, 064010 (2024).

[59] D. Schlippert, J. Hartwig, H. Albers, L.L. Richardson, C. Schubert, A. Roura, W.P. Schleich, W. Ertmer, and E.M. Rasel, Quantum test of the universality of free fall, *Phys. Rev. Lett.* **112**, 203002 (2014).

[60] L. Zhou, et al., Test of equivalence principle at $10^{-8}$ level by a dual-species double-diffraction Raman atom interferometer, *Phys. Rev. Lett.* **115**, 013004 (2015).

[61] K. Zhang, et al., Testing the universality of free fall by comparing the atoms in different hyperfine states with Bragg diffraction, *Chin. Phys. Lett.* **37**, 043701 (2020).

[62] B. Barrett, et al., Testing the universality of free fall using correlated atom interferometers, *AVS Quantum Sci.* **4**, 014401 (2022).

[63] D. S. Durfee, Y. K. Shaham, and M. A. Kasevich, Long-Term Stability of an Area-Reversible Atom-Interferometer Sagnac Gyroscope, *Phys. Rev. Lett.* **97**, 240801 (2006)

[64] B. Stray, et al., Quantum sensing for gravity cartography, *Nature* **602**, 590 (2022).
[65] A. W. Young, W. J. Eckner, W. R. Milner, D. Kedar, M. A. Norcia, E. Oelker, N. Schine, J. Ye, and A. M. Kaufman, Half-minute-scale atomic coherence and high relative stability in a tweezer clock, *Nature* **588**, 408 (2020).

[66] X. Zheng, J. Dolde, V. Lochab, B. N. Merriman, H. Li, and S. Kolkowitz, Differential clock comparisons with a multiplexed optical lattice clock, *Nature* **602**, 425 (2022).

[67] M. J. Collett, and G. J. Tee, Ellipse fitting for interferometry. Part 2: experimental realization, *Appl. Opt.* **53**, 7697 (2014).

[68] M. J. Collett, and G. J. Tee, Ellipse fitting for interferometry. Part 3: dynamic method, *J. Opt. Soc. Am. A* **32**, 491 (2015).

[69] J. K. Stockton, X. Wu, and M. A. Kasevich, Bayesian estimation of differential interferometer phase, *Phys. Rev. A* **76**, 033613 (2007).

[70] K. Eckert, P. Hyllus, D. Bruß, U. V. Poulsen, M. Lewenstein, C. Jentsch, T. Müller, E. M. Rasel and W. Ertmer, Differential atom interferometry beyond the standard quantum limit, *Phys. Rev. A* **73**, 013814 (2006).

[71] M. Gessner, L. Pezzè, and A. Smerzi, Sensitivity Bounds for Multiparameter Quantum Metrology, *Phys. Rev. Lett.* **121**, 130503 (2018).

[72] M. Landini, M. Fattori, L. Pezzè, and A. Smerzi, Phase-noise protection in quantum-enhanced differential interferometry *New J. of Phys.* **16**, 113074 (2018).

[73] M. J. Collett, and G. J. Tee, Ellipse fitting for interferometry. Part 1: static methods, *J. Opt. Soc. Am. A* **31**, 2573 (2014).

[74] H.-Y. Zhu, M. Xue, J. Li and M. Deng, Efficacy of differential phase extraction methods in dual atom interferometers under correlated noise conditions, *Results in Physics* **52**, 106855 (2023).

[75] X. Zhang, et al., Dependence of the ellipse fitting noise on the differential phase between interferometers in atom gravity gradiometers, *Opt. Express* **31**, 44102 (2023).

[76] K. Ridley and A. Rodgers, An investigation of errors in ellipse-fitting for cold-atom interferometers, *EPJ Quantum Technol.* **11**, 79 (2024).

[77] B. Yurke, S. L. McCall, J. R. Klauder, SU (2) and SU (1, 1) interferometers *Phys. Rev. A* **33**, 4033 (1986).

[78] F. T. Arecchi, E. Courtens, R. Gilmore, and H. Thomas, Atomic Coherent States in Quantum Optics, *Phys. Rev. A* **6**, 2211 (1972).

[79] An exact analytical solution to $\sigma_z^2|_{\phi=0} = \sigma_z^2|_{\phi=\pi/2}$ can be obtained in the regime $1/N \ll \tau \ll 1/\sqrt{N}$, with $N \gg 1$. We use $\sigma_z^2|_{\phi=0} \approx 1/(2N^2\tau^2) + N^2\tau^4/12$ and $\sigma_z^2|_{\phi=\pi/2} \approx N^3\tau^4/4$, giving Eq. (17).

[80] W. Gander, G.H. Golub, and R. Strebel, Least-squares fitting of circles and ellipses, *BIT* **34**, 558 (1994).

[81] P.L. Rosin, A Note on the Least Squares Fitting of Ellipses, *Pattern Recognition Letters* **14**, 799 (1993).
[82] A. Fitzgibbon, M. Pilu and R. B. Fisher, Direct least square fitting of ellipses, *Proceedings of 13th International Conference on Pattern Recognition 1*, 253 (1996).

[83] R. Halír, J. Flusser, Numerically stable direct least squares fitting of ellipses. *Proceedings of the 6th International Conference in Central Europe on Computer Graphics and Visualization. WSCG '98.* (1998).

[84] W. Gander, Least squares with a quadratic constraint, *Numer. Math.* **36**, 291 (1980).

[85] F. L. Bookstein, Fitting conic sections to scattered data, *Computer Graphics and Image Processing* **9**(1), 56 (1979).

[86] G. Taubin, Estimation of Planar Curves, Surfaces and NonPlanar Space Curves Defined by Implicit Equations, With Applications to Edge and Range Image Segmentation, *IEEE Trans. Pattern Analysis and Machine Intelligence* **13** (11), 1115 (1991).

[87] A.W. Fitzgibbon, and R.B. Fisher, A Buyer's Guide to Conic Fitting, *Proceedings of the British Machine Conference*. BMVC 1995 (1995).

[88] K. Kanatani, Statistical bias of conic fitting and renormalization, *IEEE Transactions on Pattern Analysis and Machine Intelligence* **16**(3), 320 (1994).

[89] V. Bužek, R. Derka, and S. Massar, Optimal Quantum Clocks, *Phys. Rev. Lett.* **82**, 2207 (1999).

[90] J. Hurt, Asymptotic expansions of functions of statistics, *Aplikace matematiky* **21**, 444 (1976).
5. Appendix

Appendix A. One-parameter fit

The one-parameter fit method relies on a simplified algebraic fit that consists to minimize Eq. (23) with

$$
d(z_A, z_B; C_v) = C_{\tau_B}^2 z_A^2 + C_{\tau_A}^2 z_B^2 - 2C_{\tau_A} C_{\tau_B} z_A z_B h - C_{\tau_A}^2 C_{\tau_B}^2 (1 - h^2), \quad (A.1)
$$

where $h = \cos(\delta\phi)$ is the only parameter to be estimated and where $d(z_A, z_B; C_v) = 0$ corresponds to the average ellipse. The estimate $\delta\phi_{est}$ is determined by the condition $\frac{d}{dh} \sum_{j=1}^{n} d^2(z_{A,j}, z_{B,j}) = 0$, and corresponds to solve the cubic polynomial equation given by Eq. (25), where the coefficients are $G_l = (1/n) \sum_{j=1}^{n} g_l(z_{A,j}, z_{B,j})$, $l \in \{0, 3\}$ with

$$
g_0(z_A, z_B) = (C_{\tau_B}^2 z_A^2 + C_{\tau_A}^2 z_B^2 - C_{\tau_A}^2 C_{\tau_B}^2) z_A z_B, \quad (A.2)
$$

$$
g_1(z_A, z_B) = -C_{\tau_A} C_{\tau_B} (C_{\tau_B}^2 z_A^2 + C_{\tau_A}^2 z_B^2 - C_{\tau_A}^2 C_{\tau_B}^2 + 2z_A z_B), \quad (A.3)
$$

$$
g_2(z_A, z_B) = 3 C_{\tau_A}^2 C_{\tau_B}^2 z_A z_B, \quad (A.4)
$$

$$
g_3(z_A, z_B) = -C_{\tau_A}^3 C_{\tau_B}^3. \tag{A.5}
$$

The differential phase estimate is found as

$$
\delta\phi_{\text{est}} = f(G_0, G_1, G_2, G_3) = \arccos[p_{\text{cube}}(G_0, G_1, G_2, G_3)] \quad (A.6)
$$

where $p_{\text{cube}}$ is the real solution of the cubic equation given by the Cardano formula. Since $(z_{A,j}, z_{B,j})$ are independent samples with the same distribution $\mathcal{P}(z_A, z_B | \delta\phi)$, then $g_l(z_{A,j}, z_{B,j})$ are independent samples of $g_l(z_A, z_B)$ with identical distributions. In addition, the estimator of $\delta\phi_{\text{est}}$ being a function of sample moments, Eqs. (26) and (27) remains valid in the limit of $n$ large [90]. We use the identities $\langle G_l \rangle = \langle g_l \rangle$ and $\text{Cov}(G_j, G_l) = (1/n) \text{Cov}(g_j, g_l)$ to the calculate the first and second moments of $g_l$. The mean values $\langle g_l \rangle$ are given by

$$
\langle g_l \rangle = \frac{1}{2\pi} \int_0^{2\pi} d\phi_{\text{cn}} \langle \psi_A^{\text{out}}(\phi_A) | \langle \psi_B^{\text{out}}(\phi_B) | \hat{g}_l | \psi_A^{\text{out}}(\phi_A) \rangle | \psi_B^{\text{out}}(\phi_B) \rangle, \quad (A.7)
$$

where $\phi_A = \phi_{cn} + \delta\phi/2$ and $\phi_B = \phi_{cn} - \delta\phi/2$. Here, $\hat{g}_l$ is the operator equivalent of $g_l$ given by
$\hat{g}_l = g_l(\hat{z}_A, \hat{z}_B)$ which does not require symmetrization since $\hat{z}_A$ and $\hat{z}_B$ commute. Analogous
formulas can be derived for the second moments.

Appendix A.1. Squeezed-spin-state moments

Accessing the first and second moment of $\hat{g}_l$ respectively requires to evaluate the the moment
of $\hat{z}$ up to third and sixth orders. The $k$-th order is defined as $\langle \psi^{\text{out}}(\phi) | (\hat{J}_z)^k | \psi^{\text{out}}(\phi) \rangle = \langle \psi^{\text{Coh}} | (\hat{J}_z(\phi, v, \tau))^k | \psi^{\text{Coh}} \rangle$. Here we used the relation $|\psi^{\text{out}}(\phi)\rangle = U|\psi^{\text{Coh}}\rangle$ where $U = \exp(-i\phi\hat{J}_y)\exp(-iv\hat{J}_x)\exp(-i\tau\hat{J}_z^2)$ and $U^\dagger(\hat{J}_z)^k U = (U^\dagger\hat{J}_z U)^k$. On the overall, we have

$$
\hat{J}_z(\phi, v, \tau) = \cos(\phi) [\cos v \hat{J}_z + \sin v \hat{J}_y(\tau)] - \sin(\phi) \hat{J}_x(\tau), \quad (A.8)
$$
with $\hat{J}_x(\tau) = [\hat{J}_+ \exp_+(\tau) + \exp_-(\tau) \hat{J}_-]/2$, $\hat{J}_y(\tau) = [\hat{J}_+ \exp_+(\tau) - \exp_-(\tau) \hat{J}_-]/(2i)$, $\hat{J}_\pm = \hat{J}_x \pm i\hat{J}_y$ and $\exp_\pm(\tau) = \exp[\pm 2i\tau(\hat{J}_z + 1/2)]$. The expectation value of Eq. (A.8) leads to Eq. (9) when rescaling by $2/N$. Similarly, order 2 leads to Eq. where the two extrema $\sigma_z^2|_{\phi=0}$ and $\sigma_z^2|_{\phi=\pi/2}$ correspond to the variance on $|\psi^{\text{Coh}}\rangle$ of $\cos \nu \hat{J}_z + \sin \nu \hat{J}_y(\tau)$ and of $\hat{J}_x(\tau)$, respectively, properly rescaled by $4/N^2$ (12).

## Appendix A.2. Fit bias

The method outlined above also allows to find an expression for the $\langle g_l \rangle$ and the fit bias $\mathcal{B}(\delta\phi_{est})$. The calculation further requires the evaluation of the integral in Eq. (A.7). Assuming $\tau_A = \tau_B \equiv \tau$, we obtain:

$$ \langle g_0 \rangle = C_{\tau}^3 \left[ -\frac{C_{\tau}^3}{2} + \frac{3}{4N} C_{\tau} + \frac{3\Sigma_1 + \Sigma_2}{4} \right] h, \qquad (A.9) $$

$$ \langle g_1 \rangle = - \left\{ \frac{C_{\tau}^4}{4} + \frac{C_{\tau}^2}{2} (5\sigma_0^2 + 3\sigma_{\pi/2}^2) + \frac{(\sigma_0^2 - \sigma_{\pi/2}^2)^2}{4} - 2\sigma_0^2\sigma_{\pi/2}^2 \right. \\ \left. + \left[ \frac{C_{\tau}^4}{2} - C_{\tau}^2 (\sigma_0^2 - \sigma_{\pi/2}^2) + \frac{(\sigma_0^2 - \sigma_{\pi/2}^2)^2}{2} \right] h^2 \right\} C_{\tau}^2, \qquad (A.10) $$

$$ \langle g_2 \rangle = \frac{3}{2} C_{\tau}^6 h, \qquad (A.11) $$

$$ \langle g_3 \rangle = -C_{\tau}^6. \qquad (A.12) $$

In these equations, $\sigma_0^2$ and $\sigma_{\pi/2}^2$ are shorthand notation for $\sigma_z^2|_{\phi=0}$ and $\sigma_z^2|_{\phi=\pi/2}$. In the limit $N \gg 1$, we have $\Sigma_1 = (C_{3\tau} + 3C_\tau)/4 + 3C_\tau/N$ and $\Sigma_2 = -3[\cos^2(\nu) \sin^2(\tau) C_\tau + \sin^2(\nu)(C_{3\tau} - C_\tau)/4 - \sin(\nu) \cos(\nu) \sin(2\tau) C_{2\tau}]$. A Taylor expansion of $f(\langle g_0 \rangle, \langle g_1 \rangle, \langle g_2 \rangle, \langle g_3 \rangle)$ around the unbiased point defined by $\tau = 0$ and $N \to \infty$, lead to a simplified expression for the bias. Defining $\langle g_l \rangle_\infty = \lim_{N \to \infty} \langle g_l \rangle|_{\tau=0}$, we get $\langle g_0 \rangle_\infty = h/4$, $\langle g_1 \rangle_\infty = -(1/4 + h^2/2)$, $\langle g_2 \rangle_\infty = 3h/2$ and $\langle g_3 \rangle_\infty = -1$. The solution of the cubic equation is then the unbiased phase difference, i.e. $f(g_\infty) = \delta\phi$. At first order, Eq. (29) is obtained by:

$$ \mathcal{B}(\delta\phi_{est}) = \langle\delta\phi_{est}\rangle - \delta\phi \approx \sum_{l=0}^{3} \left(\frac{\partial f}{\partial g_l}\right)_{g=g_{\infty}} (\langle g_l \rangle - \langle g_l \rangle_{\infty}), \qquad (A.13) $$

where $\partial f / \partial g_l = -[1 - (p_{cube})^2]^{-1/2} (\partial p_{cube} / \partial g_l)$ such as $(\partial f / \partial g_l)_{g=g_\infty} = -4(1 - h^2)^{-1/2}(1 + 2h^2)^{-1} h^l$.

## Appendix B. Dependence on the number of fitted points

In Fig. B1, we show the dependence of $\mathcal{B}(\delta\phi_{est})$ (empty symbols) and $\sigma_{\delta\phi_{est}}$ (filled symbols) with respect to the number of points sampled $n$ in the ellipse for two coherent states, i.e. $\tau = 0$, (circles) and two squeezed states of strength $\tau^*$ (stars). While the standard deviation shows the expected scaling $\sim 1/\sqrt{n}$, the bias saturates for sufficiently large $n$ to a value
<figure><img src="image_8.png" /><figcaption>**Figure B1.** Bias $\mathcal{B}(\delta\phi_{est})$ and standard deviation $\sigma_{\delta\phi_{est}}$ as functions of $n$. Results are relative to $N = 500$ atoms, a differential phase shift $\delta\phi = \pi/16$ and are obtained by means of an algebraic fit with trace constraint used in Figs. 3 to 6. The dashed line corresponds to the SQL for the considered differential scheme and given by Eq. 14.</figcaption></figure>

dependent on $\delta\phi$, $\tau$ and $N$. Despite a different behavior for $\tau = 0$ and $\tau^*$, we note that for $n \ge 500$ leads to saturation. In the paper we always used $n = 1000$. The situation where the overall fit error is dominated by the standard deviation can be realized with squeezed states on a wider range of $n$ values. Based on the scaling laws determined in Sec. 3.3 and in the limit $n \to \infty$, the condition $\mathcal{B}(\delta\phi_{est}) < \sigma_{\delta\phi_{est}}$ translates to (i) $n < N$ in the case of coherent states, as $\mathcal{B}(\delta\phi_{est}) \sim 1/N$ and $\sigma_{\delta\phi_{est}} \sim 1/\sqrt{nN}$, and to (ii) $n < N^{7/3}$ for squeezed states with $\tau = \tau^*$, as $\mathcal{B}(\delta\phi_{est}) \sim 1/N^{4/3}$ and $\sigma_{\delta\phi_{est}} \sim 1/(\sqrt{n}N^{1/6})$.