Phase Estimation from Atom Position Measurements

We study the measurement of the position of atoms as a means to estimate the relative phase between two Bose-Einstein condensates. First, we consider $N$ atoms released from a double-well trap, forming an interference pattern, and show that a simple least-squares fit to the density gives a shot-noise limited sensitivity. The shot-noise limit can instead be overcome by using correlation functions of order $\sqrt{N}$ or larger. The measurement of the $N\mathrm{th}$-order correlation function allows to estimate the relative phase at the Heisenberg limit. Phase estimation through the measurement of the center-of-mass of the interference pattern can also provide sub-shot-noise sensitivity. Finally, we study the effect of the overlap between the two clouds on the phase estimation, when Mach-Zehnder interferometry is performed in a double-well.


Introduction
Interferometry aims at the estimation of the relative phase between two wave-packets. In a standard optical interferometer, like the well known Mach-Zehnder setup [1], these two wave-packets correspond to the light traveling inside the two arms of the device, and the relative phase θ is acquired, for instance, as a result of different optical path length. After the phase is accumulated, the two wave-packets are recombined through a beam-splitter, and the signal at the two output ports depends on θ. The phase can be estimated by measuring the difference in intensities between these ports. Apart from photons, atoms can be employed for interferometric purposes as well [2]. The atoms present some interesting advantages with respect to light, especially due to the non-zero mass. In particular, the creation of atomic Bose-Einstein condensates (BECs) opened a new chapter in the field of interferometry. The BEC, which behaves like a macroscopic matter-wave, constitutes a coherent and well-controllable source of particles. This makes the BEC a promising system to measure the electromagnetic [3,4,5] or gravitational [6,7,8] forces. Moreover, the inter-atomic interactions in the BEC are a source of nonlinearity, which can be used to create non-classical states [9,10,11,12] useful to overcome the limit imposed by the classical physics on measurement precision [13,14].
A BEC interferometer can be implemented using a double-well trap [15,16,17,18,19,20,21,22,23], where the two wave-packets are localized about the two minima of the external potential. In such configuration, a relative phase θ can be accumulated by letting the system evolve in time in presence of an energy difference between the two potential minima, and in absence of coupling between the two wells. After this stage, one can, for example, recombine the wave-packets by implementing a beam splitter (thereby realizing a Mach-Zehnder interferometer). This will imply a further dynamical evolution during which atoms oscillate between the wells for a time which must be precisely under control, and over which interactions are negligible.
The recombination of the wave-packets can be done in a simpler way, just by releasing them form the double-well trap, so they form an interference pattern, as shown in Fig.1. In this manuscript we discuss how the information about the phase can be extracted from this pattern and derive the sensitivity for different estimation strategies. The manuscript is organized as follows. In Section 2 we formulate the problem and introduce the basic tool -the N-body correlation function, where N is the total number of atoms. In Section 3 we demonstrate that by performing a leastsquares fit to the measured density [15], the estimation sensitivity ∆θ is bounded by the shot-noise. As discussed in detail in Section 4, in order to overcome this limit, high-order spatial correlation functions must be measured, namely, of order not smaller than √ N . In particular, when estimation is performed using the N-th order correlation function, the sensitivity saturates the bound set by the Quantum Fisher Information (QFI) [24]. Then, in Section 5 we analyze an estimation scheme based on the detection of the position of the center-of-mass of the interference pattern, which can still yield sub-shot-noise sensitivity. Finally, in Section 6 we study the sensitivity of the Mach-Zehnder interferometer implemented in a double-well, and we observe that a non-zero overlap between the wave-packets dramatically reduces the sensitivity. Some details of the calculations are presented in the Appendix. The present manuscript is an extension of our previous work [25].

