Overcoming loss of contrast in atom interferometry due to gravity gradients

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 $\pi /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:

$$\begin{eqnarray}&&\hat{H}=\frac{1}{2m}{{{\bf \hat{p}}}^{{\rm T}}}{\bf \hat{p}}-\frac{m}{2}{{{\bf \hat{x}}}^{{\rm T}}}\Gamma (t){\bf \hat{x}}-m{{{\bf g}}^{{\rm T}}}(t){\bf \hat{x}}+{{V}_{0}}(t).\end{eqnarray} \tag{ 1 }$$

It can describe the effect of time-dependent gravity gradients and uniform force (or acceleration) fields, characterized respectively by the tensor $\Gamma (t)$ and the vector ${\bf g}(t)$, as well as different internal energies, characterized by the spatially independent potential ${{V}_{0}}(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 $\hbar {{{\bf k}}_{i}}$ associated with the different laser pulses, and possibly of $\Gamma (t)$ and ${\bf g}(t)$ too [38].

As shown in appendix B, for quadratic Hamiltonians like that of equation (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. (A similar construction for the free propagation between laser pulses has been considered in [33, 39, 40].) 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

$$\begin{eqnarray}&&\mid \psi (t)\rangle =\frac{1}{\sqrt{2}}\left[ {{{\rm e}}^{{\rm i}{{\Phi }_{a}}}}\hat{\mathcal{D}}({{{\boldsymbol{ \chi }} }_{a}})\mid {{\psi }_{{\rm c}}}(t)\rangle +{{{\rm e}}^{{\rm i}{{\Phi }_{b}}}}\hat{\mathcal{D}}({{{\boldsymbol{ \chi }} }_{b}})\mid {{\psi }_{{\rm c}}}(t)\rangle \right],\end{eqnarray} \tag{ 2 }$$

where $\mid {{\psi }_{{\rm c}}}(t)\rangle$ 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 operator $\hat{\mathcal{D}}({\boldsymbol{ \chi }} )$ is defined as

$$\begin{eqnarray}&&\hat{\mathcal{D}}({\boldsymbol{ \chi }} )={{{\rm e}}^{-\frac{{\rm i}}{\hbar }\,{{{\boldsymbol{ \chi }} }^{{\rm T}}}J{\boldsymbol{\hat{ \xi }}} }},\end{eqnarray} \tag{ 3 }$$

where we used a vector notation for phase-space quantities, so that ${\boldsymbol{\hat{ \xi }}} ={{({\bf \hat{x}},{\bf \hat{p}})}^{{\rm T}}}$, and introduced the symplectic form

$$\begin{eqnarray}J=\left( \begin{array}{ccccccccccccccc} 0 & {\mathbb{1}} \\ -{\mathbb{1}} & 0 \\ \end{array} \right).\end{eqnarray} \tag{ 4 }$$

The displacement vectors ${\boldsymbol{ \chi }} (t)={{({\boldsymbol{\mathcal{R}}} (t),{\boldsymbol{\mathcal{P}}} (t))}^{{\rm T}}}$ in equation (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 figure 1. Therefore, we have the following superposition at one of the exit ports (after the second beam splitter):

$$\begin{eqnarray}\begin{array}{rcl} \mid {{\psi }_{{\rm i}}}(t)\rangle & = & \frac{1}{2}\left[ {{{\rm e}}^{{\rm i}{{\Phi }_{1}}}}\hat{\mathcal{D}}({{{\boldsymbol{ \chi }} }_{1}})\mid {{\psi }_{{\rm c}}}(t)\rangle +{{{\rm e}}^{{\rm i}{{\Phi }_{2}}}}\hat{\mathcal{D}}({{{\boldsymbol{ \chi }} }_{2}})\mid {{\psi }_{{\rm c}}}(t)\rangle \right] \\ {} & = & \frac{1}{2}{{{\rm e}}^{{\rm i}{{\Phi }_{1}}}}\hat{\mathcal{D}}({{{\boldsymbol{ \chi }} }_{1}})\left[ 1+{{{\rm e}}^{{\rm i}\delta \Phi }}\hat{\mathcal{D}}(\delta {\boldsymbol{ \chi }} ) \right]\mid {{\psi }_{{\rm c}}}(t)\rangle , \\ \end{array}\end{eqnarray} \tag{ 5 }$$

with $\delta {\boldsymbol{ \chi }} ={{{\boldsymbol{ \chi }} }_{2}}-{{{\boldsymbol{ \chi }} }_{1}}$ and $\delta \Phi ={{\Phi }_{2}}-{{\Phi }_{1}}+{\boldsymbol{ \chi }} _{1}^{{\rm T}}J{{{\boldsymbol{ \chi }} }_{2}}/2\hbar$, where the extra phase arises from the composition of displacement operators:

$$\begin{eqnarray}&&\hat{\mathcal{D}}(-{{{\boldsymbol{ \chi }} }_{1}})\hat{\mathcal{D}}({{{\boldsymbol{ \chi }} }_{2}})={{{\rm e}}^{\frac{{\rm i}}{2\hbar }{\boldsymbol{ \chi }} _{1}^{{\rm T}}J{{{\boldsymbol{ \chi }} }_{2}}}}\hat{\mathcal{D}}(\delta {\boldsymbol{ \chi }} ).\end{eqnarray} \tag{ 6 }$$

The result is analogous for the state $\mid {{\psi }_{{\rm II}}}(t)\rangle$ at the second exit port.

Zoom In Zoom Out Reset image size

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 section 3.2.

We will employ the term 'open interferometer' whenever $\delta {\boldsymbol{ \chi }} \ne {\bf 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 $\langle {{\psi }_{{\rm i}}}(t)\mid {{\psi }_{{\rm i}}}(t)\rangle$ and from equation (5) one can easily obtain

$$\begin{eqnarray}&&\frac{{{N}_{{\rm i}}}}{{{N}_{{\rm i}}}+{{N}_{{\rm II}}}}=\langle {{\psi }_{{\rm i}}}(t)\mid {{\psi }_{{\rm i}}}(t)\rangle =\frac{1}{2}\left( 1+C{\rm cos} \delta \phi \right),\end{eqnarray} \tag{ 7 }$$

with

$$\begin{eqnarray}&&C=\mid \langle {{\psi }_{{\rm c}}}(t)\mid \hat{\mathcal{D}}(\delta {\boldsymbol{ \chi }} )\mid {{\psi }_{{\rm c}}}(t)\rangle \mid \leqslant 1.\end{eqnarray} \tag{ 8 }$$

Therefore, $\delta {\boldsymbol{ \chi }} =0$ implies C = 1 and full oscillations, between 0 and N, in the total atom number at each exit port. In contrast, $\delta {\boldsymbol{ \chi }} \ne 0$ causes a contrast reduction³ with smaller oscillations around $N/2$ of the atom number at each port, as illustrated in figure 2. (Completely analogous results hold for the second exit port.) Note that in general the expectation value $\langle {{\psi }_{{\rm c}}}(t)\mid \hat{\mathcal{D}}(\delta {\boldsymbol{ \chi }} )\mid {{\psi }_{{\rm c}}}(t)\rangle$ is complex and its phase gives a contribution to the phase $\delta \phi$ in equation (7) in addition to the phase shift $\delta \Phi$ in equation (5).

Zoom In Zoom Out Reset image size

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 ${{{\boldsymbol{\mathcal{R}}} }_{1}}=0$ to ease the notation, the probability density associated with the state in equation (5) becomes

$$\begin{eqnarray}\begin{array}{rcl} \mid {{\psi }_{{\rm i}}}({\bf x},t){{\mid }^{2}} & = & \frac{1}{4}\left( \mid {{\psi }_{{\rm c}}}({\bf x},t){{\mid }^{2}}+\mid {{\psi }_{{\rm c}}}({\bf x}-\delta {\boldsymbol{\mathcal{R}}} ,t){{\mid }^{2}} \right) \\ {} & {} & +\frac{1}{2}\operatorname{Re}\left( {{{\rm e}}^{{\rm i}\delta \Phi ^{\prime} }}{{{\rm e}}^{\frac{{\rm i}}{\hbar }\delta {{{\boldsymbol{\mathcal{P}}} }^{{\rm T}}}{\bf x}}}\psi _{{\rm c}}^{*}({\bf x},t){{\psi }_{{\rm c}}}({\bf x}-\delta {\boldsymbol{\mathcal{R}}} ,t) \right), \\ \end{array}\end{eqnarray} \tag{ 9 }$$

where some spatially independent phase terms (i.e. independent of ${\bf x}$) have been absorbed in the phase $\delta \Phi ^{\prime}$ 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 equation (C.3), equation (9) becomes

$$\begin{eqnarray}&&\mid {{\psi }_{{\rm i}}}({\bf x},t){{\mid }^{2}}\approx \frac{1}{4}\mid 1+{{{\rm e}}^{{\rm i}\delta {{\Phi }^{\prime\prime }}}}{{{\rm e}}^{\frac{{\rm i}}{\hbar }{{\left( \delta {\boldsymbol{\mathcal{P}}} -\frac{m}{\Delta t}\delta {\boldsymbol{\mathcal{R}}} \right)}^{{\rm T}}}{\bf x}}}{{\mid }^{2}}\mid {{\psi }_{{\rm c}}}({\bf x},t){{\mid }^{2}},\end{eqnarray} \tag{ 10 }$$

where some spatially independent phase terms have again been absorbed in the phase $\delta \Phi ^{\prime\prime}$ and $\Delta 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 $\tilde{\psi }(m\;{\bf x}/\Delta t,{{t}_{0}})$ in equation (C.3). Moreover, we have neglected the shift by $\delta {\boldsymbol{\mathcal{R}}}$ of the envelope, which is a reasonable approximation whenever the envelope $\mid {{\psi }_{{\rm c}}}({\bf x},t)\mid$ does not vary too rapidly in space and its size is much larger than $\delta {\boldsymbol{\mathcal{R}}}$. (This approximation, which is often an excellent one in typical situations of interest, will be discussed in detail in sections 3.3 and 4.2.) Equation (10) reveals the existence of a fringe pattern in the density profile of each exit port with a fringe spacing ${{\lambda }_{{\rm fr}}}$ characterized by

$$\begin{eqnarray}&&\frac{2\pi }{{{\lambda }_{{\rm fr}}}}{\bf \hat{n}}=\frac{1}{\hbar }\left( \delta {\boldsymbol{\mathcal{P}}} -\frac{m}{\Delta t}\delta {\boldsymbol{\mathcal{R}}} \right).\end{eqnarray} \tag{ 11 }$$

Integrating equation (9) over space, one gets the analog of equations (7) and (8) in position representation with

$$\begin{eqnarray}&&C=\mid \int {{{\rm d}}^{3}}x{{{\rm e}}^{\frac{{\rm i}}{\hbar }\delta {{{\boldsymbol{\mathcal{P}}} }^{{\rm T}}}{\bf x}}}\psi _{{\rm c}}^{*}({\bf x}+\delta {\boldsymbol{\mathcal{R}}} /2,t){{\psi }_{{\rm c}}}({\bf x}-\delta {\boldsymbol{\mathcal{R}}} /2,t)\mid ,\end{eqnarray} \tag{ 12 }$$

where we have shifted the integration variable ${\bf x}$ by $\delta {\boldsymbol{\mathcal{R}}} /2$ to write the result in a more convenient form. Specializing to the state in equation (C.3) and making use of the same assumptions employed when deriving equation (10), which includes neglecting the shift of the envelopes, equation (12) becomes⁴

$$\begin{eqnarray}&&C\approx \mid \int {{{\rm d}}^{3}}x{\rm exp} \left[ \frac{{\rm i}}{\hbar }{{\left( \delta {\boldsymbol{\mathcal{P}}} -\frac{m}{\Delta t}\delta {\boldsymbol{\mathcal{R}}} \right)}^{{\rm T}}}{\bf x} \right]\mid {{\psi }_{{\rm c}}}({\bf x},t){{\mid }^{2}}\mid .\end{eqnarray} \tag{ 13 }$$

It is then easy to see that when the fringe spacing is much larger than the size of the envelope, as shown in the upper row of figure 3, the port goes from almost completely bright to almost completely dark as $\delta \phi$ changes, and the integrated particle number at each port oscillates from 0 to N (with $C\approx 1$). On the other hand, for fringe spacing much smaller than the size of the envelope, as shown in the lower row of figure 3, $C\ll 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.

Zoom In Zoom Out Reset image size

2.3.1. Asymmetric pulse timing

We start by considering the simple example of an asymmetric Mach–Zehnder interferometer (MZI), where the time separations between the mirror pulse and the beam splitter pulses (initial and final) differ by $\delta T$, as shown in figure 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

$$\begin{eqnarray}&&\delta {\boldsymbol{\mathcal{R}}} =-{{{\bf v}}_{{\rm rec}}}\delta T.\end{eqnarray} \tag{ 14 }$$

This gives rise to a fringe pattern in the density profile with fringe spacing

$$\begin{eqnarray}&&{{\lambda }_{{\rm fr}}}=\frac{2\pi \hbar \;\Delta t}{m\mid \delta {\boldsymbol{\mathcal{R}}} \mid },\end{eqnarray} \tag{ 15 }$$

as observed for instance in the experiments reported in [31, 42]. 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 ${{{\bf v}}_{{\rm rec}}}$); see section 5 below for a more detailed discussion.

Zoom In Zoom Out Reset image size

2.3.2. Gravity gradient

Let us consider now the effect of the gravity gradient when the dynamics is governed by a Hamiltonian like that in equation (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 equation (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:

$$\begin{eqnarray}\begin{array}{rcl} C & = & \mid \langle {{\psi }_{{\rm c}}}({{t}_{{\rm bs}}})\mid \hat{U}_{0}^{\dagger }(t,{{t}_{{\rm bs}}})\hat{\mathcal{D}}(\delta {\boldsymbol{ \chi }} ){{{\hat{U}}}_{0}}(t,{{t}_{{\rm bs}}})\mid {{\psi }_{{\rm c}}}({{t}_{{\rm bs}}})\rangle \mid \\ {} & = & \mid \langle {{\psi }_{{\rm c}}}({{t}_{{\rm bs}}})\mid \hat{\mathcal{D}}\left( \delta {\boldsymbol{ \chi }} ({{t}_{{\rm bs}}}) \right)\mid {{\psi }_{{\rm c}}}({{t}_{{\rm bs}}})\rangle \mid , \\ \end{array}\end{eqnarray} \tag{ 16 }$$

where we took into account that $\mid {{\psi }_{{\rm c}}}(t)\rangle ={{\hat{U}}_{0}}(t,{{t}_{{\rm bs}}})\mid {{\psi }_{{\rm c}}}({{t}_{{\rm bs}}})\rangle$ and made use in the second equality of equation (B.5) applied to this case together with equation (A.11). Therefore, unless stated otherwise, from now on we will always calculate the displacement $\delta {\boldsymbol{ \chi }}$ at time ${{t}_{{\rm 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 $\mid {{\psi }_{{\rm c}}}(t)\rangle ={{\hat{U}}_{0}}(t,{{t}_{0}})\mid {{\psi }_{{\rm c}}}({{t}_{0}})\rangle$ when evaluating equation (8). Using instead the full unitary evolution operator would simply add small corrections to the result for the contrast suppressed by additional powers of $\Gamma {{T}^{2}}$.

Particularizing the results of equations (A.12) and (A.13) to the case of a MZI with half-interferometer time T and $\delta T=0$, one obtains the following position and momentum displacements at time ${{t}_{{\rm bs}}}$ and to first order in $\Gamma {{T}^{2}}$:

$$\begin{eqnarray}&&\delta {\boldsymbol{\mathcal{R}}} =\left( \Gamma \;{{T}^{2}} \right){{{\bf v}}_{{\rm rec}}}T,\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\delta {\boldsymbol{\mathcal{P}}} =(\Gamma \;{{T}^{2}})m\;{{{\bf v}}_{{\rm rec}}}.\end{eqnarray} \tag{ 18 }$$

These results are valid for an arbitrary orientation, but in the next subsection and sections 3 and 4 we will concentrate on the case of aligned gravity gradients, i.e. the case in which the direction of the laser pulses and, hence, ${{{\bf v}}_{{\rm rec}}}$ coincide with the direction of one of the principal axes of the gravity gradient tensor Γ. This in turn implies $\delta {\boldsymbol{\mathcal{R}}} \parallel {{{\bf v}}_{{\rm rec}}}$ and $\delta {\boldsymbol{\mathcal{P}}} \parallel {{{\bf v}}_{{\rm 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 section 5.

Note that the change of velocity associated with equation (18) is typically rather small (of the order of ${{10}^{-4}}\;{{{\bf v}}_{{\rm rec}}}$ or less for half-interferometer times up to $10$ s) and leads to negligible changes of the resonance condition given the typical finite time duration of the laser pulses (i.e. corresponding to energy changes much smaller than those allowed by the time–energy uncertainty relation). On the other hand, it is interesting to point out that a vanishing final displacement in phase space, $\delta {\boldsymbol{ \chi }} =0$, could be achieved by suitably adjusting the momentum transfer of the mirror pulse (this would lead for instance to trajectories symmetric under time inversion around the time of the mirror pulse in the case depicted in figure 5). Nevertheless, for long interferometer times this would require adjusting the laser frequencies by hundreds of GHz, which is not so easy to implement. An alternative strategy which is much more flexible and easier to implement will be presented in the next subsection and studied throughout the rest of the paper.

Zoom In Zoom Out Reset image size

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 $\delta 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 $\Gamma \;{{T}^{2}}$ and $\delta T$:

$$\begin{eqnarray}&&\delta {\boldsymbol{\mathcal{R}}} =-{{{\bf v}}_{{\rm rec}}}\delta T+\left( \Gamma \;{{T}^{2}} \right){{{\bf v}}_{{\rm rec}}}T,\end{eqnarray} \tag{ 19 }$$

but leaves $\delta {\boldsymbol{\mathcal{P}}}$ invariant at that order. It is, therefore, not possible to eliminate completely the phase-space displacement $\delta {\boldsymbol{ \chi }}$. Nevertheless, the key insight is that by choosing appropriately $\delta T,$ one can get a nonvanishing value for $\delta {\boldsymbol{\mathcal{R}}}$ that compensates the effect due to the nonvanishing $\delta {\boldsymbol{\mathcal{P}}}$ and causes the fringe spacing to be much larger than the size of the envelope. Indeed, we see from equation (11) that $1/{{\lambda }_{{\rm fr}}}\to 0$ when

$$\begin{eqnarray}&&\delta {\boldsymbol{\mathcal{R}}} -(t-{{t}_{0}})\delta {\boldsymbol{\mathcal{P}}} /m=0,\end{eqnarray} \tag{ 20 }$$

and this can be achieved for aligned gravity gradients by taking

$$\begin{eqnarray}&&\delta T=-(T+{{T}_{0}})\left( {{\Gamma }_{\parallel }}\;{{T}^{2}} \right),\end{eqnarray} \tag{ 21 }$$

where ${{\Gamma }_{\parallel }}$ is the eigenvalue of the Γ tensor along the direction of ${{{\bf v}}_{{\rm rec}}}$ and we have considered a total time $t-{{t}_{0}}=2T+{{T}_{0}}$, with T being the half-interferometer time and T₀ the time from t₀ 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.

Let us consider a general pure Gaussian state as the initial state of the centered wave packet, $\mid {{\psi }_{{\rm c}}}({{t}_{0}})\rangle$. Its associated Wigner function, defined in general as

$$\begin{eqnarray}&&W({\bf x},{\bf p})=\int \frac{{{{\rm d}}^{3}}\Delta }{{{\left( 2\pi \hbar \right)}^{3}}}{{{\rm e}}^{{\rm i}{{{\bf p}}^{{\rm T}}}{\boldsymbol{ \Delta }} /\hbar }}\langle {\bf x}-{\boldsymbol{ \Delta }} /2\mid {{\psi }_{{\rm c}}}\rangle \langle {{\psi }_{{\rm c}}}\mid {\bf x}+{\boldsymbol{ \Delta }} /2\rangle ,\end{eqnarray} \tag{ 22 }$$

takes the form

$$\begin{eqnarray}&&W({\bf x},{\bf p};{{t}_{0}})={{\left( 2\pi \right)}^{-3}}{{\left( {\rm det} \;\Sigma \right)}^{-\frac{1}{2}}}{{{\rm e}}^{-\frac{1}{2}{{{\boldsymbol{ \xi }} }^{{\rm T}}}{{\Sigma }^{-1}}{\boldsymbol{ \xi }} }},\end{eqnarray} \tag{ 23 }$$

where ${\boldsymbol{ \xi }} ={{({\bf x},{\bf p})}^{{\rm T}}}$ and Σ is the phase-space covariance matrix, which is directly related to the Weyl-ordered two-point functions:

$$\begin{eqnarray}{{\Sigma }_{ij}}=\frac{1}{2}\langle {{\hat{\xi }}_{i}}\ {{\hat{\xi }}_{j}}+{{\hat{\xi }}_{j}}\ {{\hat{\xi }}_{i}}\rangle ={{(\begin{array}{ccccccccccccccc} {{\Sigma }_{xx}} & {{\Sigma }_{xp}} \\ \Sigma _{xp}^{T} & {{\Sigma }_{pp}} \\ \end{array})}_{ij}}.\end{eqnarray} \tag{ 24 }$$

${{\Sigma }_{xx}},\;{{\Sigma }_{xp}}\;{\rm and}\;{{\Sigma }_{pp}}$ are $3\;\times \;3$ matrices that can be regarded as blocks of the $6\;\times \;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 [43]:

$$\begin{eqnarray}\begin{array}{rcl} \langle {{\psi }_{{\rm c}}}(t)\mid \hat{\mathcal{D}}(\delta {\boldsymbol{ \chi }} )\mid {{\psi }_{{\rm c}}}(t)\rangle & = & \int {{{\rm d}}^{3}}x^{\prime} \int {{{\rm d}}^{3}}p^{\prime} W({\bf x}^{\prime} ,{\bf p}^{\prime} ;t){{{\rm e}}^{-\frac{{\rm i}}{\hbar }\delta {{{\boldsymbol{ \chi }} }^{{\rm T}}}J{\boldsymbol{ \xi }} ^{\prime} }} \\ {} & = & \int {{{\rm d}}^{3}}x\int {{{\rm d}}^{3}}pW({\bf x},{\bf p};{{t}_{0}}){{{\rm e}}^{-\frac{{\rm i}}{\hbar }\delta {\boldsymbol{ \chi }} _{0}^{{\rm T}}J{\boldsymbol{ \xi }} }}, \\ \end{array}\end{eqnarray} \tag{ 25 }$$

where $\delta {{{\boldsymbol{ \chi }} }_{0}}\equiv \mathcal{T}({{t}_{0}},t)\delta {\boldsymbol{ \chi }}$ with the transition matrix $\mathcal{T}({{t}_{2}},{{t}_{1}})$ defined in appendix A. In the second equality we have introduced the change of variables ${\boldsymbol{ \xi }} ^{\prime} =\mathcal{T}(t,{{t}_{0}}){\boldsymbol{ \xi }}$, used the invariance of the symplectic form ${{\mathcal{T}}^{{\rm T}}}(t,{{t}_{0}})J\mathcal{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 equations (25) and (23) we get the following result for the contrast defined in equation (8):

$$\begin{eqnarray}&&C={{{\rm e}}^{-\frac{1}{2{{\hbar }^{2}}}\delta {\boldsymbol{ \chi }} _{0}^{{\rm T}}{{J}^{{\rm T}}}\Sigma J\delta {{{\boldsymbol{ \chi }} }_{0}}}}.\end{eqnarray} \tag{ 26 }$$

Up to a factor $1/2{{\hbar }^{2}}$ the exponent can be rewritten as

$$\begin{eqnarray}&&-\delta {\boldsymbol{\mathcal{R}}} _{0}^{({\rm s}){\rm T}}{{\Sigma }_{pp}}\delta {\boldsymbol{\mathcal{R}}} _{0}^{({\rm s})}-\delta {\boldsymbol{\mathcal{P}}} _{0}^{{\rm T}}\left( {{\Sigma }_{xx}}-{{\Sigma }_{xp}}\Sigma _{pp}^{-1}\Sigma _{xp}^{{\rm T}} \right)\delta {{{\boldsymbol{\mathcal{P}}} }_{0}},\end{eqnarray} \tag{ 27 }$$

with $\delta {\boldsymbol{\mathcal{R}}} _{0}^{({\rm s})}=\delta {{{\boldsymbol{\mathcal{R}}} }_{0}}-\left( \Sigma _{pp}^{-1}\Sigma _{xp}^{T} \right)\delta {{{\boldsymbol{\mathcal{P}}} }_{0}}$. In order to write $\delta {{{\boldsymbol{ \chi }} }_{0}}$ in terms of $\delta {\boldsymbol{ \chi }}$ one needs to use the transition matrix $\mathcal{T}({{t}_{0}},t)$, which is exactly given by equation (A.8) for time-independent gravity gradients. However, as explained in the introduction and in section 2.3, 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 $\mid {{\psi }_{{\rm c}}}(t)\rangle$. 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

$$\begin{eqnarray}&&\delta {{{\boldsymbol{\mathcal{R}}} }_{0}}=\delta {\boldsymbol{\mathcal{R}}} -\frac{\delta {\boldsymbol{\mathcal{P}}} }{m}(t-{{t}_{0}}),\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&\delta {{{\boldsymbol{\mathcal{P}}} }_{0}}=\delta {\boldsymbol{\mathcal{P}}} .\end{eqnarray} \tag{ 29 }$$

Thus, from equations (27)–(29) one can see that for a given $\delta {\boldsymbol{\mathcal{P}}}$ the contrast is maximized when $\delta {\boldsymbol{\mathcal{R}}} _{0}^{({\rm s})}=0$, which implies

$$\begin{eqnarray}&&\delta {\boldsymbol{\mathcal{R}}} =\frac{\delta {\boldsymbol{\mathcal{P}}} }{m}(t-{{t}_{0}})+\Sigma _{pp}^{-1}\Sigma _{xp}^{{\rm T}}\delta {\boldsymbol{\mathcal{P}}} .\end{eqnarray} \tag{ 30 }$$

Note that since the exponent in equation (26) is manifestly negative for any $\delta {{{\boldsymbol{ \chi }} }_{0}}$ and one can always choose $\delta {{{\boldsymbol{\mathcal{R}}} }_{0}}$ so that $\delta {\boldsymbol{\mathcal{R}}} _{0}^{({\rm s})}=0$ for any given $\delta {{{\boldsymbol{\mathcal{P}}} }_{0}}$, the second term in expression (27) must be negative. The interpretation and importance of this term will be discussed in section 3.3. In this respect, it is instructive to consider the one-dimensional case, for which it becomes

$$\begin{eqnarray}&&-\frac{1}{2{{\hbar }^{2}}}\frac{{\rm det} \;\Sigma }{{{\Sigma }_{pp}}}\delta \mathcal{P}_{0}^{2}\leqslant -\frac{1}{8}\frac{\delta \mathcal{P}_{0}^{2}}{{{\Sigma }_{pp}}},\end{eqnarray} \tag{ 31 }$$

which follows from the inequality $({\rm det} \;\Sigma )\geqslant {{\hbar }^{2}}/4$. This inequality is the necessary and sufficient condition for the Gaussian phase-space distribution in equation (23) to be the Wigner function of a well-defined density matrix (with ${{\hat{\rho }}^{2}}\leqslant \hat{\rho }$); the equality holds for pure states and the strict inequality for mixed ones. Moreover, it coincides with the Schrödinger uncertainty relation [45, 46], a somewhat stronger version of Heisenbergʼs uncertainty relation.

There is an alternative expression for equation (25) as the overlap of two Wigner functions with a relative phase-space displacement $\delta {\boldsymbol{ \chi }}$ (valid only for pure states):

$$\begin{eqnarray}\begin{array}{rcl} C & = & \mid \langle {{\psi }_{{\rm c}}}(t)\mid \hat{\mathcal{D}}(\delta {\boldsymbol{ \chi }} )\mid {{\psi }_{{\rm c}}}(t)\rangle \mid ={{\left( {\rm Tr}\left[ \hat{\rho }\;{{{\hat{\rho }}}_{\delta {\boldsymbol{ \chi }} }} \right] \right)}^{\frac{1}{2}}} \\ {} & = & {{(2\pi \hbar )}^{\frac{3}{2}}}\mid \int {{{\rm d}}^{3}}x{{{\rm d}}^{3}}pW({\bf x},{\bf p};t)W\left( {\bf x}-\delta {\boldsymbol{\mathcal{R}}} ,{\bf p}-\delta {\boldsymbol{\mathcal{P}}} ;t \right){{\mid }^{\frac{1}{2}}}, \\ \end{array}\end{eqnarray} \tag{ 32 }$$

where the density matrices $\hat{\rho }$ and ${{\hat{\rho }}_{\delta {\boldsymbol{ \chi }} }}$ correspond respectively to the states $\mid {{\psi }_{{\rm c}}}(t)\rangle$ and $\hat{\mathcal{D}}(\delta {\boldsymbol{ \chi }} )\mid {{\psi }_{{\rm c}}}(t)\rangle$. Making use of equation (32) 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 figure 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 $\delta \mathcal{P}$ the overlap is maximized when $\delta \mathcal{R}$ is chosen so that the two displaced Wigner functions are aligned⁵ . Taking into account that $\Sigma _{pp}^{-1}\Sigma _{xp}^{{\rm T}}$ characterizes the 'tilt' in phase space of the initial Wigner function and considering its free evolution, one can easily obtain condition (30) for maximum contrast. It is also clear that substantial overlap will be achieved as long as $\delta \mathcal{P}\ll {{({{\Sigma }_{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.

Zoom In Zoom Out Reset image size

The result in equation (30) agrees at late times with the late-time result (20) derived in section 2. Furthermore, the second term on the right-hand side of equation (30), which is proportional to Σ_xp and suppressed by a factor $1/\Delta t$ compared to the first one, would have also been obtained in section 2 if we had considered spatially dependent phase contributions from $\tilde{\psi }(m\;{\bf x}/\Delta t,{{t}_{0}})$, which are always present for ${{\Sigma }_{xp}}\ne 0$.

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 as

$$\begin{eqnarray}&&\hat{\rho }({{t}_{0}})=\hat{\mathcal{D}}({{{\boldsymbol{ \chi }} }_{0}}){{\hat{\rho }}_{{\rm c}}}({{t}_{0}}){{\hat{\mathcal{D}}}^{\dagger }}({{{\boldsymbol{ \chi }} }_{0}}),\end{eqnarray} \tag{ 33 }$$

so that ${\rm Tr}\left[ {\boldsymbol{\hat{ \xi }}} \;{{{\hat{\rho }}}_{{\rm c}}}({{t}_{0}}) \right]=0$. Next, we express ${{\hat{\rho }}_{{\rm c}}}({{t}_{0}})$ as an incoherent mixture of pure states and write the density matrix for the initial state in the following way:

$$\begin{eqnarray}&&{{\hat{\rho }}_{{\rm c}}}({{t}_{0}})=\mathop{\sum }\limits_{j}{{p}_{j}}\mid {{\psi }^{(j)}}({{t}_{0}})\rangle \langle {{\psi }^{(j)}}({{t}_{0}})\mid \end{eqnarray} \tag{ 34 }$$

with $0\leqslant {{p}_{j}}\leqslant 1$ and ${{\sum }_{j}}{{p}_{j}}=1$. One can then apply the results obtained for pure states to each state $\mid {{\psi }^{(j)}}({{t}_{0}})\rangle$ 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 ${\boldsymbol{\hat{ \xi }}}$ 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 $\mid {{\psi }^{(j)}}({{t}_{0}})\rangle$ as centered states. With this choice the term $C{\rm cos} \delta \phi$ appearing in equation (7) corresponds for each pure state of the ensemble to the real part of

$$\begin{eqnarray}&&{{C}^{(j)}}{{{\rm e}}^{{\rm i}\delta {{\phi }^{(j)}}}}\equiv \langle {{\psi }^{(j)}}(t)\mid \hat{\mathcal{D}}(\delta {\boldsymbol{ \chi }} )\mid {{\psi }^{(j)}}(t)\rangle {{{\rm e}}^{{\rm i}\delta \Phi }},\end{eqnarray} \tag{ 35 }$$

where ${{C}^{(j)}}$ is defined as the modulus of the expression on the right-hand side (recall that $\delta {{\phi }^{(j)}}$ contains a possible additional phase from the expectation value of the displacement operator). The relative displacement $\delta {\boldsymbol{ \chi }}$ is the same for all members of the ensemble, and so is $\delta \Phi$, which depends on the initial displacement ${{{\boldsymbol{ \chi }} }_{0}}$ appearing in equation (33), 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:

$$\begin{eqnarray}&&C\;{{{\rm e}}^{{\rm i}\delta \phi }}=\mathop{\sum }\limits_{j}{{p}_{j}}{{C}^{(j)}}{{{\rm e}}^{{\rm i}\delta {{\phi }^{(j)}}}}\end{eqnarray} \tag{ 36 }$$

whose real part characterizes the oscillations of the total particle number at the exit port for the whole ensemble, as given by equation (7). In particular, for the contrast one gets

$$\begin{eqnarray}&&C=\mid \mathop{\sum }\limits_{j}{{p}_{j}}{{C}^{(j)}}{{{\rm e}}^{{\rm i}\delta {{\phi }^{(j)}}}}\mid \leqslant \mathop{\sum }\limits_{j}{{p}_{j}}{{C}^{(j)}},\end{eqnarray} \tag{ 37 }$$

where one has a strict inequality unless all the phases $\delta {{\phi }^{(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 $\delta \phi$, even when there is hardly any loss of contrast for each member of the ensemble.

It is clear from equation (37) that the result for pure states in equation (8) admits the following natural generalization to mixed states:

$$\begin{eqnarray}&&C=\mid {\rm Tr}\left[ \hat{\mathcal{D}}(\delta {\boldsymbol{ \chi }} ){{{\hat{\rho }}}_{{\rm c}}}(t) \right]\mid ,\end{eqnarray} \tag{ 38 }$$

which in turn implies that equation (25), 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 $\mid {{\psi }_{{\rm c}}}\rangle \langle {{\psi }_{{\rm c}}}\mid$ with ${{\hat{\rho }}_{{\rm c}}}$ in equation (22).] In contrast, equation (32), 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 ${\boldsymbol{\hat{ \xi }}}$ for the centered states, this would have led to a different initial displacement ${\boldsymbol{ \chi }} _{0}^{(j)}$ for each pure state $\mid {{\psi }^{(j)}}({{t}_{0}})\rangle$ and, hence, a different $\delta {{\Phi }^{(j)}}$ too. However, this would be exactly compensated by the expectation value of the displacement operator with respect to the new centered states $\mid \psi _{{\rm c}}^{(j)}(t)\rangle$ so that $\delta {{\phi }^{(j)}}$ would remain unchanged.

Finally, it should be emphasized that so far the whole discussion and conclusions about mixed states, including the validity of equation (25), were completely general and not just restricted to Gaussian states. When particularizing to Gaussian density matrices, one can further conclude that equation (26), derived for pure Gaussian states in the previous subsection, can be directly applied to mixed states as well.

In order to illustrate how effective our proposed mitigation strategy can be, we consider the effects of an aligned gravity gradient with ${{\Gamma }_{zz}}=3\;\times \;{{10}^{-6}}\;{{{\rm s}}^{-2}}$, comparable to that on the surface of Earth, and a pure Gaussian state for $^{87}{\rm Rb}$ atoms with ${{\Sigma }_{zz}}={{\left( 100\;\mu {\rm m} \right)}^{2}}$ and a momentum width along this direction corresponding to an effective temperature ${{T}_{{\rm eff}}}=1\;{\rm nK}$ (implying a velocity width of $0.3\;{\rm mm}\;{{{\rm s}}^{-1}}$). This temperature, which has already been successfully implemented in interferometry experiments employing the so-called delta-kick cooling technique [31], is much less demanding than the requirement of ${{T}_{{\rm eff}}}\approx 70\;{\rm 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 figure 7 one can see how the gravity gradient leads to a substantial loss of contrast for $T=5\;{\rm s}$ of half-interferometer time, whereas the mitigation strategy enables an extension beyond $T=20\;{\rm 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\;{\rm s}$. Note that determining the gravity gradient with a 10% accuracy is quite feasible since the effects of the local topography and nearby objects are typically of that order at most. Moreover, commercially available gradiometers can provide measurements of the gravity gradient with 0.1% accuracy, and the contribution of local masses can also be modeled with sufficient accuracy, for instance the payload in a space mission.

Zoom In Zoom Out Reset image size

In the previous example Σ_zp had a nonvanishing value which was determined (up to a sign) by the requirement that $({\rm det} \;\Sigma )={{\hbar }^{2}}/4$ for pure Gaussian states. As a second example we will consider a Gaussian mixed state under the same conditions but with ${{\Sigma }_{zp}}=0$ as well as ${{\Sigma }_{zz}}={{\left( 141\;\mu {\rm m} \right)}^{2}}$ and Σ_pp corresponding to ${{T}_{{\rm eff}}}=2\;{\rm 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 figure 8, the mitigation strategy is somewhat less effective in this case and the contrast is reduced almost completely for $T\gtrsim 10\;{\rm s}$. This has a simple intuitive explanation in terms of the dephasing between the different members of the ensemble of pure states discussed in section 3.2 and leading to the inequality in equation (37). 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 (except for their centers) with the pure state in our first example [47]. 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.

Zoom In Zoom Out Reset image size

We finish this section by discussing the role and interpretation of the term quadratic in $\delta {\boldsymbol{\mathcal{P}}}$ appearing in equation (27), 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 $\delta {\boldsymbol{\mathcal{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 $\delta \mathcal{P}/\left( 2\;\sqrt{{{\Sigma }_{pp}}} \right)$, and minus one half the square of this quantity (as in equation (31)) 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 equation (18) into equation (31), 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):

$$\begin{eqnarray}&&\frac{1}{8}\frac{\delta {{\mathcal{P}}^{2}}}{{{\Sigma }_{pp}}}=\frac{1}{8}{{\left[ \left( {{\Gamma }_{zz}}{{T}^{2}} \right)\frac{{{v}_{{\rm rec}}}}{{{\sigma }_{v}}} \right]}^{2}}.\end{eqnarray} \tag{ 39 }$$

Thus, for $T=10\;{\rm s}$ and typical values of the gravity gradient on Earthʼs surface this term is very small all the way down to ${{\sigma }_{v}}\sim {{10}^{-3}}\;{{v}_{{\rm rec}}}$, which corresponds to an extremely narrow velocity distribution (and ${{T}_{{\rm eff}}}\sim 1\;{\rm pK}$ for $^{87}{\rm Rb}$). On the other hand, as seen in equation (31), for mixed states the result in equation (39) is enhanced by a factor $4\;({\rm det} \;\Sigma )/{{\hbar }^{2}}\gt 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.

Throughout this section we will consider an expanding BEC evolving through an interferometer sequence after being released from a trap [31, 48], 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 (D.1) 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 section 2 by taking an initial state $\mid {{\psi }_{{\rm c}}}({{t}_{1}})\rangle$ 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:

$$\begin{eqnarray}&&{{\psi }_{{\rm c}}}({\bf x},t)\approx {{\left( N{\rm det} \;\Lambda (t) \right)}^{-\frac{1}{2}}}{{{\rm e}}^{{\rm i}\beta (t)}}{{{\rm e}}^{{\rm i}\frac{m}{2\hbar }{{{\bf x}}^{{\rm T}}}\dot{\Lambda }{{\Lambda }^{-1}}{\bf x}}}\;{{\psi }_{{\rm TF}}}({{\Lambda }^{-1}}{\bf x}),\end{eqnarray} \tag{ 40 }$$

where the different objects appearing in this expression are defined in equations (D.5), (D.8), (D.9) and (D.11). Furthermore, it also describes correctly the free evolution at arbitrarily late times, as shown in D.3. Therefore, one can directly employ equation (40) 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 $\mid {{\psi }_{{\rm TF}}}{{\mid }^{2}}$ is normalized to the number of atoms in the BEC, we have divided by $\sqrt{N}$ so that ${{\psi }_{{\rm c}}}$ agrees with the single-particle normalization employed in sections 2 and 3.

Substituting the scaling solution (40) into equation (9), one obtains the following result for the probability density, analogous to equation (10):

$$\begin{eqnarray}&&\mid {{\psi }_{{\rm i}}}({\bf x},t){{\mid }^{2}}\approx \frac{1}{4N{\rm det} \;\Lambda }\mid 1+{{{\rm e}}^{{\rm i}\delta \phi }}{{{\rm e}}^{\frac{{\rm i}}{\hbar }{{(\delta {\boldsymbol{\mathcal{P}}} -m\dot{\Lambda }{{\Lambda }^{-1}}\delta {\boldsymbol{\mathcal{R}}} )}^{{\rm T}}}{\bf x}}}{{\mid }^{2}}\mid {{\psi }_{{\rm TF}}}({{\Lambda }^{-1}}{\bf x}){{\mid }^{2}},\end{eqnarray} \tag{ 41 }$$

where we have assumed that the size of the envelope (the rescaled Thomas–Fermi radius) is much larger than $\delta {\boldsymbol{\mathcal{R}}}$ and have neglected the relative shift between the envelopes. The mitigation strategy is completely analogous to that already discussed in section 2.4 for the free-particle case. From equation (41) we see that the limit of large fringe spacing (larger than the size of the envelope) corresponds to the condition

$$\begin{eqnarray}&&\delta {\boldsymbol{\mathcal{P}}} -m\dot{\Lambda }{{\Lambda }^{-1}}\delta {\boldsymbol{\mathcal{R}}} =0.\end{eqnarray} \tag{ 42 }$$

This can be achieved with a small shift $\delta T$ in the timing of the last beam splitter pulse, which can be easily determined by substituting equations (17) and (18) into equation (42) and gives a result similar to equation (21).

At late times ${{\Lambda }_{ii}}\approx {{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:

$$\begin{eqnarray}&&{{\left( \dot{\Lambda }{{\Lambda }^{-1}} \right)}_{ii}}=\frac{1}{\Delta t}\left[ 1-\frac{{{c}_{i}}}{{{b}_{i}}}\frac{1}{\Delta t}+O\left( 1/\Delta {{t}^{2}} \right) \right].\end{eqnarray} \tag{ 43 }$$

Thus, to lowest order $\dot{\Lambda }{{\Lambda }^{-1}}\approx (1/\Delta t){\mathbb{1}}$ and the late-time limit of equation (41) agrees with the result in equation (10). Specifically, the terms of order 1 and $1/\Delta t$ in equation (42) and the exponent of equation (41) coincide with those in equation (11) of section 2.2. Moreover, if one takes an initial time t₀ when the expansion is already in the linear regime, the term of order $1/\Delta {{t}^{2}}$ on the right-hand side of equation (43) gives rise to a contribution in equation (42) analogous to the term proportional to $\Sigma _{pp}^{-1}\Sigma _{xp}^{{\rm T}}$ in equation (30) for the Gaussian case. This term is a consequence of having $\dot{\Lambda }({{t}_{0}})\ne 0$ and would have also been obtained in section 2.2 if the existence of a spatially dependent imaginary part of $\tilde{\psi }(m\;{\bf x}/\Delta t,{{t}_{0}})$, which it implies, had been taken into account there.

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 (40) into equation (12). If one neglects $\delta {\boldsymbol{\mathcal{R}}}$ compared to the size of the envelope, which is a very good approximation in most cases of practical interest (as illustrated by figure 9 and further discussed below), one gets the following result, analogous to equation (13):

$$\begin{eqnarray}&&C\approx \frac{1}{N}\int {{{\rm d}}^{3}}\zeta {\rm cos} \left[ \left( \delta {{{\boldsymbol{\mathcal{P}}} }^{{\rm T}}}\Lambda -m\;\delta {{{\boldsymbol{\mathcal{R}}} }^{{\rm T}}}\dot{\Lambda } \right){\boldsymbol{ \zeta }} /\hbar \right]\psi _{{\rm TF}}^{2}({\boldsymbol{ \zeta }} ),\end{eqnarray} \tag{ 44 }$$

where ${\boldsymbol{ \zeta }}$ are the rescaled coordinates defined in equation (D.3) and it has been taken into account that ${{\psi }_{{\rm TF}}}({\boldsymbol{ \zeta }} )$, given by equation (D.11), is a real and even function. We have also made use of the equality ${{\Lambda }^{{\rm T}}}\dot{\Lambda }={{\dot{\Lambda }}^{{\rm T}}}\Lambda$, proved at the end of D.1.

Zoom In Zoom Out Reset image size

For ${{{\bf v}}_{{\rm rec}}}$ aligned with a principal axis of both the trapping potential and the gravity gradient tensor Γ, one can choose a Cartesian coordinate system where $\Lambda ={\rm diag}\left( {{\Lambda }_{\parallel }},\Lambda _{\bot }^{(1)},\Lambda _{\bot }^{(2)} \right)$ with the longitudinal direction ∥ corresponding to the direction of ${{{\bf v}}_{{\rm rec}}}$. Integrating over the transverse directions, equation (44) becomes

$$\begin{eqnarray}&&C\approx \frac{15}{16}\mid \int _{-R_{{\rm TF}}^{\parallel }}^{R_{{\rm TF}}^{\parallel }}\frac{{\rm d}{{\zeta }_{\parallel }}}{R_{{\rm TF}}^{\parallel }}{\rm cos} \left( \frac{2\pi }{{{{\bar{\lambda }}}_{{\rm fr}}}}{{\zeta }_{\parallel }} \right){{\left[ 1-{{\left( \frac{{{\zeta }_{\parallel }}}{R_{{\rm TF}}^{\parallel }} \right)}^{2}} \right]}^{2}}\mid ,\end{eqnarray} \tag{ 45 }$$

where $R_{{\rm TF}}^{\parallel }$ is the Thomas–Fermi radius along the longitudinal direction and $2\pi /{{\bar{\lambda }}_{{\rm fr}}}=\left( \delta \mathcal{P}\;{{\Lambda }_{\parallel }}-m\;\delta \mathcal{R}\;{{{\dot{\Lambda }}}_{\parallel }} \right)/\hbar$ with ${{\bar{\lambda }}_{{\rm fr}}}$ being the fringe spacing in the rescaled coordinates. The result in equation (45) can be easily generalized to nonaligned gravity gradients and traps. First, one determines ${{\lambda }_{{\rm fr}}}$ and the longitudinal direction from the modulus and direction of the vector contracted with ${\boldsymbol{ \zeta }}$ within the argument of the cosine in equation (44) (note that the matrix $\Lambda (t)$ is not necessarily diagonal anymore). Next, one integrates the two transverse directions, which results in the 1D marginal distribution of the density distribution $\psi _{{\rm TF}}^{2}({\boldsymbol{ \zeta }} )$. As shown in [49], the form of this marginal distribution is the same as that appearing in equation (45) even for nonaligned traps: one simply needs to replace the square of the Thomas–Fermi radius ${{(R_{{\rm TF}}^{\parallel })}^{2}}$ with the ${{\Sigma }_{\parallel \parallel }}$ component of the matrix Σ characterizing the ground state in equation (D.11). Thus, equation (45) is also applicable to the general case with this simple replacement.

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 (40).

Similarly to what was done in section 3.3 for Gaussian states, we will consider a gravity gradient with ${{\Gamma }_{\parallel }}=3\;\times \;{{10}^{-6}}\;{{{\rm s}}^{-2}}$ and aligned with a principal axis of the trap potential and with v_rec. In addition, we will characterize our BEC of $^{87}{\rm Rb}$ atoms at some time t₀ when most of the repulsive interaction energy has been converted into kinetic energy and the expansion dynamics has entered the linear regime, with $\Lambda (t)=C+B\;(t-{{t}_{0}})$. Furthermore, since the result for the contrast in equation (45) only depends on the product $R_{{\rm TF}}^{\parallel }{{\Lambda }_{\parallel }}(t)$ (and its time derivative), one can entirely parametrize the system in terms of $R_{{\rm TF}}^{\parallel }{{\Lambda }_{\parallel }}({{t}_{0}})$ and $R_{{\rm TF}}^{\parallel }{{\dot{\Lambda }}_{\parallel }}({{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 figure 9 (with and without the mitigation strategy) for an initial half-width $R_{{\rm TF}}^{\parallel }{{\Lambda }_{\parallel }}({{t}_{0}})=60\;\mu {\rm m}$ and two different widths of the velocity distribution: $0.3\;{\rm mm}\;{{{\rm s}}^{-1}}$ and $0.08\;{\rm mm}\;{{{\rm s}}^{-1}}$ (with associated effective temperatures of $1\;{\rm nK}$ and $70\;{\rm pK}$).

From equation (44) and as seen in figure 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 (40), 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 (40) into equation (12) without neglecting the shift between the two envelopes:

$$\begin{eqnarray}\begin{array}{rcl} C & \approx & \frac{1}{N}\mid \int {{{\rm d}}^{3}}\zeta {\rm cos} \left[ \left( \delta {{{\boldsymbol{\mathcal{P}}} }^{{\rm T}}}\Lambda -m\;\delta {{{\boldsymbol{\mathcal{R}}} }^{{\rm T}}}\dot{\Lambda } \right){\boldsymbol{ \zeta }} \;/\hbar \right] \\ {} & {} & \times {{\psi }_{{\rm TF}}}\left( {\boldsymbol{ \zeta }} +{{\Lambda }^{-1}}\delta {\boldsymbol{\mathcal{R}}} /2 \right){{\psi }_{{\rm TF}}}\left( {\boldsymbol{ \zeta }} -{{\Lambda }^{-1}}\delta {\boldsymbol{\mathcal{R}}} /2 \right)\mid , \\ \end{array}\end{eqnarray} \tag{ 46 }$$

which is an exact result for the centered wave function ${{\psi }_{{\rm c}}}({\bf x},t)$ given by equation (40). In deriving equation (46) we have taken again into account that ${{\psi }_{{\rm TF}}}({\boldsymbol{ \zeta }} )$ is a real and even function.

For aligned gravity gradients and traps one can proceed similarly to what was done above to obtain equation (45). While the integrals with respect to the transverse coordinates in equation (46) can also be performed analytically, in this case there is no simple analytical result for the integral with respect to the longitudinal coordinate ${{\zeta }_{\parallel }}$ due to the more complicated functional dependence of ${{\psi }_{{\rm TF}}}\left( {\boldsymbol{ \zeta }} +{{\Lambda }^{-1}}\delta {\boldsymbol{\mathcal{R}}} /2 \right){{\psi }_{{\rm TF}}}\left( {\boldsymbol{ \zeta }} -{{\Lambda }^{-1}}\delta {\boldsymbol{\mathcal{R}}} /2 \right)$ compared to $\psi _{{\rm TF}}^{2}({\boldsymbol{ \zeta }} )$. This last integral can be evaluated numerically, which was done to obtain the results represented as dots in figure 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 section 3.3. 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 figure 9 the effect is very small and can only start to be appreciated for very long half-interferometer times. Moreover, it is slightly more significant for lower ${{T}_{{\rm eff}}}$ because the size of the envelopes at late times is smaller in that case.

With respect to the Hamiltonian ${{\hat{H}}_{0}}$ the phase-space vector ${\boldsymbol{\hat{ \xi }}} ={{({\bf \hat{x}},{\bf \hat{p}})}^{{\rm T}}}$, consisting of the position and momentum operators, exhibits the following relationship between the Heisenberg and the Schrödinger pictures: ${{{\boldsymbol{\hat{ \xi }}} }_{{\rm h}}}=\hat{U}_{0}^{\dagger }(t,{{t}_{0}}){\boldsymbol{\hat{ \xi }}} {{\hat{U}}_{0}}(t,{{t}_{0}})=\mathcal{T}(t,{{t}_{0}}){\boldsymbol{\hat{ \xi }}}$, where $\mathcal{T}(t,{{t}_{0}})$ coincides with the transition matrix associated to the homogeneous part of equation (A.1), which is a consequence of linearity and the fact that the Heisenberg equations have the same form as the classical equations of motion. Such a relation can be inverted to give ${{\hat{U}}_{0}}(t,{{t}_{0}}){\boldsymbol{\hat{ \xi }}} \hat{U}_{0}^{\dagger }(t,{{t}_{0}})={{\mathcal{T}}^{-1}}(t,{{t}_{0}}){\boldsymbol{\hat{ \xi }}}$. This together with equation (A.7) implies the following result for the displacement operator:

$$\begin{eqnarray}&&{{\hat{U}}_{0}}(t,{{t}_{0}})\hat{\mathcal{D}}\left( {{{\boldsymbol{ \chi }} }_{0}} \right)\hat{U}_{0}^{\dagger }(t,{{t}_{0}})=\hat{\mathcal{D}}\left( {{{\boldsymbol{ \chi }} }_{{\rm h}}}(t) \right),\end{eqnarray} \tag{ B.5 }$$

where ${{{\boldsymbol{ \chi }} }_{{\rm h}}}(t)=\mathcal{T}(t,{{t}_{0}}){{{\boldsymbol{ \chi }} }_{0}}$ corresponds to the classical phase-space trajectory associated with the Hamiltonian H₀ and with initial condition ${{{\boldsymbol{ \chi }} }_{{\rm h}}}({{t}_{0}})={{{\boldsymbol{ \chi }} }_{0}}$.

Taking all that into account and inserting the resolution of the identity $I=\hat{U}_{0}^{\dagger }(t,{{t}_{0}}){{\hat{U}}_{0}}(t,{{t}_{0}})$ after the displacement operator on the right-hand side of equation (B.1), the evolution of that state generated by the purely quadratic part of the Hamiltonian can be written in the following suggestive way:

$$\begin{eqnarray}&&\mid \psi (t)\rangle ={{\hat{U}}_{0}}(t,{{t}_{0}})\mid \psi ({{t}_{0}})\rangle =\hat{\mathcal{D}}\left( {{{\boldsymbol{ \chi }} }_{{\rm h}}}(t) \right)\mid {{\psi }_{{\rm c}}}(t)\rangle ,\end{eqnarray} \tag{ B.6 }$$

with $\mid {{\psi }_{{\rm c}}}(t)\rangle ={{\hat{U}}_{0}}(t,{{t}_{0}})\mid {{\psi }_{{\rm c}}}({{t}_{0}})\rangle$ and

$$\begin{eqnarray}&&{{{\boldsymbol{ \chi }} }_{{\rm h}}}(t)=\mathcal{T}(t,{{t}_{0}}){{{\boldsymbol{ \chi }} }_{0}}.\end{eqnarray} \tag{ B.7 }$$

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 nth pulse) by a phase factor ${{{\rm e}}^{{\rm i}{{\varepsilon }_{n}}{{\varphi }_{n}}}}\;{{{\rm e}}^{{\rm i}\mid {{\varepsilon }_{n}}\mid \frac{\pi }{2}}}$ times a displacement operator ${{{\rm e}}^{{\rm i}{{\varepsilon }_{n}}{\bf k}_{n}^{{\rm T}}{\bf \hat{x}}}}=\hat{\mathcal{D}}\left( {{{\boldsymbol{\mathcal{F}}} }_{n}} \right)$ acting on every wave packet. The integer _n, the laser phase φ_n, the momentum transfer $\hbar {{{\bf k}}_{n}}$ and the associated phase-space displacement ${{{\boldsymbol{\mathcal{F}}} }_{n}}$ (defined in equation (A.4)) 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 form $\hat{\mathcal{D}}\left( {\boldsymbol{ \chi }} (t_{n}^{-}) \right)\mid {{\psi }_{{\rm c}}}({{t}_{n}})\rangle$ right before, the action of a laser pulse at time t_n generates the following state:

$$\begin{eqnarray}\begin{array}{rcl} \mid \psi ({{t}_{n}})\rangle & = & {{{\rm e}}^{{\rm i}{{\varepsilon }_{n}}{{\varphi }_{n}}}}\hat{\mathcal{D}}\left( {{{\boldsymbol{\mathcal{F}}} }_{n}} \right)\hat{\mathcal{D}}\left( {\boldsymbol{ \chi }} (t_{n}^{-}) \right)\mid {{\psi }_{{\rm c}}}({{t}_{n}})\rangle \\ {} & = & {{{\rm e}}^{{\rm i}{{\varepsilon }_{n}}{{\varphi }_{n}}}}{{{\rm e}}^{-\frac{{\rm i}}{2\hbar }{\boldsymbol{\mathcal{F}}} _{n}^{{\rm T}}J\,{\boldsymbol{ \chi }} (t_{n}^{-})}}\hat{\mathcal{D}}\left( {\boldsymbol{ \chi }} (t_{n}^{-})+{{{\boldsymbol{\mathcal{F}}} }_{n}} \right)\mid {{\psi }_{{\rm c}}}({{t}_{n}})\rangle , \\ \end{array}\end{eqnarray} \tag{ B.8 }$$

where we used the composition formula in equation (6). Evolution between pulses proceeds as described in the previous subsection, so that the state at some time $t\gt {{t}_{n}}$ before the action of any further pulse is given by

$$\begin{eqnarray}&&\mid \psi (t)\rangle ={{{\rm e}}^{{\rm i}\Phi (t)}}\;\hat{\mathcal{D}}\left( {\boldsymbol{ \chi }} (t) \right)\mid {{\psi }_{{\rm c}}}(t)\rangle ,\end{eqnarray} \tag{ B.9 }$$

where the displacement ${\boldsymbol{ \chi }} (t)$ corresponds to the classical trajectory with initial conditions ${\boldsymbol{ \chi }} ({{t}_{0}})={{{\boldsymbol{ \chi }} }_{0}}$ and including the kicks from the laser pulses:

$$\begin{eqnarray}\begin{array}{rcl} {\boldsymbol{ \chi }} (t) & = & \mathcal{T}(t,{{t}_{n}}){\boldsymbol{ \chi }} (t_{n}^{-})+\mathcal{T}(t,{{t}_{n}}){{{\boldsymbol{\mathcal{F}}} }_{n}} \\ {} & = & {{{\boldsymbol{ \chi }} }_{{\rm h}}}(t)+\left( {{\mathcal{T}}_{{\rm ret}}}\cdot {{{\boldsymbol{\mathcal{F}}} }_{{\rm lp}}} \right)(t), \\ \end{array}\end{eqnarray} \tag{ B.10 }$$

with ${{{\boldsymbol{\mathcal{F}}} }_{{\rm lp}}}$ defined in equation (A.4). The phase $\Phi (t)$ in equation (B.9) is in turn given by

$$\begin{eqnarray}&&\Phi (t)=\varphi -\frac{1}{2\hbar }\mathop{\sum }\limits_{n}^{j=1}{\boldsymbol{\mathcal{F}}} _{j}^{{\rm T}}J\;{\boldsymbol{ \chi }} ({{t}_{j}}),\end{eqnarray} \tag{ B.11 }$$

with $\varphi =\sum _{j=1}^{n}({{\varepsilon }_{j}}\;{{\varphi }_{j}}+\mid {{\varepsilon }_{j}}\mid \pi /2)$. In equations (B.10)–(B.11) we have already provided the general result for n pulses. It can be easily proved by induction. First, one uses equation (B.8) to prove the result for n = 1. Next, assuming that the result holds for $n-1$ pulses, one can use equation (B.8) 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 $\delta t\to 0$ of r pulses with ${{{\boldsymbol{\mathcal{F}}} }_{j}}=\delta t{\boldsymbol{\mathcal{G}}} ({{t}_{j}})$, where ${\boldsymbol{\mathcal{G}}} (t)$ is defined in equation (A.3), ${{t}_{j}}={{t}_{0}}+j\;\delta t$ with j = 1,..., r and $r\;\delta t=t-{{t}_{0}}$. Equation (B.10) is then generalized to equation (A.5), and (B.11) becomes

$$\begin{eqnarray}&&\Phi (t)=\varphi -\frac{1}{2\hbar }\int _{{{t}_{0}}}^{t}{\rm d}t^{\prime} \left( \left[ {\boldsymbol{\mathcal{F}}} _{{\rm lp}}^{{\rm T}}(t^{\prime} )+{{{\boldsymbol{\mathcal{G}}} }^{{\rm T}}}(t^{\prime} ) \right]J\;{\boldsymbol{ \chi }} (t^{\prime} )+2{{V}_{0}}(t^{\prime} ) \right),\end{eqnarray} \tag{ B.12 }$$

where we have also included the contribution from the spatially independent potential ${{V}_{0}}(t)$.

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 figure 1, the superposition is given by equation (5). The displacement $\delta {\boldsymbol{ \chi }}$ can be calculated as explained in appendix A. Here we present a compact but complete derivation of the result for the phase shift $\delta \Phi$.

First one takes the expression for $\delta \Phi$ provided after equation (5) and substitutes Φ₁ and Φ₂ by the result in equation (B.12). (For the moment we do not consider the contributions from the purely linear term of the potential and take ${\boldsymbol{\mathcal{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 ${\boldsymbol{\bar{ \chi }}} =({{{\boldsymbol{ \chi }} }_{1}}+{{{\boldsymbol{ \chi }} }_{2}})/2$ and $\delta {\boldsymbol{ \chi }} ={{{\boldsymbol{ \chi }} }_{2}}-{{{\boldsymbol{ \chi }} }_{1}}$. One also proceeds in the same way with the sets of pulse displacements ${\boldsymbol{\mathcal{F}}} _{j}^{(1)}$ and ${\boldsymbol{\mathcal{F}}} _{j}^{(2)}$ associated with the two classical trajectories as well as the phases φ₁ and φ₂. Together with several additional steps described in more detail below, this leads to the following intermediate and final result:

$$\begin{eqnarray}\begin{array}{rcl} \delta \Phi (t) & = & {{\Phi }_{2}}-{{\Phi }_{1}}+\frac{1}{2\hbar }\;{\boldsymbol{ \chi }} _{1}^{{\rm T}}J{{{\boldsymbol{ \chi }} }_{2}} \\ {} & = & \delta \varphi -\frac{1}{2\hbar }\int _{{{t}_{0}}}^{t}{\rm d}t^{\prime} \left[ \delta {\boldsymbol{\mathcal{F}}} _{{\rm lp}}^{{\rm T}}(t^{\prime} )J{\boldsymbol{\bar{ \chi }}} (t^{\prime} )+{\boldsymbol {\mathcal{\bar{F}}}} _{{\rm lp}}^{{\rm T}}(t^{\prime} )J\;\delta {\boldsymbol{ \chi }} (t^{\prime} )+2\;\delta {{V}_{0}}(t^{\prime} ) \right] \\ {} & {} & -\frac{1}{2\hbar }\;\delta {{{\boldsymbol{ \chi }} }^{{\rm T}}}(t)J{\boldsymbol{\bar{ \chi }}} (t) \\ {} & = & \delta \varphi -\frac{1}{2\hbar }\int _{{{t}_{0}}}^{t}{\rm d}t^{\prime} \left[ \delta {\boldsymbol{\mathcal{F}}} _{{\rm lp}}^{{\rm T}}(t^{\prime} )J{\boldsymbol{\bar{ \chi }}} (t^{\prime} )-\delta {{{\boldsymbol{ \chi }} }^{{\rm T}}}(t^{\prime} )J\;{{{{\boldsymbol {\mathcal{\bar{F}}}} }}_{{\rm lp}}}(t^{\prime} )+\delta {{{{\boldsymbol{\dot{ \chi }}} }}^{{\rm T}}}(t^{\prime} )J{\boldsymbol{\bar{ \chi }}} (t^{\prime} ) \right. \\ {} & {} & +\left. \delta {{{\boldsymbol{ \chi }} }^{{\rm T}}}(t^{\prime} )J\mathop{{{\boldsymbol{\bar{ \chi }}} }}\limits^{.}(t^{\prime} )+2\;\delta {{V}_{0}}(t^{\prime} ) \right] \\ {} & = & \delta \varphi -\frac{1}{\hbar }\int _{{{t}_{0}}}^{t}{\rm d}t^{\prime} \left[ \delta {\boldsymbol{\mathcal{F}}} _{{\rm lp}}^{{\rm T}}(t^{\prime} )J\;{\boldsymbol{\bar{ \chi }}} (t^{\prime} )+\delta {{V}_{0}}(t^{\prime} ) \right]. \\ \end{array}\end{eqnarray} \tag{ B.13 }$$

In writing the third equality, 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}^{{\rm T}}}=-J$. We also wrote the last term of the second equality as the time integral from t₀ to t of its time derivative, taking into account that $\delta {\boldsymbol{ \chi }} ({{t}_{0}})=0$ for t₀ 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 $\delta {\boldsymbol{\dot{ \chi }}} =\mathcal{H}\delta {\boldsymbol{ \chi }} +\delta {{{\boldsymbol{\mathcal{F}}} }_{{\rm lp}}}$ and $\mathop{{{\boldsymbol{\bar{ \chi }}} }}\limits^{.}=\mathcal{H}{\boldsymbol{\bar{ \chi }}} +{{{\boldsymbol {\mathcal{\bar{F}}}} }_{{\rm lp}}}$ that follow from equation (A.1). Second, we used the identity ${{\mathcal{H}}^{{\rm T}}}J=-J\mathcal{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 equation (A.5), including ${\boldsymbol{\mathcal{G}}}$, with ${{{\boldsymbol {\mathcal{\bar{F}}}} }_{{\rm lp}}}$ in place of ${{{\boldsymbol{\mathcal{F}}} }_{{\rm lp}}}$. Equation (B.13) agrees then with the general result obtained by Antoine and Bordé in [60, 61] (see also [8, 50–53] for earlier results including the effects of gravity gradients). Furthermore, one can easily generalize the result to the case with different values of ${\boldsymbol{\mathcal{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 caused, for example, by time-dependent magnetic field gradients in Raman-based interferometers due the quadratic Zeeman effect; furthermore, time-dependent homogeneous magnetic fields give rise in that case to nontrivial effects through $\delta {{V}_{0}}(t^{\prime} )$.) Considering branch-dependent ${\boldsymbol{\mathcal{G}}}$ and V₀ 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, as discussed in detail in [37]. The new result can be written in a compact form as

$$\begin{eqnarray}&&\delta \Phi (t)=\delta \varphi -\frac{1}{\hbar }\int _{{{t}_{0}}}^{t}{\rm d}t^{\prime} \left( \left[ \delta {\boldsymbol{\mathcal{F}}} _{{\rm lp}}^{{\rm T}}(t^{\prime} )+\delta {{{\boldsymbol{\mathcal{G}}} }^{{\rm T}}}(t^{\prime} ) \right]J\;{\boldsymbol{\bar{ \chi }}} (t^{\prime} )+\delta {{V}_{0}}(t^{\prime} ) \right),\end{eqnarray} \tag{ B.14 }$$

where ${\boldsymbol{\bar{ \chi }}} (t^{\prime} )$ is given by equation (A.5) with ${{{\boldsymbol {\mathcal{\bar{F}}}} }_{{\rm lp}}}$ and ${\boldsymbol {\mathcal{\bar{G}}}}$.

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, equation (B.14) becomes

$$\begin{eqnarray}\begin{array}{rcl} \delta \Phi (t) & = & \delta \varphi -\frac{1}{\hbar }\int _{{{t}_{0}}}^{t}{\rm d}t^{\prime} \;\delta {{V}_{0}}(t^{\prime} )-\frac{1}{\hbar }\;\delta {{{\boldsymbol{ \chi }} }^{{\rm T}}}(t)J\,\mathcal{T}(t,{{t}_{0}}){{{\boldsymbol{ \chi }} }_{0}} \\ {} & {} & -\frac{1}{\hbar }\int _{{{t}_{0}}}^{t}{\rm d}t^{\prime} \int _{{{t}_{0}}}^{t^{\prime} }{\rm d}{{t}^{\prime\prime }}\left[ \delta {\boldsymbol{\mathcal{F}}} _{{\rm lp}}^{{\rm T}}(t^{\prime} )+\delta {{{\boldsymbol{\mathcal{G}}} }^{{\rm T}}}(t^{\prime} ) \right]J \\ {} & {} & \times \mathcal{T}(t^{\prime} ,{{t}^{\prime\prime }})\left[ {{{{\boldsymbol {\mathcal{\bar{F}}}} }}_{{\rm lp}}}({{t}^{\prime\prime }})+{\boldsymbol {\mathcal{\bar{G}}}} ({{t}^{\prime\prime }}) \right], \\ \end{array}\end{eqnarray} \tag{ B.15 }$$

where we have used equation (A.5) for ${\boldsymbol{\bar{ \chi }}} (t^{\prime} )$, the identity $J\;\mathcal{T}(t^{\prime} ,{{t}_{0}})={{\mathcal{T}}^{{\rm T}}}(t,t^{\prime} )J\;\mathcal{T}(t,{{t}_{0}})$, which follows from equation (A.7), and the expression for $\delta {{{\boldsymbol{ \chi }} }^{{\rm T}}}(t)$ in terms of $\delta {\boldsymbol{\mathcal{F}}} _{{\rm lp}}^{{\rm T}}$ and $\delta {{{\boldsymbol{\mathcal{G}}} }^{{\rm T}}}$ obtained from equation (A.5).

It should be noted that although we have focused on the Hamiltonian (B.2), the derivation in this appendix can be directly applied to a general quadratic Hamiltonian (including general linear terms) as well. In particular, the results can be straightforwardly applied to calculations in rotating frames, where new quadratic terms associated with the Coriolis and centrifugal forces arise.

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 [64]. The approach can also be implemented directly at the level of the Gross–Pitaevskii equation [63, 65]. One starts by introducing the linearly rescaled coordinates

$$\begin{eqnarray}&&{\boldsymbol{ \zeta }} ={{\Lambda }^{-1}}(t){\bf x},\end{eqnarray} \tag{ D.3 }$$

where $\Lambda (t)$ is a time-dependent matrix to be determined below. Next, one considers the following rescaling of the wave function:

$$\begin{eqnarray}&&\psi ({\bf x},t)={{\left( {\rm det} \;\Lambda (t) \right)}^{-\frac{1}{2}}}{{{\rm e}}^{{\rm i}\beta (t)}}{{{\rm e}}^{{\rm i}\Phi ({\boldsymbol{ \zeta }} ,t)}}{{\psi }_{\Lambda }}({\boldsymbol{ \zeta }} ,t),\end{eqnarray} \tag{ D.4 }$$

where the prefactor with the square root of the determinant guarantees a simple normalization condition for the new wave function ${{\psi }_{\Lambda }}({\boldsymbol{ \zeta }} ,t)$ in terms of the rescaled coordinates ${\boldsymbol{ \zeta }}$, namely $\int {{{\rm d}}^{3}}\zeta \mid {{\psi }_{\Lambda }}({\boldsymbol{ \zeta }} ,t){{\mid }^{2}}=N$. The spatially dependent and independent phases $\Phi ({\boldsymbol{ \zeta }} ,t)$ and $\beta (t)$ are given respectively by

$$\begin{eqnarray}&&\beta (t)=-\frac{\mu }{\hbar }\int _{{{t}_{0}}}^{t}{\rm d}t^{\prime} {{\left( {\rm det} \;\Lambda (t^{\prime} ) \right)}^{-1}},\end{eqnarray} \tag{ D.5 }$$

$$\begin{eqnarray}&&\Phi ({\boldsymbol{ \zeta }} ,t)=\frac{m}{2\hbar }{{{\boldsymbol{ \zeta }} }^{{\rm T}}}A(t){\boldsymbol{ \zeta }} ,\end{eqnarray} \tag{ D.6 }$$

where

$$\begin{eqnarray}&&A(t)={{\Lambda }^{{\rm T}}}\frac{{\rm d}\Lambda }{{\rm d}t}.\end{eqnarray} \tag{ D.7 }$$

As discussed below, with this choice the phase $\Phi ({\boldsymbol{ \zeta }} ,t)$ will account for most of the kinetic energy of the expanding BEC in the Thomas–Fermi regime.

Substituting equation (D.4) into the original Gross–Pitaevskii equation (D.1) and demanding that $\Lambda (t)$ satisfies the ordinary differential equation

$$\begin{eqnarray}&&{{\Lambda }^{{\rm T}}}(t)\left( \frac{{{{\rm d}}^{2}}\Lambda }{{\rm d}{{t}^{2}}}+{{\Omega }^{2}}(t)\Lambda (t) \right)=\frac{{{\Omega }^{2}}({{t}_{0}})}{{\rm det} \;\Lambda (t)},\end{eqnarray} \tag{ D.8 }$$

with initial conditions

$$\begin{eqnarray}&&\Lambda ({{t}_{0}})={\mathbb{1}},\qquad \dot{\Lambda }({{t}_{0}})=0,\end{eqnarray} \tag{ D.9 }$$

one finds that the equation governing the dynamics of ${{\psi }_{\Lambda }}({\boldsymbol{ \zeta }} ,t)$ takes the following simple and convenient form:

$$\begin{eqnarray}\begin{array}{rcl} {\rm i}\hbar \frac{\partial {{\psi }_{\Lambda }}}{\partial t}{{\mid }_{\zeta }} & = & -\frac{{{\hbar }^{2}}}{2m}\nabla _{\zeta }^{{\rm T}}{{\left( {{\Lambda }^{{\rm T}}}\Lambda \right)}^{-1}}{{\nabla }_{\zeta }}{{\psi }_{\Lambda }} \\ {} & {} & +\frac{1}{{\rm det} \;\Lambda }\left[ \frac{m}{2}{{{\boldsymbol{ \zeta }} }^{{\rm T}}}{{\Omega }^{2}}({{t}_{0}}){\boldsymbol{ \zeta }} +g\;\mid {{\psi }_{\Lambda }}{{\mid }^{2}}-\mu \right]{{\psi }_{\Lambda }}. \\ \end{array}\end{eqnarray} \tag{ D.10 }$$

If one neglects the first term on the right-hand side of equation (D.10), 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:

$$\begin{eqnarray}\begin{array}{rcl} {{\psi }_{\Lambda }}({\boldsymbol{ \zeta }} ,t)\approx {{\psi }_{{\rm TF}}}({\boldsymbol{ \zeta }} ) & = & \frac{1}{\sqrt{g}}\left( \mu -{{{\boldsymbol{ \zeta }} }^{{\rm T}}}{{\Omega }^{2}}({{t}_{0}}){\boldsymbol{ \zeta }} \right)_{+}^{\frac{1}{2}} \\ {} & = & {{\left( \frac{15}{8\pi }\frac{N}{{\rm det} \;\Sigma } \right)}^{\frac{1}{2}}}\left( 1-{{{\boldsymbol{ \zeta }} }^{{\rm T}}}{{\Sigma }^{-1}}{\boldsymbol{ \zeta }} \right)_{+}^{\frac{1}{2}}, \\ \end{array}\end{eqnarray} \tag{ D.11 }$$

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 equation (D.10) 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 equation (D.1) (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 (D.6) 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 equation (D.6). Moreover, the conditions (D.8) and (D.9) fulfilled by $\Lambda (t)$ guarantee that A(t) is symmetric itself. This can be shown as follows. First, one notes that the initial conditions (D.9) imply $A({{t}_{0}})={{A}^{{\rm T}}}({{t}_{0}})=0$. Next, one can check that equation (D.8) implies $\dot{A}(t)-{{\dot{A}}^{{\rm T}}}(t)=0$. It then follows that $A(t)={{A}^{{\rm T}}}(t)\ \forall t$.

Next, we briefly discuss the 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 equation (D.4) into equation (22) and employing the time-independent approximation (D.11) for ${{\psi }_{\Lambda }}({\boldsymbol{ \zeta }} ,t)$. (A more detailed calculation can be found in [37].) In the original unrescaled coordinates it is given by

$$\begin{eqnarray}&&\begin{array}{l} W({\bf x},{\bf p},t) \\ \quad =\frac{1}{{\rm det} \;\Lambda (t)}\int \frac{{{{\rm d}}^{3}}\Delta }{{{(2\pi \hbar )}^{3}}}{{{\rm e}}^{{\rm i}{{{\bf p}}^{{\rm T}}}{\boldsymbol{ \Delta }} /\hbar }}{{{\rm e}}^{-{\rm i}\frac{m}{\hbar }{{{\bf x}}^{{\rm T}}}\dot{\Lambda }{{\Lambda }^{-1}}{\boldsymbol{ \Delta }} }}{{\psi }_{{\rm TF}}}\left( {{\Lambda }^{-1}}({\bf x}+{\boldsymbol{ \Delta }} /2) \right){{\psi }_{{\rm TF}}}\left( {{\Lambda }^{-1}}({\bf x}-{\boldsymbol{ \Delta }} /2) \right). \\ \end{array}\end{eqnarray} \tag{ D.12 }$$

Restoring the rescaled variable ${\boldsymbol{ \zeta }}$ and introducing the analogously rescaled variable ${\boldsymbol{\bar{ \Delta }}} ={{\Lambda }^{-1}}{\boldsymbol{ \Delta }}$ together with the momentum ${{{\bf p}}_{\zeta }}={{\Lambda }^{{\rm T}}}{\bf p}$, it becomes

$$\begin{eqnarray}\begin{array}{rcl} \mathcal{W}({\boldsymbol{ \zeta }} ,{{{\bf p}}_{\zeta }},t) & = & \int \frac{{{{\rm d}}^{3}}\bar{\Delta }}{{{(2\pi \hbar )}^{3}}}{{{\rm e}}^{{\rm i}{\bf p}_{\zeta }^{{\rm T}}{\boldsymbol{\bar{ \Delta }}} /\hbar }}{{{\rm e}}^{-{\rm i}\frac{m}{\hbar }{{{\boldsymbol{ \zeta }} }^{{\rm T}}}A(t){\boldsymbol{\bar{ \Delta }}} }}{{\psi }_{{\rm TF}}}({\boldsymbol{ \zeta }} +{\boldsymbol{\bar{ \Delta }}} /2){{\psi }_{{\rm TF}}}({\boldsymbol{ \zeta }} -{\boldsymbol{\bar{ \Delta }}} /2) \\ {} & = & {{W}_{{\rm TF}}}\left( {\boldsymbol{ \zeta }} ,{{{\bf p}}_{\zeta }}-m\;A(t){\boldsymbol{ \zeta }} \right). \\ \end{array}\end{eqnarray} \tag{ D.13 }$$

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 $({\bf x},{\bf 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}_{{\rm TF}}}({\boldsymbol{ \zeta }} ,{{{\bf p}}_{\zeta }})$ (which also exhibits tails with smaller amplitude oscillations and negative values) corresponds to a region of spatial width $2{{R}_{{\rm TF}}}$ and momentum width $\Delta p\sim \hbar /{{R}_{{\rm 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 section 3. 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. In fact, one can provide a unified treatment for all these cases, as explicitly shown in [37].

An interesting result within the scaling approximation, pointed out in [49, 63], 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 (D.4):

$$\begin{eqnarray}\begin{array}{rcl} \tilde{\psi }({\bf p},t) & = & {{\left( {\rm det} \;\Lambda (t) \right)}^{-\frac{1}{2}}}{{{\rm e}}^{{\rm i}\beta (t)}}\int \frac{{{{\rm d}}^{3}}x}{{{(2\pi \hbar )}^{3/2}}}\;{{{\rm e}}^{-\frac{{\rm i}}{\hbar }{{{\bf p}}^{{\rm T}}}{\bf x}}}{{{\rm e}}^{{\rm i}\frac{m}{2\hbar }\,{{{\boldsymbol{ \xi }} }^{{\rm T}}}A(t){\boldsymbol{ \xi }} }}{{\psi }_{{\rm TF}}}({\boldsymbol{ \xi }} ) \\ {} & = & {{\left( {\rm det} \;\Lambda (t) \right)}^{\frac{1}{2}}}{{{\rm e}}^{{\rm i}\beta (t)}}\int \frac{{{{\rm d}}^{3}}\xi }{{{(2\pi \hbar )}^{3/2}}}\;{{{\rm e}}^{-\frac{{\rm i}}{\hbar }{{{\bf p}}^{{\rm T}}}\Lambda (t){\boldsymbol{ \xi }} }}{{{\rm e}}^{{\rm i}\frac{m}{2\hbar }\,{{{\boldsymbol{ \xi }} }^{{\rm T}}}A(t){\boldsymbol{ \xi }} }}{{\psi }_{{\rm TF}}}({\boldsymbol{ \xi }} ), \\ \end{array}\end{eqnarray} \tag{ D.14 }$$

where in the second equality we simply performed a change of integration variables. One can then use the same procedure employed in the evaluation of equation (C.1) by rewriting equation (D.14) as

$$\begin{eqnarray}\begin{array}{rcl} \tilde{\psi }({\bf p},t) & = & {{{\rm e}}^{\frac{-\,{\rm i}}{2\,m\hbar }{{{\bf p}}^{{\rm T}}}\Lambda \,{{A}^{-1}}{{\Lambda }^{{\rm T}}}{\bf p}}}\int \frac{{{{\rm d}}^{3}}\xi }{{{(2\pi \hbar )}^{3/2}}}\ {{({\rm det} \Lambda )}^{\frac{1}{2}}}\;{{{\rm e}}^{{\rm i}\beta }}\;{{\psi }_{{\rm TF}}}({\boldsymbol{ \xi }} ) \\ {} & {} & \times {{{\rm e}}^{{\rm i}\frac{m}{2\hbar }{{({\boldsymbol{ \xi }} -{{A}^{-1}}{{\Lambda }^{{\rm T}}}{\bf p}/m)}^{{\rm T}}}A({\boldsymbol{ \xi }} -{{A}^{-1}}{{\Lambda }^{{\rm T}}}{\bf p}/m)}}, \\ \end{array}\end{eqnarray} \tag{ D.15 }$$

and using the saddle-point approximation to obtain

$$\begin{eqnarray}&&\tilde{\psi }({\bf p},t)\approx \frac{{{{\rm i}}^{3/2}}{{({\rm det} \;\Lambda )}^{1/2}}\;{{{\rm e}}^{{\rm i}\beta }}}{{{m}^{3/2}}\;{{({\rm det} \;A)}^{1/2}}}\ {{{\rm e}}^{-\frac{{\rm i}}{2m\hbar }{{{\bf p}}^{{\rm T}}}\Lambda \,{{A}^{-1}}{{\Lambda }^{{\rm T}}}{\bf p}}}{{\psi }_{{\rm TF}}}\left( {{A}^{-1}}{{\Lambda }^{{\rm T}}}{\bf p}/m \right),\end{eqnarray} \tag{ D.16 }$$

which is valid when the widths of the Gaussian in equation (D.15) are smaller than the characteristic scales $\Delta {\boldsymbol{ \xi }}$ of ${{\psi }_{{\rm TF}}}({\boldsymbol{ \xi }} )$, i.e. when $\sqrt{\hbar /(m\mid {{A}_{ii}}\mid )}\ll \Delta {{\xi }^{i}}\sim R_{{\rm TF}}^{(i)}$ 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 $\Lambda (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_{{\rm TF}}^{(i)}\gg \sqrt{\hbar /(m\mid {{B}_{ii}}\mid \mid {{\Lambda }_{ii}}\mid )}$. (Here we considered the basis, where the matrices B and $\Lambda (t)$ both diagonalize, which is always possible for nonrotating traps and coincides with the basis where the matrix ${{\Omega }^{2}}(t)$ diagonalizes.) Note, for comparison, that Heisenbergʼs uncertainty principle implies $R_{{\rm TF}}^{(i)}\mid {{\Lambda }_{ii}}\mid \;m\;R_{{\rm TF}}^{(i)}\mid {{B}_{ii}}\mid \gtrsim \hbar$, i.e. $R_{{\rm TF}}^{(i)}\gtrsim \sqrt{\hbar /(m\mid {{B}_{ii}}\mid \mid {{\Lambda }_{ii}}\mid )}$. 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 equation (D.13). 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 equation (D.16) 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.

The result of the previous subsection can be combined with that of appendix C to analyze the late-time behavior of an expanding BEC. In order to do so, we consider an intermediate time t₁ at which the density of the condensate is already low enough that we can neglect the nonlinear interactions. Hence, at this point the matrix $\Lambda (t)$ has already entered the linear regime. One can use equation (D.16) to evaluate $\tilde{\psi }({\bf p},{{t}_{1}})$, and substituting the expression for $\Lambda (t)$ linear in time, one gets:

$$\begin{eqnarray}&&\tilde{\psi }({\bf p},{{t}_{1}})\approx \frac{{{({\rm det} \;B)}^{-\frac{1}{2}}}}{{{(-{\rm i}\;m)}^{3/2}}}\ {{{\rm e}}^{{\rm i}\beta }}\;{{{\rm e}}^{-\frac{{\rm i}}{\hbar }{{{\bf p}}^{2}}\Delta {{t}_{1}}/2m}}{{\psi }_{{\rm TF}}}\left( {{B}^{-1}}{\bf p}/m \right),\end{eqnarray} \tag{ D.17 }$$

with $\Delta {{t}_{1}}=({{t}_{1}}-{{t}_{0}})$. Furthermore, for times $t\geqslant {{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 equation (C.1) applied to this case and taking equation (D.17) for the initial wave function at time t₁ in momentum representation, one gets

$$\begin{eqnarray}\begin{array}{rcl} \psi ({\bf x},t) & = & \int \frac{{{{\rm d}}^{3}}p}{{{(2\pi \hbar )}^{3/2}}}\;{{{\rm e}}^{{\rm i}\beta }}\;{{{\rm e}}^{\frac{{\rm i}}{\hbar }{{{\bf p}}^{{\rm T}}}{\bf x}}}{{{\rm e}}^{-\frac{{\rm i}}{\hbar }({{{\bf p}}^{2}}/2m)(t-{{t}_{1}})}}\tilde{\psi }({\bf p},{{t}_{1}}) \\ {} & = & \frac{{{({\rm det} \;B)}^{-\frac{1}{2}}}}{{{(-{\rm i}\;m)}^{3/2}}}\int \frac{{{{\rm d}}^{3}}p}{{{(2\pi \hbar )}^{3/2}}}\;{{{\rm e}}^{{\rm i}\beta }}\;{{{\rm e}}^{\frac{{\rm i}}{\hbar }{{{\bf p}}^{{\rm T}}}{\bf x}}}{{{\rm e}}^{-\frac{{\rm i}}{\hbar }{{{\bf p}}^{2}}\,\Delta t/2m}} \\ {} & {} & \times {{\psi }_{{\rm TF}}}({{B}^{-1}}{\bf p}/m), \\ \end{array}\end{eqnarray} \tag{ D.18 }$$

where in the second equality we combined the two phase factors with exponent proportional to ${{{\bf p}}^{2}}$ and introduced $\Delta t=(t-{{t}_{1}})+\Delta {{t}_{1}}=(t-{{t}_{0}})$. Next, rewriting equation (D.18) as done in equation (C.2) and using the same kind of saddle-point approximation that led to equation (C.3), one finally obtains:

$$\begin{eqnarray}&&\psi ({\bf x},t)\approx {{\left( {\rm det} (B\Delta t) \right)}^{-\frac{1}{2}}}{{{\rm e}}^{{\rm i}\beta }}\;{{{\rm e}}^{\frac{{\rm i}}{\hbar }\frac{m}{2\Delta t}{{{\bf x}}^{2}}}}\;{{\psi }_{{\rm TF}}}\left( {{B}^{-1}}{\bf x}/\Delta t \right),\end{eqnarray} \tag{ D.19 }$$

which coincides with the scaling solution (D.4) when C can be neglected compared to $B\;\Delta t$. The saddle-point approximation is valid provided that the typical scales of variation $\Delta {\bf p}$ of the function $\tilde{\psi }({\bf p},{{t}_{1}})$ satisfy the condition $\sqrt{m\hbar /\Delta t}\ll \Delta {{p}^{i}}$, which implies $\sqrt{\hbar /(m\mid {{B}_{ii}}{{\mid }^{2}}\Delta t)}\ll \mid {{B}_{ii}}{{\mid }^{-1}}\Delta {{p}^{i}}/m\sim R_{{\rm TF}}^{(i)}$ for $i=1,2,3$. For sufficiently late times (i.e. sufficiently large $\Delta t$) this condition will be fulfilled as long as the condition justifying the derivation of equation (D.16) 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₁ at which the expansion dynamics has already entered the linear regime, it will be valid at arbitrarily late times.