Quantum-enhanced differential atom interferometers and clocks with spin-squeezing swapping

Thanks to common-mode noise rejection, differential configurations are crucial for realistic applications of phase and frequency estimation with atom interferometers. Currently, differential protocols with uncorrelated particles and mode-separable settings reach a sensitivity bounded by the standard quantum limit (SQL). Here we show that differential interferometry can be understood as a distributed multiparameter estimation problem and can benefit from both mode and particle entanglement. Our protocol uses a single spin-squeezed state that is mode-swapped among common interferometric modes. The mode swapping is optimized to estimate the differential phase shift with sub-SQL sensitivity. Numerical calculations are supported by analytical approximations that guide the optimization of the protocol. The scheme is also tested with simulation of noise in atomic clocks and interferometers.


France
Thanks to common-mode noise rejection, differential configurations are crucial for realistic applications of phase and frequency estimation with atom interferometers. Current differential protocols with uncorrelated particles and modeseparable settings reach a sensitivity bounded by the standard quantum limit (SQL). Here we show that differential interferometry can be understood as a distributed multiparameter estimation problem and can benefit from both mode and particle entanglement. Our protocol uses a single spin-squeezed state that is modeswapped among common interferometric modes. The mode swapping is optimized to estimate the differential phase shift with sub-SQL sensitivity. Numerical calculations are supported by analytical approximations that guide the optimization of the protocol. The scheme is also tested with simulation of noise in atomic clocks and interferometers.
Atom interferometers (AIs) are extraordinarily sensitive devices [1][2][3] that find key applications in fundamental and applied science [4,5]. Thanks to the coupling of atoms with external forces and light fields, AIs allows for ultraprecise measurements in inertial sensors [6] and atomic clocks [7,8].
Absolute measurements with AIs are limited by the onset of systematic effects and/or phase noise, inherent to the interferometer itself and to the environment. Therefore, experiments aiming at precision measurements often benefit of a differential configuration. In this case, two interferometers (A and B), working simultaneously, acquire a total phase shift θ A,B = φ A,B + φ pn , respectively, where φ A,B are the signals and φ pn encompasses systematic effects or/and stochastic phase noise that are common to both sensors. Correlations between the output signals of the two interferometers can be used to infer φ A − φ B [9-14], while cancelling the common-mode noise. Differential AIs configurations have been exploited for gravity gradiometry [15][16][17][18][19], gravity-field curvature measurements [20], relativistic geodesy [21], the measurements of the gravitational constant [22], to probe the law of gravitation through tests of the universality of free-fall [23][24][25][26][27], and are also expected to find applications in the detection of gravitational waves [28][29][30][31]. Furthermore, correlation spectroscopy exploits differential frequency measurements. This technique is based on the simultaneous interrogation of two atoms [32,33] or two atomic ensembles with the same laser: it allows for the cancellation of laser fluctuations and is immune to dead times. Differential frequency comparison in atomic clocks has been exploited for the measurement of the gravitational redshift [34][35][36] by taking advantage of long interrogation times.
Current differential AIs use independent interferometers and uncorrelated atoms: a distributed sensing configuration commonly indicated as modeseparable and particle-separable (MsPs) [37]. In this case, the sensitivity of the differential measurement is inherently limited by the standard quantum limit (SQL) [15,22,[37][38][39][40]. This bound is given by ∆ 2 (θ A − θ B ) SQL = 4/N , where N is the total number of particles [37,41]. A promising research goal is to devise feasible schemes that reduce the quantum noise in interferometric measurements and, in particular, overcome the SQL with the use of quantum resources [42,43].
The experimental investigation on entanglementenhanced AIs has mainly focused on the generation and the manipulation of spin-squeezed states [42]. Spin squeezing in momentum modes, suitable for inertial sensing with AIs, can be created through atomatom interaction in Bose-Einstein condensates [44][45][46] or atom-light interaction in an optical cavity [47,48]. Recent experimental achievements include the creation of useful entangled states through the transfer of squeezing generated in internal states into well-defined and separated external momentum modes [49], the demonstration of long-lived spinsqueezed states [50], and the implementation of a full spatial interferometer scheme with a direct observation of sub-SQL sensitivity [48]. Furthermore, entanglement-enhanced atomic clocks using spinsqueezed states have been demonstrated [51][52][53][54].
In this manuscript, we propose an efficient method for quantum-enhanced differential phase estimation that exploits spin-squeezing. A spin squeezed state is first generated between the two modes of the interferometer A, see Fig. 1(a), for instance via one-axis twisting [55]. The squeezed state is then shared with the modes of interferometer B via a linear transformation -namely an optimized beam-splitter -and mixed with a coherent spin state. We indicate this operation as mode-swapping (MS): it generates entanglement between the modes of the two interferometers, thus realizing a mode-entangled and particleentangled (MePe) distributed quantum sensor [37]. This scheme extends to atomic systems a protocol proposed in Ref. [56], which uses optical squeezedvacuum and coherent states in a network of Mach-Zehnder sensors. With ab-initio numerical calculations supported by an analytical approach, we i) predict the possibility to overcome the SQL with a substantial sensitivity gain, up to ∆ 2 (θ A − θ B ) = O(1/N 3/2 ), and ii) demonstrate that the protocol benefits of common-noise suppression, while using standard measurements and estimation. The sensitivity scaling O(1/N 3/2 ) is due to the one-axis-twisting method and can be pushed to the Heisenberg limit by using other spin-squeezing generation schemes. It should be emphasized that, with the same resources (namely a coherent and a single spin-squeezed state) in a mode-separable and particle-entangled (MsPe) setting [37], see Fig. 1(b), it is only possible to overcome the SQL by a unit factor. The practical advantage of our AI protocol is that it fully exploits a single spin-squeezed state for differential measurements: we can thus envisage the realization of high sensitivity enhancement within state-of-the-art capabilities [42, 57, 58], with a variety of applications.