The model
To begin the discussion of different estimation methods based on position measurement, we introduce the two-mode field operator of a bosonic gas in a double-well potential, whereâ † /b † creates an atom in the left/right well. With the atoms trapped, the relative phase θ is imprinted between the modes. This stage is represented by a unitary evolution U (θ)=e −iθĴz of the initial state |ψ in of the double-well system. The three operatorŝ form a closed algebra of angular momentum. With the phase acquired, the trap is switched off and the two clouds described by the mode functions ψ a/b (x, t) freely expand.
The most general quantity, containing the statistical information about the positions of particles forming the interference pattern, is the conditional probability of finding N particles at positions x N = (x 1 . . . x N ). It can be expressed in terms of the N-th order correlation function p N ( Here, |ψ out denotes the state after the interferometric transformation, |ψ out = e −iθĴz |ψ in . To provide a compact and useful expression for this probability, we take following steps.
First, we decompose the initial state in the well-population basis, |ψ in = N n=0 C n |n, N−n and suppose that the expansion coefficients are real and posses the symmetry C n = C N −n . As we will argue later, such choice of C n 's is natural in context of this work. We switch from the Schrödinger to the Heisenberg representation, where the field operator evolves according to, The next setp is to introduce the basis of the coherent phase states [26] defined as (where |0 is the state with zero particles). The action of the field operator on these states can be written in a simple form, where u θ (x, ϕ; t) = ψ a (x, t)e i 2 (θ+ϕ) + ψ b (x, t)e − i 2 (θ+ϕ) . Next, we expand the Fock states in the basis of the coherent states Thus we can easily write the result of action of the field operator on the input state, Now evaluation of G N ( x N , θ) (and therefore the p N ( x N , θ) as well) is straightforward -we act N times withΨ θ on the input state and calculate the modulus square of the result. After normalization we obtain In the remaining part of the manuscript, we fix t large enough so that the interference pattern is formed. In this regime, the physical properties of the system change only by scaling ∼ √ t of the characterisitc dimensions of the system. The probability (2) is the starting point for the following discussion of various phase estimation strategies.

Estimation via the fit to the density
The simplest way of estimating the value of θ is by fitting the average density to the interference pattern, as the position of the maximum depends on the relative phase between the two clouds. Such fit is commonly employed with BECs in double-wells, in order to determine, for instance, the phase coherence in the system [6,9,15].
In the experimental realization, the interference pattern is sampled using M bins located at positions x i (i = 1 . . . M). The number of particles n i in each bin is measured m times, giving the set n i /m. If the size ∆x of each of M bins is small, n i is given by the average density The value of θ is determined from the least squares formula The fluctuations in each bin, ∆ 2 n i = lim m→∞ m k=1 (n (k) i − n i ) 2 , are calculated from the probability p(n i |θ) of detecting n i particles in the i-th bin (for details of derrivation, see Appendix A), Using Eq.(2) we obtain In Eq. (4), n i and ∆ 2 n i are assumed to be known, since in the phase estimation stage the only measured data aren i . The quantities n i and ∆ 2 n i are instead constructed during the calibration stage, preceding the phase estimation stage, by repeating the experiment with different known values of θ. If the number of experiments in the calibration is large, and in absence of thermal and technical noise, the measured n i and ∆ 2 n i will tend to the theoretical predictions given in Eq. (3) and (6), respectively. Our goal at this point is to determine how the quantum fluctuations ∆ 2 n i in the i-th bin influence the sensitivity of the phase estimation via the fit (4). To this end, we employ the concept of the maximum likelihood estimation (MLE) [27,28]. If some quantity ξ is measured, the MLE is defined as the choice of θ which maximizes the conditional probability P (ξ|θ) for the occurrence of ξ given θ. That is, the phase θ is estimated from the condition d dθ P (ξ|θ) = 0. In case of the fit discussed here, the estimation is based on the measured average occupationsn i . If the number of measurements m is large, then according to the central limit theorem the probability distribution for the averagen i in the i-th bin tends to the Gaussian p(n i |θ) 2∆ 2 n i /m . In every shot the atom counts are correlated between the bins. However, in order to link the MLE with the sensitivity of the least squares fit, we construct the likelihood function as if the measurement resultsn i andn j , with i = j, were uncorrelated. Thus the total probability of measuring the values {n} = (n 1 . . .n M ) is a product P ({n}|θ) = M i=1 p(n i |θ). We note that in this case, indeed, the condition for the MLE, d dθ P ({n}|θ) = 0 coincides with Eq.(4). The MLE sensitivity saturates the Cramer-Rao Lower Bound [27,28], ∆ 2 θ = F −1 . Here F is the Fisher information (FI), Therefore, the sensitivity for the least squares fit is given by Eq.(7) as well.
In the following, we demonstrate that this sensitivity is bounded by the shot-noise. Let us assume for the moment that the second term in the Eq.(6) -which is proportional to (∆x) 2 -can be neglected. Then, as can be seen from Eq.(6), the particle number distribution is Poissonian. The FI from (7) reads with the one-particle probability (9) We now calculate the Fisher information (8) explicitly. As the interference pattern is formed after long expansion time, the mode functions can be written as whereσ = t m ,ψ is a Fourier transform of the initial wave-packets, common for ψ a and ψ b and the separation of the wells is 2x 0 . This gives Notice that when the expansion time is long, the functionψ varies slowly, as compared to the period of oscillations of sin ϕ and cos ϕ. Therefore, in the above expression, one can substitute the oscillatory part with its average value. Using the normalization ofψ we obtain As 0 ≤ a ≤ 1, we have 0 ≤ F ≤ mN. Therefore the Fisher information for the fit is always smaller than the shot-noise, giving ∆θ ≥ ∆θ SN = 1 √ m 1 √ N for any two-mode input state (∆θ SN denotes the shot-noise sensitivity). Below we argue that the inclusion of the second term in the fluctuations in Eq. (6) does not improve the sensitivity.
In Eq.(6), the first term G 1 (x i , θ)∆x scales linearly with N while the second term, as a function of N, is a polynomial of the order not higher than two, aN 2 + bN + c. The fluctuations ∆ 2 n i must be positive, thus a ≥ 0. Otherwise, for large N, no matter how small ∆x, we would have ∆ 2 n i < 0. If a > 0, the positive second term enlarges the fluctuations and thus worsens the sensitivity. When a = 0, the first and second terms in Eq.(6) scale linearly with N, and for small ∆x the second term can be neglected, thus we end up again with Eq.(8) ‡.
In order to calculate the sensitivity ∆θ = 1 √ F in Eq.(11), we consider the ground states of the BEC in a double well potential. In the two-mode approximation, the Hamiltonian of the system readŝ We construct a family of states A by finding ground states of the above Hamiltonian for various values of the ratio γ = E C N E J . And so, for γ > 0, the elements of A are number-squeezed states and tend to the twin-Fock state |ψ in = | N 2 , N 2 with γ → ∞. For γ < 0 the elements of A are phase-squeezed states [23]. With γ → −∞, the ground state of (12) tends to the NOON state |ψ in = 1 Notice that for all |ψ in ∈ A, the coefficients C n , which were introduced in previous section, are real and symmetric. For each state in A, we calculate a = 2 N Ĵ x , and insert it into Eq. (11). The sensitivity shown in Fig.2 is clearly bounded by the shot noise.
This limitation for the sensitivity can be explained as follows. The value of the FI given by Eq. (8) is expressed in terms of the single-particle probability. We expect the useful non-classical many body correlations to decrease the value of ∆θ, but the FI (8) is insensitive to these correlations, and thus must be bounded by the shot-noise. In the next section we demonstrate that the estimation based on the measurement of position correlations can improve the phase sensitivity.

Estimation via the correlation functions
In the estimation protocol discussed in this section, the phase θ is deduced from the measurement of the k-th order correlation function G k ( x k |θ). As previously, we choose to deduce θ using the MLE: a set of k positions x k is measured, and the phase is chosen We notice that by setting k = 1, i.e. the estimator is a single-particle density, we recover the FI from Eq. (8). Therefore, the measurement of positions of N particles used as independent is, in terms of sensitivity, equivalent to fitting the average density to the interference pattern, and is limited by the shot-noise.
Let us calculate the FI for the case k = N, corresponding to the measurement of the full N-body correlation function. We represent the mode functions ψ a/b using Eq.(10). This expression allows to calculate u θ (x, ϕ), and, in turn, the probability (2), which is then inserted into Eq. (13). The integrals over space can be performed analytically (see Appendix B for details) and the outcome is Here, by F Q we denote the QFI, which is a maximal value of the Fisher information with respect to all possible measurements [24]. In the case of pure states, the value of the QFI is given by 4m times the variance of the phase-shift generator, thus in our case it reads As denoted by open circles in Fig.3, the Eq. (14) gives ∆θ = √ m∆θ SN for the coherent state (γ = 0), and overcomes this bound for all |ψ in ∈ A with γ < 0. The NOON state gives the Heisenberg limit, ∆θ HL = 1 √ m 1 N . We now discuss how the sensitivity given by the inverse of Eq.(13) changes for k < N. The space integrals for k = 1, N cannot be evaluated analytically, thus we calculate the FI numerically taking Gaussian wave-packets with the initial width σ 0 = 0.1 and half of the well separation x 0 = 1. The Fig.3 shows how the sensitivity from Eq. (13) for N = 8 atoms changes with increasing k as a function of |ψ in ∈ A. The sensitivity improves with growing k, and goes below the shot-noise limit at k min = 4.  For higher numbers of particles, we numerically checked that k min tends to √ N . Therefore, one would have to measure the correlation function of the order of at least √ N in order to beat the shot-noise limit. In a realistic experiment with cold atoms, where N ≃ 1000, it would be very difficult to use the correlation function of such high order for phase estimation. The biggest difficulty resides in the calibration stage, during which one would need to experimentally probe a function of a k dimensional variable x k (with k > √ N ) for different values of theta. In the following section we present a phase estimation scheme based on the measurement of the center-of-mass of the interference pattern. Although the probability for measuring the center-of-mass at position x is a function of just a one-dimensional variable, it can still provide the sub-shot-noise sensitivity. Nevertheless, we will demonstrate that the implementation of this estimation protocol can be challenging.

Measurement of all N atoms
In order to estimate θ from the measurement of the center-of-mass, one has to go through a relatively simple calibration stage. Positions of N atoms are recorded independently and from this data location of the center-of-mass is deduced. Many repetitions of the experiment give the function p cm (x|θ) of a one-dimensional variable. The expression for this function can be extracted from the full N-body probability (2) by  where "δ" is the Dirac delta. To provide an analytical expression for this probability, we model the mode-functions by Gaussians as in Eq. (15). Using a reasonable assumption that the initial separation of the wave-packets is much larger than their width, i.e. e −x 2 0 /σ 2 0 ≪ 1, we obtain The details of this derivation are presented in Appendix C. Notice an interesting property -the above probability depends on θ only for states with non-negligible NOON components C 0 and C N , as already noticed in [29]. When p cm (x|θ) is known, the phase can be estimated using the MLE, as used in the previous sections. Then once again the sensitivity is given by the inverse of the FI, which can be calculated analytically, where m is the number of experiments. In Fig.4 we plot the sensitivity calculated by the inverse of the FI (17) as a function of |ψ in ∈ A with γ ≤ 0. Although the estimation through the center-of-mass is not optimal (∆θ > 1 √ F Q ), the sensitivity can be better than the shot-noise, with ∆θ → ∆θ HL for |ψ in → NOON. The calibration stage is not as difficult as in the case of high-order correlations, however phase estimation based on the center-of-mass measurement demands detection of all N atoms [30,31,32,33], as we show below.

Measurement of k < N atoms
If the measurement of the center-of-mass is based on detection of k < N atoms, the probability (16) transforms into where p k ( x k |θ) = d x N −k p N ( x N |θ). The probability (18) can be calculated in a manner similar to that presented in Appendix C. The result is where a = 2 Notice that for k = N we recover the result from the previous section a = 2C 0 C N = 1 2 (C 0 +C N ) 2 (as we are using the symmetric states, C 0 = C N ). The FI for the probability (19) can be calculated analytically, Let us now evaluate a -and thus F -for various k ≃ N. For k = N and the NOON state, we have C 0 = C N = 1 √ 2 , giving a = 1 and F = mN 2 . From Eq. (20) we notice that, for any k, a is the sum of N − k terms, each depending on the coefficients C i and C i+k . And so, for k = N − 1, a will be maximal for a NOON-like state with C 0 = C N −1 = 1 2 and C 1 = C N = 1 2 . For this state we obtain a = 1 √ N , and for large N the value of the FI is F = mN. Therefore, the phase estimation using the center-of-mass of N − 1 particles gives a sensitivity bounded by the shot-noise. Each loss of an atom decreases the FI roughly by a factor of N, drastically deteriorating the sensitivity. In Fig.5 we plot the sensitivity √ m∆θ calculated with the FI from Eq.

Formulation
So far, we focused on the position measurement of atoms released from a double-well trap, and studied the phase estimation sensitivity. The fit to the density gives sensitivity limited by the shot-noise, and this bound can be overcome by phase estimation with correlation functions of the order of at least √ N . As it is difficult to measure these correlations in the experiment, it will be challenging to beat the shot-noise limit using the interference pattern. Although the sensitivity of the phase estimation based on the center-of-mass measurement can also be sub-shot-noise, the protocol is extremely vulnerable to the loss of particles.
In the above scenario, the sub-shot-noise sensitivity, which relies on non-classical particle correlations, is reached by directly measuring spatial correlations between the atoms forming the interference pattern and using the latter as estimators for the phase shift. On the other hand, it is well known that the Mach-Zehdner Interferometer (MZI) can easily provide sub-shot-noise sensitivity just by a simple measurement of the population imbalance between the two arms and a proper choice of the input state |ψ in . This is because, in the MZI, the correlations between the two modes carry the part of the information contained in the particle correlations which is useful for phase estimation. When the clouds are released from the trap and the two modes start to overlap, the correlations between the two modes are lost, since an atom detected in the overlap region cannot be told to have come from either of the two initially separated clouds. This is the reason why it is necessary then to measure directly high-order spatial correlation functions in order to reach sub-shot-noise sensitivity.
It would be thus interesting to quantify the effect of the wave-packets' overlap on the sensitivity of the MZI. This analysis has also a practical interest since, in the implementation of the atomic MZI, the precision of the population imbalance measurement can be improved by opening the trap and letting the clouds expand for a while. In this way the density of the clouds drops, facilitating the measurement of the number of particles. However, during the expansion, the clouds inevitably start to overlap, leading to loss of information about the origin of the particles, as noted above.
In this section, we show how the increasing overlap deteriorates the sensitivity of the MZI in two estimation scenarios.
The MZI consists of three stages: two beam-splitters represented by unitary evolution operators e ∓i π 2Ĵ x separated by the phase acquisition e −iθĴz . The atomic MZI can be realized as follows. Consider a two-mode system governed by the Hamiltonian (12) with E C = 0. The first beam-splitter is done by letting the atoms tunnel between the two wells for t = π 2 E J so the unitary evolution operator readsÛ 1 = e −i π 2Ĵ x . Then, an inter-well barrier is raised, in order to supress the oscillations (E J = 0) and a phase between the wells is imprinted, givingÛ 2 = e −iθĴz . The interferometric sequence is closed by another beam-splitter,Û 3 = e i π 2Ĵ x . The full evolution operator readŝ where we used commutation relations of the angular momentum operators. In order to analyze the sensitivity of the MZI, we introduce the conditional probability p N ( x N |θ) of detecting N atoms at positions x N = x 1 . . . x N . To evaluate this probability for any initial state of the double-well system |ψ in , we take the same steps as in Section II. In the Heisenberg representation, the field operator evolves aŝ Again, we express the action of the field operator on |ψ in using the basis of the coherent phase-states and obtain Eq.(1) with Therefore, the probability p N ( x N |θ) for the MZI is given by Eq.(2) with the u θ (x, ϕ; t) function defined above.

Measurement of the population imbalance
The most common phase-estimation protocol discussed in context of the MZI is the measurement of the population imbalance between the two arms of the interferometer.
In order to assess how the sensitivity of this protocol is influenced by the expansion of the wave-packets, we introduce the probability of measuring n L atoms in the left sub-space as follows This probability depends on the expansion time via ψ a,b (x, t), which enter the definition of p N ( x N |θ). The sensitivity for various expansion times, if m ≫ 1 measurements are performed, can be calculated using the error propagation formula, where is the average value of the population imbalance and are the associated fluctuations. The probability (23) resembles Eq.(5), and the above moments are calculated as in Appendix A resulting in The two lowest correlation functions for the MZI read G 1 (x|θ) = Ψ † θ (x, t)Ψ θ (x, t) and t) , with the field operator from Eq. (22), and the averages calculated with the input state. When the two wave-packets don't overlap, i.e. ψ a (x, t)ψ * b (x, t) ≃ 0 for all x ∈ R, Eq.(24) simplifies to This is the well-known expression for the sensitivity of the population imbalance between separated arms. It gives ∆θ ≤ ∆θ SN for all |ψ in ∈ A with γ ≥ 0. We investigate the impact of the overlap on the sensitivity (24) by modelling the free expansion of the wave-packets ψ a/b (x, t) by Gaussians, and take x 0 = 1 and the initial width σ 0 = 0.1. In Fig. 6 we plot the sensitivity √ m∆θ taking θ = 0 and N = 100 for three different expansion times τ . The initial sensitivity deteriorates as soon as the condensates start to overlap, and the sub-shotnoise sensitivity is lost for long expansion times. We attribute this decline to the loss of information about the correlations between the modes. Therefore, special attention has to be paid to avoid the overlap when letting the two trapped condensates spread. Although we expect that the expansion facilitates the atom-number measurement, the conclusion of this section is that any overlap of the spatial modes has a strong negative impact on the sensitivity of the MZI.

