Overcoming loss of contrast in atom interferometry due to gravity gradients

Long-time atom interferometry is instrumental to various high-precision measurements of fundamental physical properties, including tests of the equivalence principle. Due to rotations and gravity gradients, the classical trajectories characterizing the motion of the wave packets for the two branches of the interferometer do not close in phase space, an effect which increases significantly with the interferometer time. The relative displacement between the interfering wave packets in such open interferometers leads to a fringe pattern in the density profile at each exit port and a loss of contrast in the oscillations of the integrated particle number as a function of the phase shift. Paying particular attention to gravity gradients, we present a simple mitigation strategy involving small changes in the timing of the laser pulses which is very easy to implement. A useful representation-free description of the state evolution in an atom interferometer is introduced and employed to analyze the loss of contrast and mitigation strategy in the general case. (As a by-product, a remarkably compact derivation of the phase-shift in a general light-pulse atom interferometer is provided.) Furthermore, exact results are obtained for (pure and mixed) Gaussian states which allow a simple interpretation in terms of the alignment of Wigner functions in phase-space. Analytical results are also obtained for expanding Bose-Einstein condensates within the time-dependent Thomas-Fermi approximation. Finally, a combined strategy for rotations and nonaligned gravity gradients is considered as well.

Long-time atom interferometry is instrumental to various high-precision measurements of fundamental physical properties, including tests of the equivalence principle. Due to rotations and gravity gradients, the classical trajectories characterizing the motion of the wave packets for the two branches of the interferometer do not close in phase space, an effect which increases significantly with the interferometer time. The relative displacement between the interfering wave packets in such open interferometers leads to a fringe pattern in the density profile at each exit port and a loss of contrast in the oscillations of the integrated particle number as a function of the phase shift. Paying particular attention to gravity gradients, we present a simple mitigation strategy involving small changes in the timing of the laser pulses which is very easy to implement. A useful representation-free description of the state evolution in an atom interferometer is introduced and employed to analyze the loss of contrast and mitigation strategy in the general case. (As a by-product, a remarkably compact derivation of the phase-shift in a general light-pulse atom interferometer is provided.) Furthermore, exact results are obtained for (pure and mixed) Gaussian states which allow a simple interpretation in terms of the alignment of Wigner functions in phase-space. Analytical results are also obtained for expanding Bose-Einstein condensates within the time-dependent Thomas-Fermi approximation. Finally, a combined strategy for rotations and nonaligned gravity gradients is considered as well.