Results
In the following, we first introduce the interferometer scheme of Fig. 1(a), see Sec. 1.1. The sensitivity of the device is then evaluated using an analytical approach, in Sec. 1.2, as well as numerical calculations for collective spin systems, in Sec. 1.3. Finally, in Sec. 1.4, we provide explicit examples of differen-tial phase and frequency estimation in the presence of common noise, discussing the advantages with respect to mode-separable approaches.

Interferometer scheme
The differential scheme of Fig. 1(a) starts with N atoms prepared in mode a, while the other modes (b, c and d) are empty at this stage. A 50-50 beam splitting pulse between modes a and d [yellow box in Fig. 1(a)] generates a binomial coherent particle distribution, We then apply a two-mode spinsqueezing transformation to generate entanglement between particles in modes a and b [green box in Fig. 1(a)]. We consider one-axis twisting dynamics [55] followed by an opportune rotation of the state, In Eq. 2, the term (Ĵ {a,b} x ) 2 is associated to a non-linear evolution due to atom-atom interaction in a Bose-Einstein condensate [44][45][46] or atomlight interaction in a cavity [47,48,59], and generates a two-mode spin-squeezing dynamics, where τ E denotes the effective squeezing time. The term exp{−iν EĴ {a,b} z } rotates the state in the Bloch sphere by the angle ν E . It is used to optimize the orientation of the quadrature of the quantum state to be sensitive to a phase measurement. Experimentally, such transformation can be generated by manipulating the phase of the laser coupling mode a and b. Here,Ĵ are collective pseudo-spin operators for the modes a and b (c and d). The operatorâ (b,ĉ andd) annihilates a particles in mode a (b, c, and d), whileâ † (b † ,ĉ † andd † ) is the corresponding bosonic creation operator.
The key ingredient of our proposal is the MS [blue box in Fig. 1(a)], corresponding to a linear coupling between modes b and c. We indicate as |ψ inp = Figure 1: Differential interferometry within (a) mode-entangled particle-entangled (MePe) and a (b) mode-separable particle-entangled (MsPe) settings. Panel (a): A single input atomic ensemble of N atoms is initially coherently split between modes a and d by a first π/2 pulse (yellow box). A spin-squeezed state is generated in modes a and b via one-axis twisting evolution of strength τ E and orientation angle ν E (green box, see Eq. 2). The MS (blue box) distributes the spin squeezing between the modes b and c through the parameters ϑ MS and ϕ MS representing the strength and the phase of the coupling laser (see Eq. 3). Two separated Mach-Zehnder interferometers (orange boxes) are used to estimate θ A and θ B . The inset shows a schematic of the quantum state at different stages of the protocol on the (x b + x c )-(p b + p c ) quadrature plane, obtained within a Holstein-Primakoff transformation (see text). The distribution corresponding to the state |ψ 2 has variances (1 + e −2r )/2 and (1 + e +2r )/2 along the x b + x c and p b + p c direction, respectively (green dashed line), where r = N τ E /4. These are wider than the isotropic vacuum variance equal to 1 (dotted black line). In contrast, the state after MS, |ψ inp has variance (blue solid line) e −2r along x b + x c and e +2r along p b + p c , below the vacuum limit. Panel (b): Two separated interferometers, using spin squeezed states of N/2 atoms (τ A,B S indicates the squeezing strength -green boxes). In both schemes, the θ A and θ B are evaluated locally by counting the number of atoms at each output detectors.
U MS |ψ 2 the quantum state at this stage, witĥ where ϑ MS and ϕ MS are, respectively, the strength and the phase of the laser coupling mode b an c. As schematically shown in the inset of Fig. 1(a) and explained in more details below, the MS swaps the spin squeezing from the modes a and b (where the squeezing is generated) to a combination of suitable interferometer's modes. Without the MS, e.g. coupling strength ϑ MS = 0, the spin squeezing generated in the a and b modes is only useful to estimate θ A . In this study, the combination of a non-linear dynamic of strength τ E with an adequate orientation of the four-mode state provided by ν E and ϕ M S enables the creation of a quantum-enhanced sensitive states for differential measurements.
The preparation of |ψ inp is followed by two independent Mach-Zehnder interferometers, indicated as A (for modes a and b) and B (for modes c and d). The scheme ends with the measurement of the number of particles in each output port, the output state being |ψ out = exp(−iθ AĴ For two uncoupled interferometers, we have Γ AB = 0 and Eq. (4) recovers the sum of uncertainties In particular, this opens the possibility to overcome the SQL for the estimation of the differential phase. Through this paper, quantifies the sensitivity gain over the SQL [41].

Analytical model
The analytical description of the interferometer is based on two approximations, see App. A.1 for details. First, we write the state |ψ 1 , Eq. (1), as . This approximation assumes N 1, see App. A.1. Experimentally, one can easily generate Bose-Einstein condensate with N approximately ranging from 10 2 to 10 6 atoms such that the condition N 1 is easily satisfied. The phases of the two coherent states in modes a and d are equal and set to zero without loss of generality. Furthermore, we describe the one-axis twisting transformation by using a Holstein-Primakoff transformation replacing the operatorsâ with the classical numberâ ∼ N/2. We obtain where |S(ξ) b =Ŝ b (ξ)|0 b is a single-mode squeezed-vacuum state andŜ b (ξ) = e (ξ * bb −ξb †b † )/2 is the squeezing operator [61] with squeezing parameter ξ = re iϕ . Here, r = N τ E /4, ϕ = π/2 + ν E and n s = sinh 2 r is the average number of particles in the state. Equation (7) is expected to be accurate if the squeezing dynamics is applied for a sufficiently short time such that the mode a remains essentially undepleted, namely N/2 n s (corresponding to τ E 8/N ). Taking the quadrature operatorx b (λ) = (be −iλ +b † e iλ )/ √ 2, the state |S(ξ) b is squeezed along the direction λ = π/4 + ν E /2 in the quadrature plane, namely ∆ 2x b (π/4 + ν E /2) = e −2r /2, and anti-squeezed along the orthogonal direction, , we predict a spin-squeezing parameter (8) This expression agrees with the exact result [55], up to the third order in N τ E [62], namely e −2r = 3 . Finally, the MS operation consists of a linear coupling between modes b and c. It is described by the unitary matrix withÛ MS given by Eq. (3) [63]. The state at the input ports of the two interferometers writes (10) where the unitary transformation,Û MS , applies to the state |S(ξ) b |0 c .
Based on the approximations discussed above, we have thus established a relation, in the appropriate regime, between an atomic system using spinsqueezed states and a quantum optical system using squeezed vacuum states, the MS playing the role of a quantum circuit. We can thus apply the results of Ref.
[56], which discussed the sensitivity of a distributed quantum sensor of Mach-Zehnder interferometers using coherent and single-mode squeezed states. Taking Eq. (10) as sensor's input state, we can obtain an analytical expression for Eq. (4) [56], which is lengthy and explicitly given in App. A.2. In particular, it reads where Q ≥ 0 is a second-degree polynomial function in the variables cot θ A and cot θ B . In the rest of this section, we restrict to the optimal working point θ A = θ B = π/2, where Q = 0.
To clarify the role played by the MS for highly sensitive differential estimation, we rely again on a Holstein-Primakoff transformation. This time, it is applied to the operatorsâ andd relative to the output modes, that are replaced by the classical numberŝ a ∼ e −iν E /2 N/2 andd ∼ N/2, respectively. As detailed in App. A.3, Eq. (5) rewrites in terms of the quadrature variance According to Eq. (11), the phase uncertainty of the differential measurement is equivalent to the variance of the quadrature operatorx b +x c . Without mode swapping, namely when considering the state c . Hence, whatever amount of squeezing we generate in the x b quadrature, it will provide limited benefit to the sensitivity of the differential estimation, which is indeed dominated by the vacuum fluctuations in mode c: ∆ 2 (x b +x c ) = (1+e −2r )/2. In contrast, it is possible to reduce the variance ofx b +x c below the vacuum fluctuations by exploiting the MS transformation Eq. (9). The smallest quadrature variance -and, thus, the best sensitivity -is clearly reached when both the rotation angle ν E in Eq. (2) and the parameters of Eq. (9) are optimized. We find ν E = −π/4 + lπ + 4mπ and δ cb = −π/8 + lπ/2, with l, m ∈ Z, and |u bb | = |u cb | = 1/ √ 2 [see App. B.1]. In the following, we choose l = −1, m = 1, getting ν E = 11π/4 and δ cb = −5π/8. Since δ cb = ϕ M S − π/2, this corresponds to an optimal MS operation in Eq. (3) with ϕ MS = −π/8, and the condition |u bb | = |u cb | = 1/ √ 2 translates to ϑ MS = π/2. It is shown in App. A.3 that such optimal choice of parameters yields ∆ 2 (x b +x c ) = e −2r : the squeezing initially generated in the quadrature planex b -p b has been "swapped" to thex b +x c quadrature [see the inset of Fig. 1

(a)] by means of a suitable MS operation.
According to Eq. (11),

which predicts a sensitivity gain
below the SQL. Therefore, through MS operation, we are able to convert a squeezing useful for singleparameter estimation to a squeezing useful to beat the SQL in the differential measurement. Without relying to a Holstein-Primakoff transformation, we obtain which is valid for N n s . The detailed derivation of Eq. (13) can be found in App. B.1. For sufficiently small values of τ E , such that N e 2rn s ≈ 4n 2 s , the termn s /N in Eq. (13) can be neglected, recovering the gain factor Eq. (12). Instead, at the optimal working pointn s = √ N /2, which still satisfies the conditionn s N for N 1, Eq. (13) predicts a maximum gain corresponding to the minimum uncertainty

Numerical results
In the following, we relieve the approximations used in the previous section and perform numerical computation of collective spin-coupling operations. In particular, we provide a detailed analysis of Eqs. (4)-(6) as well as a comparison with the MsPe scheme of Fig. 1(b). For the sake of clarity, we distinguish the squeezing strength τ E of the MePe strategy of Fig. 1(a), to the one of MsPe strategy, τ A S and τ B S of Fig. 1(b). The squeezing strengths are normalized over the value τ ref = 1.2(N/2) −2/3 providing the optimal amount of spin squeezing within the one-axistwisting scheme [55] with N/2 atoms.

Gain and bandwidth
In Fig 2( while, for the MS operation, we have used the optimized laser strength, ϑ MS = π/2, corresponding to a 50-50 beam-splitter, with a fixed (not-optimized) laser phase, ϕ MS = π/2, in Eq. (3). As expected from Eqs. (11)-(13), we observe a substantial amount of gain over the SQL, despite the phase of the MS, ϕ MS , is not optimized here (see discussion below). In addition, at the optimal working point (as a function of τ E ), we obtain a scaling G 2 max ∼ √ N , in agreement with the analytical predictions of Eq. (14). The dot-dashed line in the figure shows the MsPe case where the MS is not applied (ϑ MS = 0). This situation is equivalent to τ A S = τ E and τ B S = 0 in Fig. 1(b): only the estimation of θ A benefits from squeezing (∆ 2 θ A < 2/N ), while θ B is estimated with uncertainty ∆ 2 θ B = 2/N . In this case, we have only a modest gain in the estimation of θ A − θ B reaching, at best, approximately a factor 2 over the SQL (horizontal dotted line, Besides the sensitivity gain, another relevant figure of merit is the relative width of the phase domain We indicate this quantity as the bandwidth, R, defined as  values of N , is shown in Fig. 2(b). The different lines correspond to the same parameters as the ones in Fig. 2(a). Although the analytical predictions (shown as by dashed lines) preserve the qualitative trend, in this case, they do not reproduce the numerical findings. The difference here is due to the slight difference of the statistic in the initial states due to the approximations made: details can be found in Appendix A.1 and C. Overall, we see that increasing the number of particles from N = 100 to 2000 leads to a larger optimal sensitivity gain, and more surprisingly, to an improved robustness. The MS clearly improves the performance of the differential interferometer, when compared to the mode-separable case, using a single spin-squeezed state.
Comparing the two panels in Fig 2, we observe that the optimal gain point does not corresponds to the point of maximal bandwidth. To quantify the tradeoff between gain and sensitivity bandwidth, we define the average gain in a squared box of total length 2Λ pn centered a φ A = φ B = π/2, where Λ pn denotes the maximum phase noise amplitude common to both interferometer shifting the detection from the optimal working point, θ A,B = π/2, to a non-optimal one. Such situation can be obtained experimentally due to common vibration of the plateform for instance or other bias phase-shift changing randomly shot-to-shot the working point of the dif-  (a), τ E is chosen such that G 2 (π/2, π/2) = 10 dB; in panel (b), we considered the value of τ E that maximizes G 2 (π/2, π/2). These two configurations are respec- The average gain decreases with Λ pn but it remains above 1 for all values of Λ pn . We also compare the MePe strategy with the MsPe one (green dotdashed line) using two squeezed states τ A S = τ B S = τ S with (a) τ S /τ ref = 0.24 corresponding to a maximum gain of 10 dB, and (b) the optimal τ S /τ ref = 1 value. In the case (a), the MsPe shows a slightly higher average gain than the MePe scheme, but the former uses two squeezed states instead of one. In (b), the average gain of the MsPe case is higher for small Λ pn , but it quickly degrades with Λ pn as a consequence of the small bandwidth of highly spinsqueezed states. In contrast, the MePe strategy provides a higher average gain for large values of Λ pn . Finally, taking the total squeezing time as a resource, we compare the MePe strategy with the MsPe one with τ A S = τ B S = τ E /2. The corresponding average gain is shown as the black dashed line in Fig 3(a) and (b). In the optimal-squeezing case, Fig 3(b), the MsPe slightly outperforms the mode-entangled case for sufficiently large Λ pn . It should be noticed, nevertheless, that higher dynamical range could be obtained at Λ pn = π/2 for lower value of τ E in the MePe strategy (not shown in the graph).
For further clarity, we plot, in (colored regions correspond to G 2 ≥ 1). Overall, these results show that the mode swapping method allows for a good trade-off between high sensitivity gain for small common phase shift and high bandwidth for large common phase shift, at the expense of generating only a single entangled state.

Optimal mode-swapping parameters
The numerical results shown in the Figs. 2 and 3 have been obtained in the case of a 50-50 beam splitter, i.e. ϑ MS = π/2, for a specific phase value of the MS operation, namely ϕ MS = π/2. We are now interested in the optimization of ϕ MS , still keeping ϑ MS = π/2. In Fig. 4(a), we plot G 2 (π/2, π/2) as a function of ϕ MS for two different squeezing strength τ E and for optimized rotations of the squeezing strength, ν E . Here, we observed a π/2-periodic sinusoidal behavior of the sensitivity gain (predicted analytically in App. B.1). Comparing the two curves, we note that the optimal sensitivity gain is obtained for slightly different optimal phase value ϕ MS which tends, in the regime of low squeezing, to ϕ MS = −π/8, a result anticipated in the discussion of the analytical model, Sec. 1.2.
In Fig. 4(b) we show G 2 (π/2, π/2) (color plot) as a function of the orientation of the spin squeezed state,   Fig. 4(b) should be compared with Fig. 6 in App. B.2. Two conclusions can be drawn here: i) the optimal orientation of the squeezed state and the phase of the mode swapping are linearly linked to each other, and ii) the optimal values of ν E and ϕ MS are robust against experimental imperfection where phase shift as big as π/8 still enable large sensitivity gain.

Differential phase estimation
To illustrate our results, we consider the estimation of the differential phase shift φ A −φ B in the presence of a common-mode dephasing noise. The phase shifts in the two interferometers of Fig. 1 are θ A = φ A + φ pn and θ B = φ B + φ pn , respectively, where φ pn is a stochastic phase noise (fluctuating randomly from shot to shot) with a Gaussian distribution centered to zero and with width σ pn . We use the method of  Fig. 2(a)]. The strong anisotropy of the distribution is due to the correlations between θ axis, as desired. As a reference, the black circle in Fig. 5(a) has a SQL radius. The quantity of interest, φ A − φ B , is estimated by taking θ B . In Fig. 5(b), we plot the uncertainty ∆ 2 (φ A − φ B ) as a function of the noise width σ pn . As we see, the uncertainty is well below the SQL (solid black line), up to values of σ pn /π ≈ 0.2, in this case. The increase of ∆ 2 (φ A − φ B ) is directly associated to the finite bandwidth of the spin-squeezed state due to the unavoidable bending in the Bloch sphere. For comparison, we also plot the sensitivity obtained with a mode-separable strategy with i) two coherent spin states (τ A S = τ B S = 0, black circles); ii) a coherent state and an optimal spin-squeezed state ( τ A S /τ ref = 1 and τ B S = 0, blue triangles); and iii) two optimal spin-squeezed states (τ A S = τ B S = τ ref , green diamonds). The simulations of the differential phase estimation confirm the advantage of the MePe strategy over the MsPe one when using the same resources.