Estimation via the center-of-mass measurement for the MZI
In Section 5.1, we have demonstrated that when the two wave-packets overlap and form an interference pattern, phase estimation based on the center-of-mass measurement can give sub-shot-noise sensitivity. Here we study the same estimation strategy applied in the MZI case. Again, we start with the probability p cm (x|θ) of measuring the center-ofmass at position x, Using p N ( x N |θ) for the MZI, one can analytically calculate the center-of-mass probability only in the limit of small θ, where With this probability, we can again calculate the sensitivity using the error propagation formula [27,28], These two moments can be easily calculated with Eq.(26), giving, in the limit θ → 0, Notice that when the initial size of the Gaussians tends to zero, we recover the sensitivity from Eq.(25) (in the limit of θ → 0). This is not surprising, as when the mode-function are point-like, the measurements of the center-of-mass and the measurement of the population imbalance are equivalent, and related by x cm = 2x 0 n N . Therefore, for small σ, the measurement of the center-of-mass yields sub-shot-noise sensitivity for all |ψ in ∈ A with γ > 0. However, for non-zero σ, the second term in Eq.(27) spoils the sensitivity. This is because N 4 Ĵ x 2 ≥ 1 N is always satisfied. Even if the first term scales better than at the shot-noise limit, the other one does not, and will dominate for large N.
¿From what we presented in this Section, we conclude that both the population imbalance and the center-of-mass measurements can give sub-shot-noise sensitivity for the MZI, but both are very sensitive to the growing size of the wave-packets.