I. INTRODUCTION
In this article we analyze in detail the loss of integrated contrast in open atom interferometers, where classical phase-space trajectories do not close, with particular emphasis on the effects of gravity gradients. A simple but effective mitigation strategy to overcome such loss of contrast is presented.
The potential of light-pulse atom interferometry [1,2] for high-precision measurements has been amply demonstrated with its successful implementation in extremely sensitive inertial sensors, including gyroscopes [3][4][5], gradiometers [6] and the currently most precise absolute gravimeters [7,8]. It has already found applications in accurate measurements of fundamental constants [9][10][11][12][13] and tests of fundamental properties [14,15], and it is a key ingredient in plans for future tests of the equivalence principle in space [16,17], next-generation satellite geodesy missions [18] or even alternative proposals for gravitational-wave detection [19]. Achieving higher sensitivities requires extended interferometer times (the phase shift generated by accelerations, for instance, grows with the square of the interferometer time). Remarkable progress in this direction has recently been made with atomic fountains on the ground [20]. Nevertheless, the natural environment for reaching even longer interferometer times are the microgravity conditions provided by space missions [16,21] and being preliminarily investigated in parabolic flights [22], drop towers [23,24] and sounding rockets [25].
Measurements based on atom interferometry like those mentioned above typically involve keeping track of the oscillations in the total number of atoms at each exit port as certain parameters, such as the half-interferometer time T , are changed. Besides other challenges of practical nature, the long interferometer times planned for future space missions face a more fundamental limitation due to gravity gradients. These cause a relative shift in position and momentum between the pair of wave packets contributing to the superposition at each exit port, which leads to a reduction of the integrated contrast (i.e. a smaller amplitude of the oscillations in the integrated atom number at each port) and a corresponding decrease in sensitivity. As an example, plans for the European Space Agency's STE-QUEST mission [16] have estimated a 40% loss of contrast due to gravity gradients despite a required effective temperature for the atom cloud below 70 pK, and this is the main driver behind the need for such a narrow momentum distribution (although it has other added benefits on diffraction efficiencies or systematics associated with the Rabi frequency).
Such a loss of contrast is, in fact, a generic feature of open interferometers, where the classical trajectories characterizing the motion of the wave packets do not close in phase space and the resulting relative displacements in position and momentum, δR and δP, cause a reduction of the quantum overlap between the interfering states in each exit port. In position representation this reduction is due to two effects: first, a decreasing spatial overlap of the two envelopes; second, spatial oscillations in the product of one wave function and the other's complex conjugate, whose integral over space gives the quantum overlap of the two states. As we will see in forthcoming sections, in situations relevant for atom interferometry (at least before implementing our proposed mitigation strategy) the second effect typically becomes important well before the spatial displacement becomes comparable to the size of the envelopes and the first effect starts to play a role.
The spatial oscillations described in the previous paragraph give rise to a fringe pattern in the spatial density profile of each exit port which is closely related to arXiv:1401.7699v1 [physics.atom-ph] 29 Jan 2014 the contrast loss. This suggests devising a simple mitigation strategy that consists in eliminating the fringe pattern by means of a small but properly chosen change δT in the timing of the last beam-splitter pulse (i.e. the time at which it is applied). Doing so hardly alters the the momentum displacement δP, but it can be used to change δR in such a way that its effect essentially cancels that due to δP and the fringe spacing becomes much larger than the size of the envelope. In this paper we will mainly focus on gravity gradients aligned with the direction of the laser beams, but will express our results using a vector notation that makes them directly applicable to the general case (when the principal axes of the gravity gradient tensor are no longer aligned with the lasers) as well. However, since nonaligned gravity gradients produce also displacements along the transverse directions, adjusting the timing of the last beam-splitter, which only generates displacements along the longitudinal direction, is not enough. As explained in Sec. V, in that case one can deal with the transverse displacements by using the same kind of scheme based on a tip-tilt mirror which has already been employed by several groups to compensate the effects of rotations [20,[26][27][28][29].
Even though in many of the calculations that will be presented in this article an exact treatment of the gravity gradients is possible (at least for time-independent ones), the results are simpler and more transparent when treating the gravity gradient tensor Γ perturbatively. This is, in fact, an excellent approximation for the calculation of the contrast because the parameter Γ T 2 controlling the perturbative expansion (where Γ denotes the largest eigenvalue in absolute value) is typically much smaller than unity: even for half-interferometer times of T = 10 s it is still of the order of 10 −4 . Therefore, one can perform a perturbative calculation of the classical trajectories that characterize the motion of the wave packets to determine the relative displacement between the corresponding pair of wave packets at each exit port, which is of order Γ T 2 (although somewhat enhanced by a factor v rec T ). Furthermore, when determining how the shape of the wave packets evolves in time (or, more precisely, how the centered wave packets defined in Appendix B evolve), it is sufficient to consider their free evolution (neglecting the gravity gradients). This is because in the absence of relative displacement the states of the two wave packets coincide up to a phase and have maximum quantum overlap (hence, independent of Γ). Therefore, since the relative displacement is itself of order Γ, the effect of the gravity gradients on the evolution of the centered wave packet would generate corrections to the quantum overlap of higher order in Γ.
Although the main goal of this article is to study the loss of integrated contrast in open interferometers and how to overcome it, it is worth pointing out that many of the techniques and results presented here can also be directly used to explain how gravity gradients, rotations and changes in pulse timing modify the density profile observed at the exit ports of an atom interferometer [30], as measured in the experiments reported in Refs. [20,24,29]. Moreover, as a side result we provide in Appendix B a remarkably compact derivation of the general formula for the phase shift in presence of time-dependent gravity gradients and uniform force fields.
The paper is organized as follows. The general framework, the reasons for the loss of contrast in open interferometers and the mitigation strategy are presented in Sec. II, which contains our central findings. A detailed analysis for pure Gaussian states, with exact results and quantitative examples, is provided in Sec. III. They illustrate the phase-space interpretation of the loss of contrast and the mitigation strategy in terms of the alignment of squeezed Wigner functions. In addition, in Sec. III B we discuss the new features connected with mixed states in general and their particularization to the Gaussian case. Interferometers based on expanding Bose-Einstein condensates (BECs) are investigated in Sec. IV, where we give analytical results using the scaling approach and the time-dependent Thomas-Fermi approximation as well as quantitative examples. The case of nonaligned gravity gradients is considered in Sec. V. Finally, our results and future perspectives are discussed in Sec. VI. A number of technical aspects have been included in the appendices. The classical phase-space trajectories (and the relative displacements) in the presence of gravity gradients, uniform force fields and taking into account the laser kicks are computed in Appendix A. A representation-free description for quadratic Hamiltonians of the state evolution in the atom interferometer is presented in Appendix B, including a very compact derivation of the general result for the phase shift. In Appendix C we obtain the late-time result for the free evolution of a wave packet. Finally, the Gross-Pitaevskii equation and its solution for expanding BECs within the framework of the scaling approach and the time-dependent Thomas-Fermi approximation are briefly reviewed in Appendix D and its connection with free evolution at late times is discussed in some detail.

II. LOSS OF CONTRAST AND MITIGATION STRATEGY
A. Representation-free description for atom interferometers Throughout the paper we will focus on light-pulse atom interferometers, based on the diffraction of matter waves by time-modulated laser pulses which induce transitions with a fixed momentum transfer (and possibly changing the internal state). By selecting their duration and intensity, one can have π and π/2 pulses acting respectively as mirrors and beam splitters. For simplicity we will assume idealized beam splitters and mirrors with 100% efficiency, no dispersive effects within the relevant momentum range for the atoms (which would otherwise change the shape of the wave packet) and sufficiently short pulse duration so that the evolution of the external degrees of freedom can be neglected during the pulses. Furthermore, we will assume a perfect read-out of the two exit ports thanks to good spatial separation or state labeling. Finally, when considering quantum degenerate gases, their density during the interferometer sequence should be low enough so that nonlinear interactions can be neglected and their dynamics can be accurately described in terms of the Schrödinger equation for noninteracting particles.
Besides the interaction with the laser fields during the short duration of the laser pulses, the dynamics of the atoms will be governed by the following quadratic Hamiltonian:Ĥ It can describe the effect of time-dependent gravity gradients and uniform force (or acceleration) fields, characterized respectively by the tensor Γ(t) and the vector g(t). It can also be employed to describe the effect of rotations from the point of view of the non-rotating frame: one needs then to take into account the rotation of the momentum transfers k i associated with the different laser pulses, and possibly of Γ(t) and g(t) too [31]. As shown in Appendix B, for quadratic Hamiltonians like that of Eq. (1) there is a useful representation-free description for the time evolution of the atoms within the interferometer in terms of the evolution of a centered wave packet plus classical trajectories which include the kicks from the laser pulses and characterize the motion of each wave packet contributing to the superposition. For interferometers consisting of an initial and a final beam splitter together with several π pulses in between (the Mach-Zehnder configuration being a typical example) the state of the system during the interferometer phase (between the first and last beam-splitters) is given by where |ψ c (t) is a centered wave packet evolving in a purely quadratic potential and we have omitted the time dependence of the phases and displacement vectors to avoid an unnecessarily cumbersome notation. The displacement operatorD(χ) is defined aŝ where we used a vector notation for phase-space quantities, so thatξ = (x,p) T , and introduced the symplectic form The displacement vectors χ(t) = R(t), P(t) T in Eq. (2) correspond to the phase-space vectors for the classical trajectories associated with each branch, a and b, of the interferometer. Provided that the two exit ports can be perfectly distinguished, due to absence of spatial overlap or different internal states, only two wave packets contribute to each port, as depicted in Fig. 1. Therefore, we have the following superposition at one of the exit ports (after the second beam-splitter): with δχ = χ 2 − χ 1 and δΦ = Φ 2 − Φ 1 + χ T 1 J χ 2 /2 , where the extra phase arises from the composition of displacement operators: The results is analogous for the state |ψ II (t) at the second exit port. While we will focus on pure states throughout this section, mixed sates can be treated as an incoherent ensemble of pure states, as explained in detail in Sec. III B.

B. Loss of contrast in open interferometers
We will employ the term "open interferometer" whenever δχ = 0 at the time of detection. This will lead in general to a loss of contrast when considering the integrated number of particles at each exit port as a function of the phase shift between the interferometer branches. Indeed, the fraction of atoms in the first exit port is given by ψ I (t) ψ I (t) and from Eq. (5) one can easily obtain with Therefore, δχ = 0 implies C = 1 and full oscillations, between 0 and N , in the total atom number at each exit port. In contrast, δχ = 0 causes a contrast reduction with smaller oscillations around N/2 of the atom number at each port, as illustrated in Fig. 2. (Completely analogous results hold for the second exit port.) Note that in general the expectation value ψ c (t) D (δχ) ψ c (t) is complex and its phase gives a contribution to the phase δφ in Eq. (7) in addition to the phase shift δΦ in Eq. (5). In order to gain further insight on the loss of contrast in open interferometers and help us devise a suitable mitigation strategy, it is instructive to consider the situation in position representation. Choosing the origin of coordinates so that R 1 = 0 to ease the notation, the probability density associated with the state in Eq. (5) becomes where some spatially independent phase terms (i.e. independent of x) have been absorbed in the phase δΦ to make the expression more compact. It is particularly illuminating to consider the result for free evolution at sufficiently late times, which is analyzed in Appendix C and is relevant for typical long-time interferometry applications. Using the state in Eq. (C3), Eq. (9) becomes where some spatially independent phase terms have again been absorbed in the phase δΦ and ∆t = t−t 0 . For simplicity we have considered here the case in which there is no spatially dependent contribution to the phase from the factorψ(m x/∆t, t 0 ) in Eq. (C3). Moreover, we have neglected the shift by δR of the envelope, which is a reasonable approximation whenever the envelope ψ c (x, t) does not vary too rapidly in space and its size is much larger than δR. (This approximation, which is often an excellent one in typical situations of interest, will be discussed in detail in Secs. III C and IV B.) Eq. (10) reveals the existence of a fringe pattern in the density profile of each exit port with a fringe spacing λ fr characterized by Integrating Eq. (9) over space, one gets the analog of Eqs. (7)-(8) in position representation with where we have shifted the integration variable x by δR/2 to write the result in a more convenient form. Specializing to the state in Eq. (C3) and making use of the same assumptions employed when deriving Eq. (10), which includes neglecting the shift of the envelopes, Eq. (12) be- (13) It is then easy to see that when the fringe spacing is much larger than the size of the envelope, as shown in Fig. 3a, the port goes from almost completely bright to almost completely dark as δφ changes, and the integrated particle number at each port oscillates from 0 to N (with C ≈ 1). On the other hand, for fringe spacing much smaller than the size of the envelope, as shown in Fig. 3b, C 1 and the integrated particle number at each port exhibits small amplitude oscillations around N/2. Besides providing an intuitive understanding for the loss of contrast in open interferometers, having identified the existence of a fringe pattern in the density profile at each exit port as the culprit can serve as a useful guide in the search for a mitigation strategy, which should be based on trying to eliminate those fringes.

Asymmetric pulse timing
We start by considering the simple example of an asymmetric Mach-Zehnder interferometer (MZI), where 1 The relative displacement between the envelopes as well as the possibility of having spatially dependent contributions to the phase fromψ(m x/∆t, t 0 ) and the corresponding corrections to the fringe spacing (11) will be taken into account when studying the evolution of Gaussian states and of expanding BECs in Secs. III and IV respectively.
FIG. 3. The two upper pictures depict the density profile at the exit port for δφ = 0 and δφ = π when the size of the fringes (dashed dark line) is much larger than the envelope (dashed light line); the port changes from almost completely "bright" to almost completely "dark" as far as the integrated particle number is concerned. The lower ones are analogous pictures when the fringe spacing is much smaller than the size of the envelope; the integrated particle number at the exit port changes very little as δφ varies.
the time separations between the mirror pulse and the beam-splitter pulses (initial and final) differ by δT , as shown in Fig. 4. In the absence of a quadratic potential there is only a position displacement between the two wave packets at the exit port and it is simply given by This gives rise to a fringe pattern in the density profile with fringe spacing as observed for instance in the experiments reported in Refs. [24,32]. The situation is somewhat analogous to what happens for rotations, but the displacement in that case is along the transverse directions (i.e. perpendicular to v rec ); see Sec. V below for a more detailed discussion.

Gravity gradient
Let us consider now the effect of the gravity gradient when the dynamics is governed by a Hamiltonian like that in Eq. (1). Before doing so we should, however, make a couple of important remarks. Firstly, although some time will elapse in general between the last beam splitter and detection (this time should be long enough when no state labeling is employed and good spatial separation between the two ports is, thus, required), the calculation of the contrast using Eq. (8) can be equivalently carried out right after the last beam splitter. This can be seen by rewriting the expression for the contrast as follows: where we took into account that ψ c (t) = U 0 (t, t bs ) ψ c (t bs ) and made use in the second equality of Eq. (B5) applied to this case together with Eq. (A11). Therefore, unless stated otherwise, from now on we will always calculate the displacement δχ at time t bs , right after the last beam splitter. Secondly, as mentioned in the introduction, we will neglect the effect of the gravity gradient on the evolution of the centered wave packet and use the free-particle unitary evolution operator for Using instead the full unitary evolution operator would simply add small corrections to the result for the contrast suppressed by additional powers of Γ T 2 .
Particularizing the results of Eqs. (A12)-(A13) to the case of a MZI with half-interferometer time T and δT = 0, one obtains the following position and momentum displacements at time t bs and to first order in Γ T 2 : These results are valid for an arbitrary orientation, but in the next subsection and Secs. III-IV we will concentrate on the case of aligned gravity gradients, i.e. the case in which the direction of the laser pulses and, hence, v rec coincide with the direction of one of the principal axes of the gravity gradient tensor Γ. This in turn implies δR v rec and δP v rec , which holds to first order but also for the exact results derived in Appendix A. On the other hand, the case of nonaligned gravity gradients will be addressed in Sec. V.
FIG. 5. Plot illustrating the effect of a quadratic potential on the wave-packet trajectories in a Mach-Zehnder interferometer (purple lines) including the momentum kicks from the laser pulses (green arrows) -the trajectories in the absence of quadratic potential terms (dashed blue lines) are also shown for comparison. As a consequence, the two trajectories at each exit port differ in position and momentum (they have different slopes in the plot).

D. Mitigation strategy
In order to mitigate the loss of contrast due to gravity gradients, one can try to combine the two examples described in the two previous subsections. Changing the timing of the last beam-splitter by a small amount δT in the presence of a gravity gradient leads for the MZI of the previous subsection to the following position displacement to first order in Γ T 2 and δT : but leaves δP invariant at that order. It is, therefore, not possible to eliminate completely the phase-space displacement δχ. Nevertheless, the key insight is that by choosing appropriately δT , one can get a nonvanishing value for δR that compensates the effect due to the nonvanishing δP and causes the fringe spacing to be much larger than the size of the envelope. Indeed, we see from Eq. (11) that 1/λ fr → 0 when and this can be achieved for aligned gravity gradients by taking where Γ is the eigenvalue of the Γ tensor along the direction of v rec and we have considered a total time t − t 0 = 2T + T 0 , with T being the half-interferometer time and T 0 the time from t 0 till the first beam splitter. The implementation of this mitigation strategy will be analyzed in more detail in the next two sections for freely evolving Gaussian states and for expanding BECs.

III. GAUSSIAN WAVE PACKETS FOR FREE PARTICLES
In this section we will provide a more detailed analysis specialized to Gaussian states of the questions discussed in Sec. II. This is especially interesting because one can obtain exact results which are fairly simple and can provide additional insight on certain aspects of the loss of contrast in open interferometers and the proposed mitigation strategy. In particular, they illustrate the intuitive interpretation of our findings within a phase-space description of quantum mechanics based on the Wigner function [33,34]. Furthermore, the study of Gaussian states is also of practical interest because they can be directly applied to the description of long-time interferometry experiments using cold thermal atoms [20]. The case of BECs, on the other hand, will be analyzed in the next section.

A. Phase-space description
Let us consider a general pure Gaussian state as the initial state of the centered wave packet, |ψ c (t 0 ) . Its associated Wigner function, defined in general as takes the form where ξ = (x, p) T and Σ is the phase-space covariance matrix, which is directly related to the Weyl-ordered twopoint functions: These are symmetric 3 × 3 matrices that can be regarded as blocks of the 6 × 6 covariance matrix Σ, which is symmetric and positive definite. In order to calculate the contrast, it is useful to note that the expectation value of a displacement operator can be obtained by computing the inverse symplectic Fourier transform of the Wigner function as follows [33]: where δχ 0 ≡ T (t 0 , t) δχ with the transition matrix T (t 2 , t 1 ) defined in Appendix A. In the second equality we have introduced the change of variables ξ = T (t, t 0 ) ξ, used the invariance of the symplectic form T T (t, t 0 ) J T (t, t 0 ) = J and taken into account that for quadratic potentials the Wigner function evolves exactly in the same way as a classical phase-space distribution. From Eqs. (27) and (23) we get the following result for the contrast defined in Eq. (8): Up to a factor 1/2 2 the exponent can be rewritten as In order to write δχ 0 in terms of δχ one needs to use the transition matrix T (t 0 , t), which is exactly given by Eq. (A8) for timeindependent gravity gradients. However, as explained in the Introduction and in Sec. II C, in typical cases of practical interest it is perfectly justified to treat Γ perturbatively and neglect its effect on the evolution of the centered wave packet |ψ c (t) . We will do so here, even though the exact result could be straightforwardly obtained, because the expressions are simpler and more transparent and the corrections would anyway be very small. For free evolution we have Thus, from Eqs. (29)-(31) one can see that for a given δP the contrast is maximized when δR Note that since the exponent in Eq. (28) is manifestly negative for any δχ 0 and one can always choose δR 0 so that δR (s) 0 = 0 for any given δP 0 , the second term in expression (29) must be negative. The interpretation and importance of this term will be discussed in Sec. III C. In this respect, it is instructive to consider the one-dimensional case, for which it becomes which follows from the inequality (det Σ) ≥ 2 /4. This inequality is the necessary and sufficient condition for the Gaussian phase-space distribution in Eq. (23) to be the Wigner function of a well-defined density matrix (witĥ ρ 2 ≤ρ); the equality holds for pure states and the strict inequality for mixed ones. Moreover, it coincides with the Schrödinger uncertainty relation [35,36], a somewhat stronger version of Heisenberg's uncertainty relation.
There is an alternative expression for Eq. (27) as the overlap of two Wigner functions with a relative phasespace displacement δχ (valid only for pure states): where the density matricesρ andρ δχ correspond respectively to the states ψ c (t) andD(δχ) ψ c (t) . Making use of Eq. (34) to interpret the contrast as the overlap between the two shifted Wigner functions can be quite illuminating and provides an intuitive explanation, within the phase-space formulation of quantum mechanics, for the loss of contrast and the reason why the mitigation strategy works. Indeed, as shown in Fig. 6, a nonvanishing relative displacement leads in general to decreasing overlap between the displaced Wigner functions (for simplicity we focus again on the one-dimensional case). Furthermore, it is clear that for fixed δP the overlap is maximized when δR is chosen so that the two displaced Wigner functions are aligned 2 . Taking into account that Σ −1 pp Σ xp characterizes the "tilt" in phase space of the initial Wigner function and considering its free evolution, one can easily obtain condition (32) for maximum contrast. It is also clear that substantial overlap will be achieved as long as δP (Σ pp ) 1/2 . Finally, it should be noted that for generic (non-Gaussian) pure states evolution for a sufficiently long time will lead to significantly squeezed Wigner functions and arguments analogous to those just presented for Gaussian states will still hold at the qualitative level.
The result in Eq. (32) agrees at late times with the late-time result (20) derived in Sec. II. Furthermore, the second term on the right-hand side of Eq. (32), which is proportional to Σ xp and suppressed by a factor 1/∆t compared to the first one, would have also been obtained in Sec. II if we had considered spatially dependent phase contributions fromψ(m x/∆t, t 0 ), which are always present for Σ xp = 0. . The loss of contrast simply corresponds to the reduced overlap of these two phase-space distributions. Provided that δP 2 Σpp, one can adjust δR with a slight timing asymmetry δT of the laser pulses so that the two Wigner functions (A and C) are properly aligned and almost complete overlap is recovered. For simplicity a one-dimensional example is depicted.

B. Pure versus mixed states
There is an interesting generalization of the previous results to the case of mixed states, as we explain below. First, one writes the initial density matrix aŝ so that Tr ξρ c (t 0 ) = 0. Next, we expressρ c (t 0 ) as an incoherent mixture of pure states and write the density matrix for the initial state in the following way: with 0 ≤ p j ≤ 1 and j p j = 1. One can then apply the results obtained for pure states to each state ψ (j) (t 0 ) and finally take the average over the whole ensemble. In doing so, it is convenient not to impose the condition that the expectation value of the operatorξ vanishes for centered states (this was done to specify them completely, but we needn't have done so and all the results derived for pure states would still hold since no further use of the condition was made in any of the calculations and derivations). One can then use this freedom to choose ψ (j) (t 0 ) as centered states. With this choice the term C cos δφ appearing in Eq. (7) corresponds for each pure state of the ensemble to the real part of where C (j) is defined as the modulus of the expression on the right-hand side (recall that δφ (j) contains a possible additional phase from the expectation value of the displacement operator). The relative displacement δχ is the same for all members of the ensemble, and so is δΦ, which depends on the initial displacement χ 0 appearing in Eq. (35), but that is also the same for every member. After taking the average over the whole ensemble, one is naturally led to the following quantity: whose real part characterizes the oscillations of the total particle number at the exit port for the whole ensemble, as given by Eq. (7). In particular, for the contrast one gets where one has a strict inequality unless all the phases δφ (j) are equal. Thus, we see that dephasing between the different members of the ensemble can lead to a further reduction of the contrast for the oscillations of the total particle number at each port as a function of δφ, even when there is hardly any loss of contrast for each member of the ensemble. It is clear from Eq. (39) that the result for pure states in Eq. (8) admits the following natural generalization to mixed states: which in turn implies that Eq. (27), in terms of the Wigner function, can be directly applied to mixed states as well.
[The Wigner function for mixed states is defined by replacing ψ c ψ c withρ c in Eq. (22).] In contrast, Eq. (34), which depends nonlinearly on the Wigner function, only holds for pure states. It should also be noted that if we had kept the condition of vanishing expectation value of the operatorξ for the centered states, this would have led to a different initial displacement χ (j) 0 for each pure state ψ (j) (t 0 ) and, hence, a different δΦ (j) too. However, this would be exactly compensated by the expectation value of the displacement operator with respect to the new centered states ψ Finally, it should be emphasized that so far the whole discussion and conclusions about mixed states, including the validity of Eq. (27), were completely general and not just restricted to Gaussian states. When particularizing to Gaussian density matrices, one can further conclude that Eq. (28), derived for pure Gaussian states in the previous subsection, can be directly applied to mixed states as well.

C. Quantitative examples
In order to illustrate how effective our proposed mitigation strategy can be, we consider the effects of an aligned gravity gradient with Γ zz = 3 · 10 −6 s −2 , comparable to that on the surface of Earth, and a pure Gaussian state for 87 Rb atoms with Σ zz = (100 µm) 2 and a momentum width along this direction corresponding to an effective temperature T eff = 1 nK (implying a velocity width of 0.3 mm/s). This temperature, which has already been successfully implemented in interferometry experiments employing the so-called delta-kick cooling technique [24], is much less demanding than the requirement of T eff ≈ 70 pK in current plans for STE-QUEST, mainly driven by the loss of contrast due to gravity gradients when no mitigation strategy is employed. (For simplicity we assume vanishing nondiagonal elements involving the z direction and the perpendicular directions within the matrix Σ, so that those transverse directions can be trivially integrated out and the example reduces to a one-dimensional problem. Moreover, we assume that the initial conditions are specified right before the first beam splitter, i.e. we take T 0 = 0.) In Fig. 7 one can see how the gravity gradient leads to a substantial loss of contrast for T = 5 s of half-interferometer time, whereas the mitigation strategy enables an extension beyond T = 20 s with virtually no loss of contrast. This requires perfect alignment and perfect knowledge of the gravity gradient, but even the rather conservative assumption that this can only be achieved at the 10% level, so that a residual 10% of the original gravity gradient remains, extends the possible half-interferometer time to close to T = 10 s.  7. Integrated contrast as a function of the halfinterferometer time for pure Gaussian states with Σzz = (100 µm) 2 , Σpp corresponding to T eff = 1 nK (blue curves) and the corresponding nonvanishing value of Σzp. One can see a substantial loss of contrast caused by an aligned gravity gradient with Γzz = 3 · 10 −6 s −2 (continuous blue line) and how this improves dramatically when the mitigation strategy is employed (dashed line). The result for a 10% residual gravity gradient is also shown for comparison (black curve).
In the previous example Σ zp had a nonvanishing value which was determined (up to a sign) by the requirement that (det Σ) = 2 /4 for pure Gaussian states. As a second example we will consider a Gaussian mixed state under the same conditions but with Σ zp = 0 as well as Σ zz = (141 µm) 2 and Σ pp corresponding to T eff = 2 nK. (For this choice the mixed state can be regarded as an incoherent ensemble of pure states like that in the first example.) As can be seen in Fig. 8, the mitigation strategy is somewhat less effective in this case and the contrast is reduced almost completely for T 10 s. This has a simple intuitive explanation in terms of the dephasing between the different members of the ensemble of pure states discussed in Sec. III B and leading to the inequality in Eq. (39). In fact, the Wigner function associated with a mixed Gaussian state can be written as a convolution of identical Wigner functions for pure Gaussian states but centered at different points in phase space, and whose centers follow a Gaussian distribution. Such a decomposition is not unique (there are infinitely many of them), but for the second example the Wigner functions in this convolution can be chosen to coincide with the pure state in our first example [37]. The loss of contrast and the effect of the mitigation strategy for each member of the ensemble are then the same as in that example. However, there is an additional reduction of contrast due to relative dephasing between the members of the ensemble and this remains even when the mitigation strategy is employed.
FIG. 8. Same as Fig. 7 but for a mixed Gaussian state with Σxp = 0, Σzz = (141 µm) 2 and Σpp corresponding to T eff = 2 nK, which can be understood as an incoherent ensemble of pure states like the one considered there. The faster decrease of the integrated contrast compared to Fig. 7 (especially when the mitigation strategy is employed) is due to the dephasing between the different members of the ensemble.
We finish this section by discussing the role and interpretation of the term quadratic in δP appearing in Eq. (29), which is entirely responsible for any remaining contrast reduction when the mitigation strategy is used. For pure states this is due to the lack of full overlap between the two displaced Wigner functions, even after they have been aligned, because δP is fixed. This can be easily illustrated for the one-dimensional case, where the ratio between the relative displacement of the aligned Wigner functions and their width along the alignment direction equals δP/ 2 Σ pp , and minus one half the square of this quantity [as in Eq. (33)] gives, when exponentiated, the overlap between two such Gaussian Wigner functions. In position representation this is directly related to the partial overlap of the envelopes of the two wave packets, whereas the phase-space alignment essentially corresponds to the lack of oscillations in the probability density of the wave-packet superposition. Substituting Eq. (18) into Eq. (33), one gets a relation for pure states from which one can immediately see how small this term is in typical situations of practical interest (the smaller, the better since it appears as a negative exponent): Thus, for T = 10 s and typical values of the gravity gradient on Earth's surface this term is very small all the way down to σ v ∼ 10 −3 v rec , which corresponds to an extremely narrow velocity distribution (and T eff ∼ 1 pK for 87 Rb). On the other hand, as seen in Eq. (33), for mixed states the result in Eq. (41) is enhanced by a factor 4 (det Σ)/ 2 > 1 compared to a pure state with the same Σ pp , and this characterizes the strength of the relative dephasing between the members of the ensemble mentioned above.

IV. EXPANDING BECs
A. Analytical results within the scaling approach Throughout this section we will consider an expanding BEC evolving through an interferometer sequence after being released from a trap [24,38], and will assume that it is dilute enough by the time the first beam-splitter pulse is applied so that the effect of interatomic interactions can be ignored when describing the evolution throughout the interferometer. This means that during that period the nonlinear terms in the Gross-Pitaevskii equation (D1) governing the dynamics of the condensate can be neglected and it reduces to the Schrödinger equation for single atoms. One can then directly apply the methods and results of Sec. II by taking an initial state |ψ c (t 1 ) which is obtained after evolving the wave function of the condensate with the Gross-Pitaevskii equation up to a time when nonlinearities can be neglected. In this context, the generalization of the Thomas-Fermi approximation to time-dependent potentials which naturally arises within the so-called scaling approach, briefly reviewed in Appendix D, provides a simple analytical result for the wave function of an expanding BEC after release from the trap: where the different objects appearing in this expression are defined in Eqs. (D5), (D8)-(D9) and (D11). Furthermore, it also describes correctly the free evolution at arbitrarily late times, as shown in Appendix D 3. Therefore, one can directly employ Eq. (42) without the need to match the solution for the early nonlinear regime to the free evolution governed by the Schrödinger equation at late times. Note that since |ψ TF | 2 is normalized to the number of atoms in the BEC, we have divided by √ N so that ψ c agrees with the single-particle normalization employed in Secs. II-III.
Substituting the scaling solution (42) into Eq. (9), one obtains the following result for the probability density, analogous to Eq. (10): where we have assumed that the size of the envelope (the rescaled Thomas-Fermi radius) is much larger than δR and have neglected the relative shift between the envelopes. The mitigation strategy is completely analogous to that already discussed in Sec. II D for the free-particle case. From Eq. (43) we see that the limit of large fringe spacing (larger than the size of the envelope) corresponds to the condition This can be achieved with a small shift δT in the timing of the last beam-splitter pulse, which can be easily determined by substituting Eqs. (17)-(18) into Eq. (44) and gives a result similar to Eq. (21). At late times Λ ii ≈ c i + b i (t − t 0 ) for i = 1, 2, 3 in the basis where the original trap diagonalizes, and we get the following expansion: Before providing several quantitative examples in the next subsection, we derive some useful analytic expressions for the integrated contrast of interferometers based on expanding BECs. It can be calculated by substituting the scaling solution (42) into Eq. (12). If one neglects δR compared to the size of the envelope, which is a very good approximation in most cases of practical interest (as illustrated by Fig. 9 and further discussed below), one gets the following result, analogous to Eq. (13): where ζ are the rescaled coordinates defined in Eq. (D3) and it has been taken into account that ψ TF (ζ), given by Eq. (D11), is a real and even function. We have also made use of the equality Λ TΛ =Λ T Λ, proved at the end of Appendix D 1.
For v rec aligned with a principal axis of both the trapping potential and the gravity gradient tensor Γ, one can choose a Cartesian coordinate system where Λ = diag Λ , Λ where R TF is the Thomas-Fermi radius along the longitudinal direction and 2π/λ fr = δP Λ − m δRΛ / withλ fr being the fringe spacing in the rescaled coordinates. The result in Eq. (47) can be easily generalized to nonaligned gravity gradients and traps. First, one determines λ fr and the longitudinal direction from the modulus and direction of the vector contracted with ζ within the argument of the cosine in Eq. (46) (note that the matrix Λ(t) is not necessarily diagonal anymore). Next, one integrates the two transverse directions, which results in the 1-D marginal distribution of the density distribution ψ 2 TF (ζ). As shown in Ref. [39], the form of this marginal distribution is the same as that appearing in Eq. (47) even for nonaligned traps: one simply needs to replace the square of the Thomas-Fermi radius (R TF ) 2 with the Σ component of the matrix Σ characterizing the ground state in Eq. (D11). Thus, Eq. (47) is also applicable to the general case with this simple replacement.

B. Quantitative examples
In this subsection we provide several quantitative examples of contrast loss due to gravity gradients and the outcome of the mitigation strategy for the case of expanding BECs. All our results are based on the approximate scaling solution (42).
Similarly to what was done in Sec. III C for Gaussian states, we will consider an aligned gravity gradient and trap with Γ = 3·10 −6 s −2 . In addition, we will characterize our BEC of 87 Rb atoms at some time t 0 when most of the repulsive interaction energy has been converted into kinetic energy and the expansion dynamics has entered the linear regime, with Λ(t) = C + B (t − t 0 ). Furthermore, since the result for the contrast in Eq. (47) only depends on the product R TF Λ (t) (and its time derivative), one can entirely parametrize the system in terms of R TF Λ (t 0 ) and R TFΛ (t 0 ), which determine the size of the condensate and the width of its velocity distribution at that time. The corresponding results are plotted in Fig. 9 (with and without the mitigation strategy) for an initial half-width R TF Λ (t 0 ) = 60 µm and two different widths of the velocity distribution: 0.3 mm/s and 0.08 mm/s (with associated effective temperatures of 1 nK and 70 pK).
FIG. 9. Integrated contrast as a function of the halfinterferometer time for an expanding BEC of 87 Rb atoms with T eff = 1 nK (blue) and T eff = 70 pK (red). The loss of contrast is due to a gravity gradient Γ = 3·10 −6 s −2 . The curves for the cases with mitigation strategy (dashed lines) and without (continuous lines) correspond to the analytic computation neglecting the shift between the envelopes. For comparison, the dots are the result of numerically evaluating the integral along the longitudinal direction when the relative shift of the envelopes is taken into account.
From Eq. (46) and as seen in Fig. 9, it is clear that the loss of contrast due to an aligned gravity gradient is entirely removed by the mitigation strategy when the shift between the envelopes is neglected. Indeed, when using the approximate scaling solution (42), the only source of contrast loss is the relative shift between the envelopes of the two wave packets in the superposition at each exit port. Thus, in order to capture this effect, which is small but eventually appears at sufficiently late times, one needs to consider the result of substituting the scaling solution (42) into Eq. (12) without neglecting the shift between the two envelopes: which is an exact result for the centered wave function ψ c (x, t) given by Eq. (42). In deriving Eq. (48) we have taken again into account that ψ TF (ζ) is a real and even function. For aligned gravity gradients and traps one can proceed similarly to what was done above to obtain Eq. (47). While the integrals with respect to the transverse coordinates in Eq. (48) can also be performed analytically, in this case there is no simple analytical result for the integral with respect to the longitudinal coordinate ζ due to the more complicated functional dependence of ψ TF ζ + Λ −1 δR/2 ψ TF ζ − Λ −1 δR/2 compared to ψ 2 TF (ζ). This last integral can be evaluated numerically, which was done to obtain the results represented as dots in Fig. 9. Given their agreement with the continuous curves, the figure shows that neglecting the relative shift of the envelopes is indeed a very good approximation for situations of practical interest, as already discussed at the end of Sec. III C. One can also see that (within the time-dependent Thomas-Fermi approximation) this relative shift is the only source of contrast loss due to gravity gradients when the mitigation strategy is employed. For the examples plotted in Fig. 9 the effect is very small and can only start to be appreciated for very long halfinterferometer times. Moreover, it is slightly more significant for lower T eff because the size of the envelopes at late times is smaller in that case.

V. ROTATIONS AND NONALIGNED GRAVITY GRADIENTS
Besides asymmetric pulse timing and gravity gradients, another effect leading to open MZIs are rotations. Given an otherwise closed MZI in a uniform (and possibly timedependent) force field, changes in the direction of the wave vector k i for the different laser pulses will lead in general to Since we are working in the nonrotating frame as far as the dynamics of the atoms is concerned, the effect of rotations is to rotate the effective momentum transfer associated with each laser pulse [31]. (It can also cause a rotation of g(t), but that does not alter the displacement, which is independent of g(t) as explained in Appendix A.) For example, for a MZI a uniform angular velocity Ω produces the following displacements after the last beam splitter: Thus, to lowest order in Ω T we have a position displacement perpendicular to v rec and no momentum displacement. Once again this leads to the appearance of a fringe pattern in the density profile at each exit port and a loss of contrast in the oscillations of the integrated particle number as described in Sec. II. A possible way of dealing with this problem is to use a tip-tilt mirror that can change the direction of the momentum transfer associated with each laser pulse and compensate the effect of rotations, so that δR vanishes to first order in Ω T . This scheme was proposed in Ref. [26] and has been implemented by a number of groups to overcome the effects of Earth's rotation in interferometry experiments with atomic fountains [20,[27][28][29]. For nonaligned gravity gradients, i.e. when v rec does not coincide with a principal direction of the gravity gradient tensor Γ(t), Eqs. (17)- (18) imply that δR ∦ v rec and δP ∦ v rec . Fortunately, the mitigation strategy presented in the previous sections can be extended fairly easily to this new situation. The basic idea is to use the same kind of tip-tilt mirror scheme described above to generate an additional position displacement along the transverse directions (perpendicular to v rec ) so that, together with the displacements generated by rotations and gravity gradients, the transverse projection of condition (20) is fulfilled: On the other hand, analogously to what was done for aligned gravity gradients, one can change slightly the timing of the last beam-splitter pulse to guarantee that the projection of Eq. (20) along the longitudinal direction, is also satisfied. (In this case, in addition to those due to the gravity gradient and the pulse timing asymmetry, one may also need to include possible contributions to longitudinal displacements due to higher-order terms generated by rotations and the tip-tilt mirror.) Finally, it should be noted that although illustrated here with condition (20), the strategy can be straightforwardly applied to more general situations corresponding to conditions (32) and (44) as well.

VI. CONCLUSIONS
We have analyzed in detail the reduction of the integrated contrast in open atom interferometers, with particular emphasis on the case of gravity gradients, and presented a simple mitigation strategy. For aligned gravity gradients it merely involves a slight change in the timing of the last beam-splitter pulse, whereas in the nonaligned case it can be combined with schemes based on the use of tip-tilt mirrors which have been employed so far to deal with transverse relative displacements caused by rotations. As we have shown, the strategy is very effective and in principle it enables extending the halfinterferometer time beyond T = 20 s with hardly any loss of contrast due to gravity gradients. Its basic limitation is the decreasing overlap between the envelopes of the two wave packets in the superposition, which eventually becomes important for sufficiently long interferometer times. Nevertheless, this should not be a major concern in the short or mid-term future since even for very narrow momentum distributions (leading to a smaller size of the envelope at long times) with σ v ∼ 10 −3 v rec , which corresponds to T eff ∼ 1 pK for 87 Rb, the effect due to a decreasing overlap of the envelopes is still very small for half-interferometer times of T = 10 s.
An important advantage of our proposed strategy is its ease of implementation. For aligned gravity gradients it only requires a control of the pulse timing which is well within the capabilities of current experimental set-ups. (Moreover, if that were not the case, the imperfections in controlling the pulse timing would lead themselves to a comparable loss of contrast anyway.) On the other hand, for nonaligned gravity gradients one needs to use a tip-tilt mirror scheme in order to deal with the transverse displacements, but the required degree of control has already been demonstrated in experiments that employ such a scheme for compensation of rotations (and any limitations in this respect would result, if the set-up is also being used to compensate rotational effects, in a loss of contrast due to gravity gradients comparable to those caused anyway by the remaining non-compensated rotations). Thus, we see that the implementation is indeed much simpler (with essentially no new requirements) than for possible alternative approaches such as trying to engineer the initial state while keeping all the desirable features for long-time and high-precision interferometry. Furthermore, the relative change in the pulse timing given by Eq. (21), which is of order Γ T 2 , gives rise to additional terms in the expression for the phase shift δφ of the same kind and magnitude as terms already present for δT = 0. Hence, these contributions should be small enough or be taken care of by any mitigation techniques already in place, so that the required δT does not lead to any significant reduction in sensitivity nor accuracy.
In order to apply our mitigation strategy, one needs of course to know the gravity gradient beforehand. This can be measured by other means, determined from geodetic calculations and simulations of the mass distribution near the experimental set-up or even from a calibration of δT which eliminates the fringe pattern in the density profile of the exit ports (provided that the gravity gradient does not change much over time or it does so slowly enough).
Fortunately, one does not need to know the gravity gradient very precisely: even with a 10% inaccuracy (which is a rather conservative assumption) the method can effectively deal with the loss of contrast due to gravity gradients up to half-interferometer times close to T = 10 s, as shown in Fig. 7. At this point it should be mentioned that by considering multiloop (e.g. with three intermediate π pulses) and an appropriate choice of the time between the different pulses (independent of Γ) one can achieve a vanishing final displacement δχ up to a certain order in Γ T 2 (and Ω T ) for time-independent (or very slowly varying) gravity gradients [40]. However, for such configurations the lowest-order contributions to the phase shift from time-independent g (and Γ) cancel out too. These are the relevant ones for accelerometers and gravimeters as well as for tests of the weak equivalence principle (the universality of free fall) and gradiometers, where one performs differential acceleration measurements for different species or at different spatial locations. Hence, although such configurations exhibit no significant loss of contrast due to (time-independent) gravity gradients, the typical signals of interest are severely suppressed as well and they end up being of little use for that kind of applications.
Our results are valid for atom interferometers based on Raman or Bragg pulses, and they can be straightforwardly extended to double diffraction schemes [41,42]: in most of the results one simply needs to double the effective momentum transfer associated with each laser pulse. In addition, they can be applied to MZIs with several intermediate π pulses and easily generalized to configurations with more than two beam-splitter pulses (such as Ramsey-Bordé interferometers) as long as only two wave packets contribute to the superposition at each exit port and these are completely distinguishable (thanks to good spatial separation or internal-state labeling). Furthermore, as a side result we have also provided a remarkably compact derivation of the general expression for the phase shift in Appendix B, where a representationfree [43,44] description of the state evolution in a lightpulse atom interferometer has been derived. This has been obtained for quadratic Hamiltonians including timedependent uniform forces (or accelerations) and gravity gradients, characterized by g(t) and Γ(t). However, it can be easily generalized to anharmonic potentials which can be locally approximated (within regions of the wavepacket size) by quadratic ones [30]. Moreover, although the effects of rotations have been described in a nonrotating frame [31], the Hamiltonian in a rotating frame would still be quadratic and the results in Appendix B would also be applicable.
We close this section by discussing a recent proposal made in Refs. [20,29] to extract useful information for precision interferometry directly from the features of the fringe pattern in the density profile at the exit ports. The basic idea is that by fitting the observed density profile with a formula like Eq. (10) and an appropriate envelope (typically a Gaussian), one can in principle extract information on the displacement δχ and the phase shift δΦ even for a single shot [24,29,32,45]. In more realistic situations various factors will lead to a reduction of the contrast in the fringe pattern, which will be better described by the following formula for the density profile: with C f ≤ 1. Interestingly, by fitting the density profile, one can estimate the phase shift and the fringe contrast C f for every single shot, which is not possible if one considers instead the integrated atom number at each exit port. As emphasized in Refs. [20,29], this can be particularly advantageous for portable inertial sensors in dynamical platforms, which may experience significant contrast fluctuations. The approach has been demonstrated with a measurement of Earth's rotation rate with a precision of 200 nrad/s [20] and a gyrocompass with 10 mdeg precision [29]. In both sets of experiments the measured quantities were derived from the observed fringe spacings; moreover, in the latter controlled changes of the final displacement were introduced to modify the fringe spacing and facilitate its read-out. Since it allows extracting information on the phase shift for nonvanishing displacements, fitting the fringe pattern of the density profiles also offers an alternative approach in those situations where gravity gradients induce a significant displacement. A careful analysis taking into account the specific details of each particular experimental implementation may be necessary in order to determine what is the preferable method for measuring the phase shift in each case. For instance, the simpler data analysis involved in measuring the integrated particle number and the fact that is a well-proven technique which has been employed in numerous high-precision atom-interferometry experiments to date are desirable features for a space mission (and this was the choice made in the design of STE-QUEST). On the other hand, the additional information contained in the density profile offers the possibility of a more thorough and continuous monitoring of (possibly unexpected) sources of noise and systematics.
The experiments reported in Refs. [20,29] make use of cold thermal atoms and can be described in terms of the mixed Gaussian states studied in Sec. III B. In order to discuss the point-source-interferometry (PSI) approximation considered in those references, it is convenient to decompose the initial Gaussian thermal state with Σ xp = 0 as an ensemble of pure Gaussian states with Σ xp = 0 and the same Σ xx as the thermal state but a smaller momentum width Σ pp . These pure states have different central momenta weighed according to a Gaussian distribution, but are otherwise identical. As we plan to discuss in more detail in a separate publication, the treatment of Refs. [20,29] naturally corresponds to the regime where the density profile at the exit ports for each member of the ensemble exhibits no fringe pattern (although a fringe pattern may arise for the whole ensemble due to the dephasing between the different members). In general, however, the PSI approximation has a number of limitations and can miss relevant effects associated with δP 0 which can, in contrast, be fully accounted for with the formalism that we have presented here [37].
For a time-independent gravity gradient tensor the transition matrix can be straightforwardly obtained by exponentiating the matrix H(t): (A8) where we have introduced γ ≡ √ Γ. In order to calculate explicitly the transition matrix in Eq. (A8) one needs to diagonalize the symmetric tensor Γ, which is always possible with an orthogonal transformation (a rotation of the coordinate axes). In this new coordinate system the motion along each one of the three axes (known as principal axes) decouples and the dynamics reduces to that of three independent one-dimensional systems. In general, some of the eigenvalues of Γ will be negative since its trace vanishes outside gravitational sources. For those directions with negative eigenvalues one needs to write γ jj = i ω jj in Eq. (A8) and the dynamics along that direction corresponds then to that of a one-dimensional harmonic oscillator, whereas for positive eigenvalues one can directly apply Eq. (A8), which describes the dynamics of an inverted harmonic oscillator. In typical cases of interest in atom interferometry the condition |Γ jj |T 2 1 is amply satisfied for the three principal axes and it is an excellent approximation to consider the following perturbative expansion of Eq. (A8) up to linear order in Γ (the expansion involves only even powers of γ): (A9) where we neglected terms of higher order in Γ(t − t ) 2 . A similar expression (but involving a time integral) can be obtained for a general time-dependent Γ(t) by solving perturbatively Eq. (A1).
Note that the first two terms on the right-hand side of Eq. (A5) are common for the two interferometer branches 3 , so that only the third term contributes to the relative displacement between the classical trajectories for the center of the wave packets in Fig. 1: where δF lp simply corresponds to the difference between the laser kicks, as given by Eq. (A4), for the two trajectories. This also implies that the evolution of the displacement after the last beam splitter takes the following simple form: where we used the definition of T ret in terms of the transition matrix and the composition of transition matrices.
This together with Eq. (A7) implies the following result for the displacement operator: where χ h (t) = T (t, t 0 ) χ 0 corresponds to the classical phase-space trajectory associated with the Hamiltonian H 0 and with initial condition χ h (t 0 ) = χ 0 . Taking all that into account and inserting the resolution of the identity I =Û † 0 (t, t 0 )Û 0 (t, t 0 ) after the displacement operator on the right-hand side of Eq. (B1), the evolution of that state generated by the purely quadratic part of the Hamiltonian can be written in the following suggestive way: with |ψ c (t) =Û 0 (t, t 0 ) |ψ c (t 0 ) and

Purely linear part and laser pulses
Let us consider the effect of the laser pulses first. Focusing on the external degrees of freedom (the internal state can be easily taken into account by appropriately labeling the segments of the trajectories after each pulse and considering only the interference between wave packets with the same internal state at the exit ports), the effect of idealized beam-splitter and mirror pulses can be characterized (for the n-th pulse) by a phase factor e iεnϕn times a displacement operator e iεnk T nx =D F n acting on every wave packet. The integer ε n , the laser phase ϕ n , the momentum transfer k n and the associated phase-space displacement F n [defined in Eq. (A4)] depend on both the pulse number n and the interferometer branch (the classical trajectory assigned to each wave packet). In addition, beam-splitter pulses generate a superposition corresponding to two of such operators with different values of ε n acting on the same wave packet. Given a state of the formD χ(t − n ) |ψ c (t n ) right before, the action of a laser pulse at time t n generates the following state: where we used the composition formula in Eq. (6). Evolution between pulses proceeds as described in the previous subsection, so that the state at some time t > t n before the action of any further pulse is given by where the displacement χ(t) corresponds to the classical trajectory with initial conditions χ(t 0 ) = χ 0 and including the kicks from the laser pulses: with F lp defined in Eq. (A4). The phase Φ(t) in Eq. (B9) is in turn given by with ϕ = n j=1 ε j ϕ j . In Eqs. (B10)-(B11) we have already provided the general result for n pulses. It can be easily proved by induction. First, one uses Eq. (B8) to prove the result for n = 1. Next, assuming that the result holds for n − 1 pulses, one can use Eq. (B8) again to show that it then holds for n pulses too.
Let us now turn our attention to the purely linear term in the Hamiltonian. The quickest way of generalizing the previous results in order to include the effect of the linear term is by regarding it as equivalent to the continuum limit δt → 0 of r pulses with F j = δt G(t j ), where G(t) is defined in Eq. (A3), t j = t 0 + j δt with j = 1, . . . , r and r δt = t − t 0 . Eq. (B10) is then generalized to Eq. (A5) and Eq. (B11) becomes

Phase shift
At the exit ports one has a superposition of several wave packets, which can each be evolved as described in the previous subsection. For interferometers where two wave packets contribute to each exit port, as depicted for instance in Fig. 1, the superposition is given by Eq. (5). The displacement δχ can be calculated as explained in Appendix A. Here we present a compact but complete derivation of the result for the phase shift δΦ.
First one takes the expression for δΦ provided after Eq. (5) and substitutes Φ 1 and Φ 2 by the result in Eq. (B12). (For the moment we do not consider the contributions from the purely linear term of the potential and take G = 0 so that the expressions are shorter, but they can be easily added at the end, as done in the previous subsection.) Next, one rewrites the result in terms of the semisum and difference variablesχ = (χ 1 + χ 2 )/2 and δχ = χ 2 − χ 1 . One also proceeds in the same way with the sets of pulse displacements F (1) j and F (2) j associated with the two classical trajectories as well as the phases ϕ 1 and ϕ 2 . Together with several additional steps described in more detail below, this leads to the following intermediate and final result: In writing the second line, we took the transpose of the second term of the integrand (it is always possible since it is a scalar), which led to a change of sign because J T = −J. We also wrote the last term on the first line as the time integral from t 0 to t of its time derivative, taking into account that δχ(t 0 ) = 0 for t 0 before the first beam splitter. Finally, the following steps were carried out in the last equality. First, the time derivatives were substituted using the relations δχ = H δχ + δF lp and δχ = H δχ +F lp that follow from Eq. (A1). Second, we used the identity H T J = −J H, which holds for any quadratic Hamiltonian system. Of the two remaining terms one cancels the second term in the integrand, whereas the other is identical to the first term in the integrand, which gets doubled.
We can now include the contributions from the purely linear term of the Hamiltonian by proceeding as done at the end of the previous subsection. All that one needs to do is to use the full solution in Eq. (A5), including G, withF lp in place of F lp . Eq. (B13) agrees then with the general result obtained by Antoine and Bordé in Refs. [47,48]. Furthermore, one can easily generalize the result to the case with different values of G for each branch. This can be relevant for instance when there is a uniform force whose value depends on the internal state of the atoms acting in an interferometer with different internal states in each branch. It is also useful for describing the effects of anharmonic potentials which can be locally approximated (within regions of the size of each wave packet) by quadratic ones [30]. The new result can be written in a compact form as whereχ(t ) is given by Eq. (A5) withF lp andḠ. Finally, it is sometimes convenient to write the result in an alternative way which explicitly displays all the dependence on the initial conditions. When doing so, Eq. (B14) becomes where we have used Eq. (A5) forχ(t ), the identity J T (t , t 0 ) = T T (t, t ) J T (t, t 0 ), which follows from Eq. (A7), and the expression for δχ T (t) in terms of δF T lp and δG T obtained from Eq. (A5).

Appendix C: LATE-TIME FREE EVOLUTION
It is well known that the free evolution of a wave packet in position representation takes a simple form at late times, as we will see. One starts by noting that the free evolution of a wave packet |ψ(t) can be straightforwardly calculated in momentum representation, i.e. for ψ(p, t) = p|ψ(t) . The evolution in position representation can then be directly obtained by an inverse Fourier transform: which can be conveniently rewritten as with ∆t = (t − t 0 ). The integral can be exactly evaluated in certain cases, e.g. for Gaussian wave packets. Moreover, for sufficiently late times one can use a saddle-point approximation (also known as stationary-phase approximation) and obtain the following result: This is a good approximation provided that the characteristic scales ∆p over whichψ(p, t 0 ) changes significantly are much larger than the width of the Gaussian peak in the saddle-point approximation, i.e. that ∆p i (m /∆t) 1/2 for i = 1, 2, 3. It is worth mentioning that the result just derived can also be easily obtained within a phase-space formulation in terms of the Wigner function. One simply needs to take into account that the Wigner function evolves like a classical phase-space distribution, following the trajectories of a free particle in phase space, and that this leads to a shear (squeezing plus tilting) of the initial distribution. The result is equivalent to Eq. (C3) once one takes into account the rescaling of the spatial coordinates on the right-hand side of that equation and the fact that the spatially dependent phase factor leads to a sheared Wigner function, as illustrated by Eqs. (D12)-(D13) below. Further details can be found in Ref. [30].

Appendix D: DESCRIPTION OF AN EXPANDING BEC
For most purposes the expansion dynamics of a BEC released from a (possibly time-dependent) trap can be satisfactorily described in terms of a mean-field approximation based on the Gross-Pitaevskii equation, which has the following form: (D1) where m denotes the mass of the atoms and the coupling constant g is related to the s-wave scattering length a s according to g = 4π 2 a s /m. Moreover, for (locally) quadratic potentials a separation of the dynamics of the center of the wave packet analogous to that described in Appendix B for non-interacting particles is also possible [30,39,49] and it matches quite naturally with the noninteracting case when the BEC becomes dilute enough and the effects of the non-linearities can be neglected. We will, therefore, focus throughout the rest of this appendix on the dynamics of a centered condensate evolving in a purely quadratic potential of the form where Ω 2 (t) is a symmetric and initially positive-definite matrix.

Scaling approach
For purely quadratic potentials there is a very useful approach for extending the Thomas-Fermi approximation to time-dependent potentials. Within the hydrodynamic description of the condensate, regarded as a superfluid, it essentially amounts to considering a rescaled coordinate system locally comoving with the expanding fluid, i.e. where every fluid element is at rest [50]. The approach can also be implemented directly at the level of the Gross-Pitaevskii equation. One starts by introducing the linearly rescaled coordinates where Λ(t) is a time-dependent matrix to be determined below. Next, one considers the following rescaling of the wave function: where the prefactor with the square root of the determinant guarantees a simple normalization condition for the new wave function ψ Λ (ζ, t) in terms of the rescaled coordinates ζ, namely d 3 ζ |ψ Λ (ζ, t)| 2 = N . The spatially dependent and independent phases Φ(ζ, t) and β(t) are given respectively by where As discussed below, with this choice the phase Φ(ζ, t) will account for most of the kinetic energy of the expanding BEC in the Thomas-Fermi regime. Substituting Eq. (D4) into the original Gross-Pitaevskii equation (D1) and demanding that Λ(t) satisfies the ordinary differential equation with initial conditions one finds that the equation governing the dynamics of ψ Λ (ζ, t) takes the following simple and convenient form: If one neglects the first term on the right-hand side of Eq. (D10), its solution for an expanding condensate initially in the ground state of the trapping potential becomes time independent and coincides with the Thomas-Fermi approximation of the ground-state wave function: where the constant µ is determined by the total number of atoms in the condensate and the subindex + has been introduced to denote that the function vanishes when the argument of the square root becomes negative. Having neglected the first term on the right-hand side of Eq. (D10) corresponds to neglecting the contribution to the kinetic term due to the gradient of the envelope of the wave function, which constitutes a natural generalization of the Thomas-Fermi approximation to the time-dependent case and is equivalent to neglecting the so-called quantum pressure term in the hydrodynamic description. In contrast with the ground state case, for an expanding BEC there is a significant contribution from the kinetic term in Eq. (D1) (one can easily understand from energy conservation that all the initial interaction energy associated with the nonlinear term is eventually transformed into kinetic energy), but most of it is accounted for by the spatially dependent phase (D6) provided that the de Broglie wavelength associated with the characteristic momenta of the expanding superfluid is much shorter than the size of the envelope (the rescaled Thomas-Fermi radius).
It is clear that only the symmetric part of the matrix A(t) contributes to Eq. (D6). Moreover, the conditions (D8)-(D9) fulfilled by Λ(t) guarantee that A(t) is symmetric itself. This can be shown as follows. First, one notes that the initial conditions (D9) imply A(t 0 ) = A T (t 0 ) = 0. Next, one can check that Eq. (D8) implieṡ A(t) −Ȧ T (t) = 0. It then follows that A(t) = A T (t) ∀ t.
Next, we briefly discuss he phase-space description for an expanding BEC, which provides useful additional insight. The Wigner function for the approximate scaling solution can be obtained by substituting Eq. (D4) into Eq. (22) and employing the time-independent approximation (D11) for ψ Λ (ζ, t). In the original unrescaled coordinates it is given by Restoring the rescaled variable ζ and introducing the analogously rescaled variable∆ = Λ −1 ∆ together with the momentum p ζ = Λ T p, it becomes This result has a very suggestive form and shows that the Wigner function at later times is given by a sheared version of the initial Wigner function for the Thomas-Fermi ground state with the shearing degree (the amount of tilting) controlled by the matrix A(t), which vanishes at the initial time. This shear is directly related to the spatially dependent phase factor in the wave function, with a quadratic dependence on the position of the exponent, and it follows from the conversion of the nonlinear interaction energy into kinetic energy. Furthermore, given the late-time behavior of A(t), linear in time, and rescaling back to the original coordinates {x, p}, one finds that the result is compatible with free evolution at late times and with the result of Appendix C. Finally, since most of the support of the Thomas-Fermi Wigner function W TF ζ, p ζ ) (which also exhibits smaller amplitude oscillations with negative values) corresponds to a region of spatial width 2R TF and momentum width ∆p ∼ /R TF , the time evolution leads to a squeezing and tilting of this region which looks qualitatively similar to the evolution of the Gaussian states considered in Sec. III. Thus, the interpretation of the loss of contrast and the mitigation strategy in terms of the alignment of squeezed distributions in phase-space also applies to the scaling solution for expanding BECs. And the same can be said about the generic late-time behavior described in Appendix C.

Momentum distribution at intermediate times
An interesting result within the scaling approximation, pointed out in Refs. [39,49], is the simple connection existing between the initial wave function in position representation for the atoms in an expanding BEC and the wave function in momentum representation at later times. This relationship is obtained by Fourier transforming the scaling solution (D4): ψ(p, t) = det Λ(t) and using the saddle-point approximation to obtaiñ which is valid when the widths of the Gaussian in Eq. (D15) are smaller than the characteristic scales ∆ξ of ψ TF (ξ), i.e. when /(m|A ii |) ∆ξ i ∼ R (i) TF for i = 1, 2, 3 in the basis where the matrix A ij diagonalizes.
It is interesting to consider the linear regime for the scaling matrix, where Λ(t) = C + B (t − t 0 ). Whenever the time-dependent Thomas-Fermi approximation is valid, this corresponds to the situation where almost all the nonlinear interaction energy has already been converted into kinetic energy. In this case the condition for the validity of the saddle-point approximation becomes R (i) TF /(m|B ii ||Λ ii |). [Here we considered the basis where the matrices B and Λ(t) both diagonalize, which is always possible for nonrotating traps and coincides with the basis where the matrix Ω 2 (t) diagonalizes.] Note, for comparison, that Heisenberg's uncertainty principle implies R /(m|B ii ||Λ ii |). In terms of the phase-space description of the condensate's mean field this amounts to having a sufficient tilt of the associated Wigner function, as can be seen from Eq. (D13). The tilt develops when the nonlinear interaction energy is converted into kinetic energy and it is always large enough for a BEC in the Thomas-Fermi regime which is released nonadiabatically. On the other hand, when using delta-kick cooling or a similar technique to reduce the momentum spread by rotating the Wigner function in phase space, the tilt is reduced (or even entirely removed) and one gets close to the saturation of Heisenberg's uncertainty relation. If this is done efficiently enough, the condition for the validity of the saddle-point approximation will no longer hold and the accuracy of Eq. (D16) will no longer be guaranteed. In that case, in fact, even the scaling approximation will cease to be valid if the condensate is dilute enough at that time because the expansion will then be mainly driven by the gradient of the envelope rather than by the kinetic energy that results from the conversion of the nonlinear interaction energy.

Late-time scaling approach versus free evolution
The result of the previous subsection can be combined with that of Appendix C to analyze the late-time behav-ior of an expanding BEC. In order to do so, we consider an intermediate time t 1 at which the density of the condensate is already low enough that we can neglect the nonlinear interactions. Hence, at this point the matrix Λ(t) has already entered the linear regime. One can use Eq. (D16) to evaluateψ(p, t 1 ), and substituting the expression for Λ(t) linear in time, one gets: ψ(p, t 1 ) ≈ (det B) − 1 2 m 3/2 e iβ e − i p 2 ∆t1/2m ψ TF B −1 p/m , (D17) with ∆t 1 = (t 1 − t 0 ). Furthermore, for times t ≥ t 1 the atoms in the condensate behave as free particles (since the interactions can be neglected) and the method of Appendix C can be employed. Making use of Eq. (C1) applied to this case and taking Eq. (D17) for the initial wave function at time t 1 in momentum representation, one gets ψ(x, t) = d 3 p (2π ) 3/2 e iβ e i p T x e − i (p 2 /2m)(t−t1)ψ (p, t 1 ) where in the second equality we combined the two phase factors with exponent proportional to p 2 and introduced ∆t = (t−t 1 )+∆t 1 = (t−t 0 ). Next, rewriting Eq. (D18) as done in Eq. (C2) and using the same kind of saddle-point approximation that led to Eq. (C3), one finally obtains: ψ(x, t) ≈ det(B ∆t) ψ TF B −1 x/∆t , (D19) which coincides with the scaling solution (D4) when C can be neglected compared to B ∆t. The saddle-point approximation is valid provided that the typical scales of variation ∆p of the functionψ(p, t 1 ) satisfy the condition m /∆t ∆p i , which implies /(m|B ii | 2 ∆t) |B ii | −1 ∆p i /m ∼ R (i) TF for i = 1, 2, 3. For sufficiently late times (i.e. sufficiently large ∆t) this condition will be fulfilled as long as the condition justifying the derivation of Eq. (D16) holds.
The result found in this subsection is particularly interesting because it means that the validity of the scaling approximation can be automatically extended to late times during the free expansion: if it is valid at some time t 1 at which the expansion dynamics has already entered the linear regime, it will be valid at arbitrarily late times.