Differential frequency estimation
We consider here relative frequency estimation within an atomic clock scheme. In this case, we have two atomic ensembles interrogated by the same local oscillator (LO). The phase shift in the Ramsey interferometer A and B, accumulated during the interrogation time T , is where ω A,B is the atomic frequency of the atoms in the interferometer A and B, respectively, assumed time independent. The offset π/2 in Eq. (17) places the phase shifts around the optimal working point.
The accumulated frequency noise due to LO fluctuations is common mode and cancels when evaluating . Within a MsPs strategy, using two independent interferometers with N uncorrelated atoms in total, the fractional frequency variance is given by the standard quantum limit (corresponding to the smallest quantum Cramer-Rao bound for the MsPs strategy [37]), where T tot is the total effective interrogation time, given by a multiple of the Ramsey time and ω 0 is a reference value, given by ω 0 = (ω A + ω B )/2, for instance. Using the MePe scheme, we can overcome Eq. (18) for relatively short interrogation times. We find for an optimal squeezing time. Equation (19) holds for interrogation times T limited by the coherence of the LO, namely when θ A and θ B remain sufficiently close to the optimal working point θ A,B = π/2.  Fig. 5(b). In particular, the black circles are obtained for the MsPs setting with two coherent spin state of N/2 particles each. The red dotted line is Eq. (18). The red stars show the results of our MS scheme. For sufficiently short T , our scheme follows Eq. (19) and provides a substantial advantage with respect to the SQL: in other words we can reach the same absolute sensitivity as the fully separable case but with a O(1/ √ N ) decrease of averaging time, which can be an important advantage for some applications. The figure assumes that the atomic coherence lifetime is longer than the LO coherence time. In the opposite limit, the bending of the variance happens at interrogation shorter times, in the regime of validity of Eqs. (18) and (19). In this case, the mode-entangled strategy can improve the overall absolute fractional frequency variance by a factor √ N . Differently from LO stabilization in atomic clock, differential frequency estimation does not suffer the Dick effect associated to dead times [65], even when using squeezed states.