Conclusions
In this manuscript we have discussed in detail how the measurement of positions of atoms forming an interference pattern can be useful in context of atom interferometry. We showed that the phase estimation based on the fit to the density gives sensitivity limited by the shot-noise, because the FI is expressed in terms of the single-particle probability only. The sensitivity can be improved below the shot-noise limit by estimating the phase using correlation functions of order at least √ N. Moreover, we demonstrated that the information contained in the N-th order correlation function allows to perform an optimal detection strategy, reaching Heisenberg-limited sensitivity when NOON states are used. We also showed that the measurement of the position of the center-of-mass of the interference pattern gives sub-shot-noise sensitivity for all states with non-negligible NOON component. Both the measurement of high-order correlations and the centerof-mass position are difficult to perform. The former requires the construction of a function of highly-dimensional variables, while the latter works well only if all N atoms forming the interference pattern are detected. We attribute the difficulty to obtain the sub-shot-noise sensitivity to the fact that, after formation of the interference pattern, the modes cannot be distinguished, and the information useful for interferometry is only contained in the correlations between the particles. These correlations are very difficult to extract from the experimental data, therefore reaching sub-shot-noise sensitivity with two interfering BECs might prove very challenging. In the final part of this work, we turned our attention to the MZI, which is known to provide sub-shot-noise sensitivity for the simpler measurement of the population imbalance between the two clouds. We have shown that the sensitivity of the MZI is strongly influenced by a non-zero overlap between the two wave-packets, both in case of the population imbalance and the centerof-mass measurement.

Appendix A. Bin fluctuations
To derive the expression for the average number and the fluctuations of the atom count in a bin, we use the probability from Eq.(5). With help of Eq.(2) we obtain Since the distribution of n i inside the integrals is binomial, one can easily calculate the average value as follows Notice, that using the definition of a θ and b θ , this can be rewritten as t) is a one-particle density. In the same manner, one can demonstrate that the fluctuations of the number of particles in the i-th bin are given by When the bin size is much smaller than the characteristic variation of the mode functions, the above integrals are approximated by Appendix B. The FI for k = N In this Appendix we prove the Eq. (14). We start from the expression for the Nparticle probability (2) for the wave-packets after long expansion time, as in Eq. (10). A straightforward evaluation shows that in this case, the probability can be written as where f is real and reads x ĩ σ 2 cos Moreover, we assume that the initial separation of the wells of the trapping potential is much larger than the width of the mode functions. Under this assumption, the probability (B.1) is normalized. This probability is inserted into the definition of the Fisher information. As f is real, we immediately obtain Now, the order of the integration can be reversed, and the space integrals are performed first. As the mode-functions are normalized, the result is following The phase integral can be now easily evaluated, giving which is the Eq. (14).

Appendix C. Evaluation of the center-of-mass probability
In this Appendix we derive the expression for the probability of detecting the center-ofmass at position x, as in Eq. (16). The definition of p cm (x|θ) relates it to the full N-body probability by To calculate this probability we assume a long expansion time and use Eqs (B.1) and (B.2). Then, we notice that the Dirac delta can be represented as the Fourier transform, This equation can be rewritten in a following way is the Fourier transform of the probabilityp cm (x|θ). To provide an analytical expression for this probability, we assume that the initial wave-packets are Gaussian as in Eq. (15). The integration over space is performed giving p cm (k|θ) = with k 0 = 2N x 0 σ 2 and w = 2N σ 0 σ 2 . The function I(k, ϕ, ϕ ′ ) consists of three peaks, located at k = ±k 0 and k = 0. When the well separation 2x 0 is large compared to the initial width of the trapped wave-packets σ, then these three peaks are separated and do not overlap. Therefore, [I(k, ϕ, ϕ ′ )] N can be approximated by the sum of N-th powers of its three components