Possible generalizations
In the manuscript, we have focused on the use of spin-squeezed states generated via one-axis twisting dynamics. However, the protocol can be immediately extended to various squeezing procedures. In fact, within a Holstein-Primakoff transformation, we have seen that the differential sensitivity is directly related to the quadrature variance ∆ 2 (x b + x c ), see Eq. (11): any state that realizes ∆ 2 (x b + x c ) < 1 can be used to reduce the differential uncertainty below the SQL. The case of cavity squeezing and two-axis counter-twisting are briefly discussed in App. D. An interesting case is that of the two-mode squeezed vacuum, |ψ = e −ir(b †ĉ † +bĉ) |0 b |0 c , that has been realized in atomic system via spin-mixing dynamics in spinor Bose-Einstein condensates [66][67][68]. This state is characterized by ∆ 2 (x b +x c ) = 2e −2r [61] and thus satisfies ∆ 2 (x b +x c ) < 1 without requiring the MS operation, provided that r > (log 2)/2 ≈ 0.35. Interestingly, the inequality r > (log 2)/2 is also the condition for the two-mode squeezed vacuum state to fulfills the Einstein-Podolsky-Rosen (EPR) criterion [66,[69][70][71]. Recently, EPR entanglement has been proposed as a resource for optical differential homodyne measurements [72,73]. In contrast, within our MS scheme, any single-mode squeezed states with positive value of the squeezing parameter, r > 0 (or, equivalently, any one-axis twisted state with arbitrary small spin-squeezing), can realize a sub-SQL sensitivity. As single mode squeezed states do not fulfill the EPR criterion, we thus conclude that EPR entanglement is a useful resource for differential phase estimation with sub-SQL sensitivity (in agreement with Refs. [72,73]), although it is not a strictly necessary resource.
Besides differential phase estimation, our scheme can also be extended to the sub-SQL estimation of an arbitrary linear combination of the two phase shifts, Finally, it is worth emphasizing that the MePe protocol studied in this manuscript can be combined with existing protocols to extend the squeezed state sensitivity bandwidth, for instance using adaptive schemes [74,75], nonlinear [76,77] or nondemolition [78][79][80] measurements.

Possible experimental implementation with Bose-Einstein condensates
On the overall, the input four-mode entangled state proposed in Fig. 1a could be experimentally generated with a suitable combination of Raman and Bragg diffraction, on ground fountain experiment [26,81] or in space environment [82] for example. Based on recent experimental development and new methods in the field, we highlight a possible experimental sequence with Bose-Einstein condensate (BEC).
Considering the mode a ≡ | ↑, 0 , b ≡ | ↓, 2 k , c ≡ | ↓, 0 and d ≡ | ↑, 2 k with 2 k the momentum transfer, the four-mode state can be engineer as follow: (i) a BEC, initially in mode a is diffracted onto a quantum superposition of same internal state with different external momentum, e.g. coupling mode a and d, though a first Bragg pulse generating the state |Ψ 1 [yellow box]. (ii) a Raman laser couple mode a and b and entanglement dynamic is generated with the method of the Delta-kick Squeezing [46]. This method allow us to access various value of squeezing strength, τ E , through the control of atomic collision. The orientation of the quantum state, ν E , can then be manipulated through the phase of the laser transition creating the state |Ψ 2 [green box]. (iii) The MS beam-splitter can be engineered with a second Bragg pulse coupling mode b and c. The parameter ϕ MS is controlled by the phase of the laser, creating the entangled four-mode state |Ψ inp [blue box]. Bringing modes a and b and modes c and d to spatially overlap, with possible momentum transfer [49], the state |Ψ inp can be used for a differential measurement.

Conclusions
We have demonstrated a protocol for differential phase measurements with sensitivity enhanced by spin squeezing. A single two-mode spin-squeezed state is optimally manipulated via a mode-swapping operation that, although consisting essentially in a simple beam splitter, is enough to generate useful mode entanglement. Our mode-entangled and particle-entangled protocol overcomes the SQL and benefits from the cancellation of common-mode noise.
Our scheme represents a practical example of distributed quantum sensing [38,[83][84][85][86] where mode entanglement plays a crucial role. In particular, the corresponding particle-entangled mode-separable protocol using a single spin-squeezed state and a coherent state in independent interferometers cannot overcome a sensitivity scaling 1/N , with only a unit gain factor over the SQL. In contrast, our scheme achieves, with the same resources, a scaling 1/N 3/2 with oneaxis-twisted states and may even reach the Heisenberg scaling 1/N 2 with highly squeezed states.
To conclude, we notice that, after the completion of this work, a related sub-SQL differential scheme appeared [87]. There, an atomic distributed quantum sensing protocol is demonstrated experimentally by squeezing the total spin of the two interferometers, e.g.Ĵ a,b x +Ĵ c,d x , see Eq. (5), via atom-light quantum non-demolition measurements in a cavity. In our proposal, we squeeze only one collective spin, e.g.Ĵ a,b x , and then transfer the squeezing to the sumĴ a,b x +Ĵ c,d x via mode swapping. [ , valid for independent interferometers, and taking coherent spin states of N A and N B particles, respectively, such that ∆ 2 θ A,B = 1/N A,B , independently from the value of θ A,B . Finally, the optimal separable configuration is obtained for [42] L. Pezzè [63] The following relations hold between the coefficients θ MS , ϕ MS of Eq. (3) and |u bb |, |u cb |, δ cb in Eq. (9): |u bb | = cos θ MS , |u cb | = sin θ MS , δ cb = ϕ MS − π/2.
[64] We take a entangled state of N A particles and an coherent spin state of N B = N − N A particles in the interferometers A and B, respectively. For the mode-separable case we have The optimization of ∆ 2 (θ A − θ B ) with respect to N A , gives ∆ 2 (θ A − θ B ) ∼ 1/N . Instead, if two interferometers have the same number of particles, [65] M. Schulte, C. Lisdat, P. O. Schmidt, U.

A.1 Interferometer operations
Using the Stirling formula N ! ∼ √ 2πN N N e −N , valid for N 1, and absorbing the factor √ 2πN in the normalization of the state, the state Eq. (1) writes The delta function δ na+n b ,N expresses the fixed total number of particles. Neglecting this term allows for coherences between total number of particles which are, in turns, eliminated when measuring the total number of particles in output. Neglecting the delta function, the state |ψ 1 writes as where α a = α d = N/2. As a generalization, the single-mode coherent states |C(α a ) and |C(α d ) with α a = |α a |e iφ a,0 and α d = |α d |e iφ d,0 , respectively, will be used. For |α a | = |α d | = N/2 and φ a,0 = φ d,0 = 0, we recover the discussion in the main text.
The second step of the protocol is given by the one-axis twisting evolution exp{−iτ E (Ĵ {a,b} x ) 2 } followed by the rotation exp{−iν EĴ {a,b} z }. Within the Holstein-Primakoff transformation, we haveâ ∼ e iφ a,0 N/2, and thus is a single-mode squeezed-vacuum state, with squeezing parameter ξ 0 = e i(π/2−φ a,0 ) N τ E /4. In the derivation of the moment matrix, we will use ξ 0 = re iϕ 0 as a general squeezing parameter. The case discussed in the main text corresponds to the specific choice r = N τ E /4 and ϕ 0 = π/2 − φ a,0 = π/2. After squeezing, the rotation exp{−iν EĴ {a,b} z } = exp{−i(ν E /2)â †â } exp{i(ν E /2)b †b } is applied to modes a and b, thus affecting the phases φ a,0 and ϕ 0 that become φ a = φ a,0 − ν E /2 and ϕ = ϕ 0 + ν E , respectively. At this stage, we have A mode-swapping operation is then applied between modes b and c. Here, we assume a general unitary mode-swapping transformation: with δ cc = δ bb + δ bc − δ cb . The more specific case of Eq. (9) can be recovered by setting δ bb = δ cc = 0, which implies δ bc = δ cb .

A.2 Phase sensitivity equation
The uncertainty in the estimation of a generic linear combination v · θ = ν A θ A + ν B θ B using the method of moments [56, 60] is As a last assumption, we take N to be sufficiently large to apply the Holstein-Primakoff transformation: as explained in the main text, this amounts to takingâ ∼ e iφa N/2 andd ∼ e iφ d N/2. Being in the regime of large N will also allow us to neglect the two expectation values ψ inp |b †b |ψ inp and ψ inp |ĉ †ĉ |ψ inp found in the denominators. These two terms represent the average number of particles in modes b and c, respectively, after mode-swapping, so that each of them is smaller thann s and can be safely neglected if N/2 n s . Overall, we find We and thus recover Eq. (11). Next, we evaluate ∆ 2 (x b +x c ) on the mode-swapped squeezed stateÛ MS |S(ξ) b |0 c : and minimize the uncertainty using the MS transformationÛ MS . We calculate 2 are generalized quadrature operators. Note the similarity of this equation with Eq. (38) (next section), for λ = 0 and in the limit of large N , which is consistent with our previous discussion. Applying no mode-swapping to the state Even if the state |S(ξ) b was optimized through a rotation in thex b -p b quadrature plane, setting χ A = φ a − ϕ/2 = 0, that would correspond, for λ = 0, to an estimation uncertainty ∆ 2 (θ A −θ B ) θ A =θ B =π/2 ≈ 2e −2r /N + 2/N , implying a sensitivity gain G 2 = 2, at best, in the limit r → +∞. Hence, the high sensitivity gain typical of one-parameter estimation cannot be achieved here by a simple rotation in thex b -p b plane. To that aim, the effect of that rotation needs to be combined with a mode-swapping operation. To understand how modeswapping manages to enhance the sensitivity of the differential measurement, consider the chain of equalities: As shown in the next section, optimization of both the state |S(ξ) b and the mode-swapping parameters, in the large N regime, requires χ A = 0, χ B = 0 and |u bb | = |u cb | = 1/ √ 2. Application of these optimal conditions leads to Eq. (36), in which the differential quadrature operatorx b (λ) +x c (λ) already found above is recovered. This proves that, if evaluated on the mode-swapped stateÛ opt MS |S(ξ opt ) b |0 c , the variance of x b (λ) +x c (λ) reduces to that ofx b (λ − δ bb ) on |S(ξ opt ) b . Thus, equivalent results to those obtained in oneparameter estimation can be reproduced in a differential measurement using a single squeezed, mode-swapped state. More specifically, one gets For λ = 0 and using Eq.
with a sensitivity gain now given by G 2 = e 2r . Of course, the gain is not expected to increase indefinitely with squeezing and, as shown in the next section, a more accurate calculation predicts the existence of a maximum achievable gain:

B.1 Simultaneous optimization of theĴ z rotation and the mode-swapping parameters
Our goal here is to determine the rotation angle ν E in exp{−iν EĴ {a,b} z } and the mode-swapping parameters, Eq. (22), that minimize the uncertainty of our estimation method, Eq. (24). In particular, we consider the case v A = 1, v B = −1 (differential measurement case) and |α a | 2 = |α d | 2 = N/2. As a first step, we minimize that function with respect to the phases χ A and χ B defined in Eqs. (26) and (27), respectively. Given that Q(cot θ A , cot θ B ) is independent of these two parameters, we just focus on where we have usedn s to represent the average number of particles in mode b after squeezing and before the mode-swapping stage:n s = sinh 2 r = sinh 2 (N τ E /4). We notice that, by choosing sin χ A = sin χ B = 0 (which implies cos χ A = ±1 and cos χ B = ±1), we minimize the first squared term in Eq. (38) and simultaneously maximize the second one by a suitable choice of the sign of the two cosines. This yields the smallest possible value of the uncertainty, or the absolute minimum of Eq. (38) with respect to χ A and χ B . In order to maximize the second squared term in that equation, the sum of two terms of equal sign is clearly needed between round brackets. In the regime N n s , this requires cos χ A = cos χ B = ±1, which, by virtue of Eq. (26) and (27), translates to Moreover, since φ a,0 = φ d,0 = 0, ϕ 0 = π/2 and δ bb = 0 in the actual interferometer scheme, we end up with which gives k, m = 0, ±1, . . . , or δ cb = − π 8 + l π 2 , The second parametrization of the solutions clearly reveals the 4π periodicity of the sensitivity as function of ν E ; this is consequence of two simple facts: the 2π periodicity of the sensitivity as function of χ A and χ B , Eq. (38), and the fact that a 4π rotation in ν E is needed to get a full 2π variation in χ B , Eq. (27). By varying l and m, we get all the pairs (ν E , δ cb ) that minimize the function in Eq. (38) and, hence, the total uncertainty of the differential measurement. Making use of the optimal values of the parameters, from Eq. (38) one gets As a function of the amplitudes |u bb | and |u cb |, this is clearly minimized when the sum |u bb | + |u cb | attains its maximum value under the constraint |u bb | 2 + |u cb | 2 = 1, which corresponds to |u bb | = |u cb | = 1/ √ 2. Thus, the arbitrary choice of a balanced mode-swapping operation for our scheme is ultimately justified at least in the range of validity of the present analytical treatment. In the balanced configuration we obtain This is a function ofn s with a minimum at the pointn s = √ N /2, given by Thus, the gain is found to have a peak corresponding to the value G 2 max = √ N .

B.2 Optimization of theĴ z rotation for given values of the mode-swapping parameters
In this section, we consider Eq. (38) and assume φ a,0 = φ d,0 = 0, ϕ 0 = π/2, δ bb = 0 and |u bb | = |u cb | = 1/ √ 2 right from the start, this time, getting: where we have replaced χ A and χ B with their explicit expressions and set ν E ≡ ν, δ cb ≡ δ for the sake of a lighter notation. We aim to study the absolute minimum of the function with respect to the variable ν, for any fixed value of δ and of r. More specifically, we aim to determine the function ν min (δ, r) giving the minimum point ν min for any desired value of δ and r. Notice that, due to the 4π periodicity of Eq. (46) as a function of ν, infinitely many values can actually be attributed to ν min (δ, r) for any fixed δ and r. Such degeneracy can be lifted, nonetheless, by arbitrarily restricting the range of ν min (δ, r) to any given interval of length 4π. We will do it below. Let us first highlight two important symmetry properties of f (ν, δ, r).  62) and (66), respectively (all plotted as continuous black lines). Thus, we can provide an analytical justification to the linear relation observed in Fig. 4 between the MS parameter ϕ MS and the corresponding optimal rotation angle ν E .
have to be attributed to the first approximation A1. Let us develop this point further. To start with, we consider Eq. (45) completed with the addition of its phase-dependent part Q(cot θ A , cot θ B ): +4 N +n s (N −n s ) 2 + 2 N +n s (N −n s ) 2 cot 2 θ A + cot 2 θ B + (2n s + 1)n s (N −n s ) 2 (cot θ A − cot θ B ) 2 , The uncertainty obtained for a vanishing squeezing time τ E = 0 (corresponding to both r = 0 andn s = 0) is By comparison, the outcome of the numerical method, which can be exactly predicted for the case of no squeezing, is found to be: With no squeezing applied, the analytical approach generally overestimates the uncertainty of our scheme, as This discrepancy is caused by the second term in Eq. (70), namely (cot θ A − cot θ B ) 2 /N . For θ A = θ B one has (cot θ A − cot θ B ) 2 = 0, the extra term is minimum (vanishing) and ∆ 2 (θ A − θ B )| num = ∆ 2 (θ A − θ B )| th . For θ B = π − θ A , the difference ∆ 2 (θ A − θ B )| th − (cot θ A − cot θ B ) 2 /N turns out to be independent of θ A and θ B , giving a constant estimation uncertainty ∆ 2 (θ A − θ B )| num = 4/N that attains its minimum, at the shot-noise limit. All these details are visualized in Fig. 7.
The gain functions G 2 (θ A , θ B ) = 4∆ 2 (θ A − θ B )/N corresponding to Eqs. (69) and (70) have been plotted in panel (a) and (b), respectively. Looking at the figure, it becomes clear that the main effect of the term (cot θ A − cot θ B ) 2 /N is to break the left-right and up-down symmetries of Eq. (69), panel (a), and to enhance the sensitivity gain along the direction θ B = π − θ A . Finally, it is worth mentioning what the effect of a nonvanishing squeezing is for each method. As can be seen from Eq. (??), a non-vanishing squeezing time τ E 0 (implyingn s 0 in the sensitivity equation) is again related to the emergence of a term (cot θ A − cot θ B ) 2 which is, this time, not subtracted from but added to the sensitivity. As a consequence, the gain is now mainly enhanced along the θ B = θ A direction. In this case, the general behavior of the gain function resulting from the numerical method is similar to that obtained from the analytical prediction [compare panel (c) and (d) in Fig.  7].

D Two-Axis Twisting and Quantum Non-Demolition measurements
In Fig 8, we plot the sensitivity gain, G 2 , at the optimal working point and the bandwidth, R, as a function of the squeezing time τ E /τ ref in the case where the squeezing dynamic is generated through Two-Axis Twisting (TAT) dynamic [Fig 8a,b] and Quantum Non-Demolition (QND) measurements [Fig 8c,d]. Following Ref. [55], the state in Eq. 2 is now replace by and by a Gaussian state [75] in the case where entanglement is generated by QND protocol. The results are shown in the case where the number of particles is fixed to N = 100. The two protocols reproduce a similar behavior than the one shown in Fig.2 with One-Axis-Twisting dynamic. Because the orientation of the state is always optimized, QND measurements is characterized by two maximum. Here we conclude that the mode-swapping method can be extended to various squeezing procedures and is not limited to One-Axis-Twisting dynamic.   Figure 8: Maximum sensitivity gain, G 2 (π/2, π/2), and bandwidth, R, as a function of the squeezing strength τ E for different entanglement strategies. In each panels, the thin blue solid line are numerical results obtained for total N = 100 particles. The blue dashed lines refer to the separable differential scheme using a coherent spin state (τ B S = 0) a spin squeezed state (τ A S = τ E ) and no MS.