Tools for designing atom interferometers in a microgravity environment

TL;DR

This paper introduces a variational model using Gaussian clouds for rapid design and analysis of atom interferometers in microgravity, aiding in optimizing experiments to measure gravitational constants.

Contribution

The paper develops a flexible variational model based on Gaussian clouds and applies it to design and optimize atom interferometers in microgravity environments.

Findings

Model constrains experimental parameter space.

Demonstrates how to extract G from interference patterns.

Provides methods to improve measurement sensitivity.

Abstract

We present a variational model suitable for rapid preliminary design of atom interferometers in a microgravity environment. The model approximates the solution of the 3D rotating--frame Gross--Pitaevskii equation (GPE) as the sum of Nc Gaussian clouds. Each Gaussian cloud is assumed to have time--dependent center positions, widths, and linear and quadratic phase parameters. We applied the Lagrangian Variational Method (LVM) with this trial wave function to derive equations of motion for these parameters that can be adapted to any external potential. We also present a 1D version of this variational model. As an example we apply the model to a 1D atom interferometry scheme for measuring Newton's gravitational constant, G, in a microgravity environment. We show how the LVM model can (1) constrain the experimental parameter space size, (2) show how the value of G can be obtained from the experimental conditions and interference pattern characteristics, and (3) show how to improve the sensitivity of the measurement and construct a preliminary error budget.

Equations

$$i\hslash \frac{\partial \Phi }{\partial t}$$

$$\mathcal{L}[{\Phi }^{\ast }]$$

$$\sum _{\eta =x,y,z,t}\frac{\partial }{\partial \eta }\left(\frac{\partial \mathcal{L}}{\partial {\Phi }_{\eta }^{\ast }}\right)-\frac{\partial \mathcal{L}}{\partial {\Phi }^{\ast }}=0,\mathrm{where}{\Phi }_{\eta }\equiv \frac{\partial \Phi }{\partial \eta }.$$

$${\Phi }^{\mathrm{trial}}(\mathbf{r},t)={\Phi }^{\mathrm{trial}}(q_1(t),\dots ,q_n(t);\mathbf{r})$$

$$L(q_1(t),\dots ,q_n(t))=\int d^{3}r\mathcal{L}[{\left({\Phi }^{\mathrm{trial}}\right)}^{\ast }]$$

$$\frac{d}{dt}\left(\frac{\partial L}{\partial {\dot{q}}_k}\right)-\frac{\partial L}{\partial q_k}=0,k=1,\dots ,n$$

$$E_0\equiv \frac{{\hslash }^2}{2M{L}_0^2}\mathrm{and}{T}_0\equiv \frac{\hslash }{E_0}=\frac{2M{L}_0^2}{\hslash }.$$

$$\bar{x}\equiv \frac{x}{L_0}\bar{y}\equiv \frac{y}{L_0}\bar{z}\equiv \frac{z}{L_0}\mathrm{and}\bar{t}\equiv \frac{t}{T_0}.$$

$$\Phi (\mathbf{r},t)=L_0^{-3/2}\Psi (\bar{\mathbf{r}},\bar{t}).$$

$$i\frac{\partial \Psi }{\partial \bar{t}}$$

$${\bar{\mathcal{L}}}^{(3D)}\left[{\Psi }^{\ast }\right]$$

$$\sum _{\eta =x,y,z,t}\frac{\partial }{\partial \bar{\eta }}\left(\frac{\partial {\bar{\mathcal{L}}}^{(3D)}}{\partial {\Psi }_{\bar{\eta }}^{\ast }}\right)-\frac{\partial {\bar{\mathcal{L}}}^{(3D)}}{\partial {\Psi }^{\ast }}=0.$$

$${\bar{L}}^{(3D)}\left[\mathbf{x},\mathbf{w},\mathbf{\alpha },\mathbf{\beta }\right]$$

$$i\frac{\partial \psi }{\partial \bar{t}}=-\frac{{\partial }^2\psi }{\partial {\bar{x}}^2}+{\bar{V}}_{\mathrm{ext}}(\bar{x},\bar{t})\psi +{\bar{g}}_{1D}N{\left|\psi \right|}^2\psi .$$

$${\bar{\mathcal{L}}}^{(1D)}\left[{\psi }^{\ast }\right]$$

$$\frac{\partial }{\partial \bar{x}}\left(\frac{\partial {\bar{\mathcal{L}}}^{(1D)}}{\partial {\psi }_{\bar{x}}^{\ast }}\right)+\frac{\partial }{\partial \bar{t}}\left(\frac{\partial {\bar{\mathcal{L}}}^{(1D)}}{\partial {\psi }_{\bar{t}}^{\ast }}\right)-\frac{\partial {\bar{\mathcal{L}}}^{(1D)}}{\partial {\psi }^{\ast }}=0.$$

$${\bar{L}}^{(1D)}\left[\mathbf{x},\mathbf{w},\mathbf{\alpha },\mathbf{\beta }\right]={\int }_{-\infty }^{\infty }d\bar{x}{\bar{\mathcal{L}}}^{(1D)}\left[{\psi }^{\ast }\right].$$

$$\Psi (\bar{\mathbf{r}},\bar{t})=\frac{1}{\sqrt{N_c}}\sum _{j=1}^{N_c}{A}_j(\bar{t})e^{f_j(\bar{\mathbf{r}},\bar{t})+i{\bar{\mathbf{k}}}_j\cdot \bar{\mathbf{r}}}$$

$$f_j(\bar{\mathbf{r}},\bar{t})$$

$$\psi (\bar{x},\bar{t})=\frac{1}{\sqrt{N_c}}\sum _{j=1}^{N_c}{A}_j(\bar{t})e^{f_j(\bar{x},\bar{t})+i{\bar{k}}_j\bar{x}}$$

$$f_j(\bar{x},\bar{t})=-\frac{(\bar{x}-{\bar{x}}_j(\bar{t}){)}^2}{2{\bar{w}}_j^2(\bar{t})}+i{\bar{\alpha }}_j(\bar{t})\bar{x}+i{\bar{\beta }}_j(\bar{t}){\bar{x}}^2$$

$$1$$

$$A_j^2(\bar{t}){\pi }^{3/2}{\bar{w}}_{jx}(\bar{t}){\bar{w}}_{jy}(\bar{t}){\bar{w}}_{jz}(\bar{t})=1,j=1,\dots ,N_c.$$

$$A_j^2(\bar{t}){\pi }^{1/2}{\bar{w}}_j(\bar{t})=1,j=1,\dots ,N_c.$$

$${\bar{L}}^{(3D)}$$

$${\bar{U}}_{\mathrm{ext}}^{(3D)}(\mathbf{x},\mathbf{w})$$

$${\bar{U}}_{\mathrm{int}}^{(3D)}(\mathbf{x},\mathbf{w})$$

$${\bar{U}}_{\mathrm{ext}}^{(1D)}(\mathbf{x},\mathbf{w})$$

$${\bar{U}}_{\mathrm{int}}^{(1D)}(\mathbf{x},\mathbf{w})$$

$${\bar{L}}^{(1D)}$$

$${\ddot{\bar{x}}}_j$$

Full text

I Introduction

Atom interferometers (AIs) have been used in many applications Tino and Kasevich (2014); Cronin et al. (2009). In particular this is the case for gravitational measurements, such as the determination of Earth’s gravitational acceleration, $g$, Peters et al. (1999) and its gradients McGuirk et al. (2002); Sorrentino et al. (2010). Atom interferometers have also been proposed for use in Sagnac gyroscopes Gustavson et al. (2000); Canuel et al. (2006), in testing general relativity Fray et al. (2004); Dimopoulos et al. (2007); Tarallo et al. (2014), in searching for dark energy Burrage et al. (2015), and measuring gravity waves Dimopoulos et al. (2008); Graham et al. (2013). The Newtonian constant of gravitation, $G$, has also been measured with AIs Rosi et al. (2014); Lamporesi et al. (2008).

In a typical atom interferometer involving a gaseous Bose–Einstein condensate (BEC), the basic sequence consists of splitting the condensate into two or more clouds that then experience different environments while separated. Some differing property of these environments generates a phase difference between the clouds while they are separated. The phase difference produces an interference pattern when the clouds are reunited and this pattern can be used to infer the value of this property.

An emerging area of AI applications is to systems that are freely falling in a gravitational field. Such systems offer the prospect of much longer interrogation times than are possible in fixed terrestrial laboratories, where atoms are subject to constant vertical acceleration. Such “microgravity” environments have been created in drop towers van Zoest et al. (2010), sounding rockets Seidel et al. (2016), and recently aboard the International Space Station (ISS).

The National Aeronautics and Space Administration (NASA) has deployed its Cold Atom Laboratory (CAL) Aveline et al. (2017) to the ISS enabling future atom–interferometry experiments to be performed on Bose–Einstein condensates there. This will provide a user facility in which atom–interferometry experiments can be conducted in a microgravity environment.

The basic AI sequence of splitting, reuniting, and splitting the BEC again to produce an interference pattern can be viewed as a sequence of sudden splits of the condensate alternating with periods where the condensate clouds evolve freely according to their local environments. In general, the performance of AIs of this type depends on several factors including the number of interferometer pathways, the area of the interferometer, and the interrogation time (the time the clouds spend in the different environments).

AI performance can be enhanced both by increasing any or all of these quantities. Specific performance enhancements can also be achieved by designing different AI sequences as was done recently with neutron interferometers Pushin et al. (2009). Differential AI experiments, where the same AI sequence is carried out on identical condensates at the same time, can also be used in making precision measurements.

Enhancing AI performance by increasing interrogation time is especially possible in a microgravity environment. The extended interrogation times available in microgravity are also valuable for precision measurement applications. Microgravity environments in general and the CAL in particular are important platforms for making precision measurements. To further this end, designing new AI sequences that take advantage of microgravity would be highly desirable.

Designing new AI sequences involving BECs suitable for precision measurements conducted in microgravity environments presents a demanding numerical challenge. Detailed modeling of condensate behavior for AI sequences having long interrogation times and/or multiple pathways requires large amounts of computer time and storage assuming that the condensate obeys the Gross-Pitaevskii equation (GPE). Many different possible sequences would also need to be simulated and the results evaluated. Thus a method for rapid simulation of a given AI sequence would be an extremely useful tool for AI design. Such a tool would also be helpful in determining how the quantity measured by the AI could be extracted from the resulting interference patterns.

This paper presents a tool suitable for preliminary design of candidate AI sequences. The tool can provide rapid approximate modeling of condensate behavior between the splits of an AI sequence where this evolution is assumed to be governed by the GPE. This model is based on the Lagrangian Variational Method (LVM) where the assumed trial wave function approximates the GPE solution as a set of $N_{c}$ separate but possibly overlapping Gaussian clouds. We emphasize that the final choice of an AI sequence would have to be validated by simulating it with the GPE.

In section II we describe the general Lagrangian Variational Method, introduce a set of scaled units suitable for computation, and present the trial wave functions for both the one-dimensional and three-dimensional versions of the method. In section III we derive the equations of motion for the variational parameters and show that these can be cast in terms of the space and width gradients of a “variational potential” which is just the expectation value of the external potential over the trial wave function.

In section IV we describe an example of an AI scheme that could be used to measure $G$ in microgravity. In section V we apply the LVM model to this AI scheme. We compare a GPE simulation of this scheme with the LVM simulation. We also use the model to obtain an expression for the LVM interference pattern, show how $G$ can be extracted from this pattern, and how the model can facilitate preliminary design of an AI sequence. Finally, we present a summary of the work in section VI.

II The Lagrangian Variational Method and Trial Wave Function

II.1 The general Lagrangian variational method

The Lagrangian Variational Method provides approximations to the solutions of the time–dependent Gross–Pitaevskii equation Perez-Garcia et al. (1997); Edwards et al. (2012); Infeld et al. (2006); Salasnich et al. (2002); Salasnich and Malomed (2009). In three dimensions and in the rotating frame, this equation has the form Pitaevskii and Stringari (2003):

[TABLE]

where $\Phi({\bf r},t)$ is the condensate wave function, $M$ is the mass of a condensate atom, $N$ is the number of atoms in the condensate, $g_{3D}=4\pi\hbar^{2}a_{s}/M$ measures the strength of the atom–atom scattering with $a_{s}$ being the scattering length, $V_{\rm ext}({\bf r},t)$ is the potential exerted on a condensate atom by external fields, and ${\bf\Omega}$ is the angular velocity of the rotating frame.

The LVM is based on the fact that the GPE can be derived as the Euler–Lagrange equation of motion produced by the following Lagrangian density:

[TABLE]

where $\hat{\bf L}$ is the single–particle angular momentum operator.

This Lagrangian density along with the following Euler–Lagrange equation of motion produces the GPE:

[TABLE]

The Lagrangian Variational Method consists of devising a trial wave function,

[TABLE]

where the $\{q_{i}(t)\}$, $i=1,\dots,n$ are variational parameters that only depend on the time, $t$. The equations of motion of these variational parameters are derived by computing the ordinary Lagrangian:

[TABLE]

and then using the standard Euler–Lagrange equations,

[TABLE]

to produce an equation of motion associated with each variational parameter.

II.2 Scaled units

We can simplify the equations produced by the above method by introducing a set of units appropriate to the problem and a set of scaled variables (both independent and dependent). The scaled variables are defined by first choosing a length unit, $L_{0}$, and then defining energy and time units in terms of $L_{0}$:

[TABLE]

We then introduce scaled variables which are generally denoted by barred quantities. These consist of scaled space and time coordinates:

[TABLE]

We also introduce the scaled condensate wave function for the solution of the GPE:

[TABLE]

We can express the original GPE in terms of scaled quantities and this can be done for the Lagrangian density and its associated Euler–Lagrange equation as well.

In terms of scaled quantities, the rotating–frame GPE becomes:

[TABLE]

where we are considering the form of the rotating–frame GPE for the special case of rotation around the $z$ axis. In the above we have $g_{3D}\equiv\bar{g}_{3D}E_{0}L_{0}^{3}$, $V_{\rm ext}({\bf\bar{r}},\bar{t})=\bar{V}_{\rm ext}({\bf r},t)E_{0}$ and $\Omega_{z}=\bar{\Omega}_{z}/T_{0}$. The scaled Lagrangian density for this version of the GPE becomes

[TABLE]

and the scaled Euler–Lagrange equation is given by

[TABLE]

It is straightforward to insert Eq. (11) into Eq. (12) and obtain Eq. (10). The 3D Lagrangian is computed by integrating over all 3D space:

[TABLE]

In the above have denoted the dependence of $\bar{L}$ on the different types of variational parameters that will appear in the trial wave function defined in the next section. These are: the cloud center coordinates by $\bm{x}\equiv\left(\bar{x}_{1},\dots,\bar{z}_{N_{c}}\right)$, the cloud widths denoted by $\bm{w}\equiv\left(\bar{w}_{1x},\dots,\bar{w}_{N_{c}z}\right)$, the linear phase coefficients by $\bm{\alpha}\equiv\left(\alpha_{1x},\dots,\alpha_{N_{c}z}\right)$ and the quadratic phase coefficients by $\bm{\beta}\equiv\left(\beta_{1x},\dots,\beta_{N_{c}z}\right)$.

The LVM formulation is straightforward to apply in one–dimension where we denote the 1D GPE solution by $\phi(x,t)$. For scaled units we will write $\phi(x,t)\equiv L_{0}^{-1/2}\psi(\bar{x},\bar{t})$. In terms of scaled quantities, the 1D GPE becomes:

[TABLE]

where $g_{1D}\equiv\bar{g}_{1D}E_{0}L_{0}$ and $V_{\rm ext}(\bar{x},\bar{t})=\bar{V}_{\rm ext}(x,t)/E_{0}$. The 1D scaled Lagrangian density becomes

[TABLE]

and the 1D scaled Euler–Lagrange equation is given by

[TABLE]

We can calculate the Lagrangian associated with this trial wave function by integrating $\bar{\cal L}^{(1D)}$ over all 1D space:

[TABLE]

This Lagrangian will depend only on the variational parameters contained in the trial wave function. The 1D equations of motion are found with the standard Euler–Lagrange equations. Next we will present the 1D and 3D LVM trial wave functions.

II.3 $N_{c}$ 3D LVM trial wave function

In the 3D, $N_{c}$–gaussian model we take the trial wave function to be a sum of $N_{c}$ three–dimensional Gaussian clouds. The $j^{th}$ cloud has its own initial momentum, ${\bf\bar{k}}_{j}$ and set of variational parameters. These parameters are the cartesian coordinates of the cloud center: $\bar{x}_{j}$, $\bar{y}_{j}$, and $\bar{z}_{j}$; the widths along the $x$, $y$, and $z$ directions: $\bar{w}_{jx}$, $\bar{w}_{jy}$, and $\bar{w}_{jz}$; the linear phase coefficients along the $x$, $y$, and $z$ directions: $\bar{\alpha}_{jx}$, $\bar{\alpha}_{jy}$, and $\bar{\alpha}_{jz}$; and the quadratic phase coefficients along the $x$, $y$, and $z$ directions: $\bar{\beta}_{jx}$, $\bar{\beta}_{jy}$, and $\bar{\beta}_{jz}$. The $j^{th}$ cloud also has its own normalization coefficient, $A_{j}$, which will be eliminated by fixing the number of atoms in each cloud.

The mathematical form (in scaled units) of the trial wave function is the following:

[TABLE]

where

[TABLE]

This trial wave function will be used below in deriving the 3D Lagrangian function from which the 3D equations of motion will be obtained.

II.4 $N_{c}$ 1D LVM trial wave function

In the 1D, $N_{c}$–gaussian–cloud model we take the trial wave function to be a sum of $N_{c}$ one–dimensional Gaussian clouds. The $j^{th}$ cloud has its own initial momentum, $\bar{k}_{j}$ and set of variational parameters. These parameters are the cartesian coordinate of the cloud center, $\bar{x}_{j}$, the cloud width, $\bar{w}_{j}$, the linear phase coefficient, $\bar{\alpha}_{j}$, and the quadratic phase coefficient, $\bar{\beta}_{j}$. The $j^{th}$ cloud also has its own normalization coefficient, $A_{j}$, which will be eliminated by fixing the number of atoms in each cloud.

The mathematical form (in scaled units) of the trial wave function is the following:

[TABLE]

where

[TABLE]

III The LVM Equations of Motion

In this section we derive the 3D and 1D LVM equations of motion. We first describe our assumed constraints placed on the trial wave function, then we sketch the derivation of the 3D and 1D Lagrangians for the given trial wave functions and finally sketch the derivation of the final LVM equations of motion. Most of the details are relegated to appendices. The terms in these Lagrangians containing the external potential $\bar{V}_{\rm ext}$ only depend on the cloud centers, ${\bf x}$, and widths, ${\bf w}$, no matter the form of the potential. This is due solely to our choice of trial wave function. Thus the 1D and 3D Lagrangians can be written in terms of a “variational” external potential, $U_{\rm ext}({\bf x},{\bf w})$. This enables the terms in the final equations of motion that account for the external potential to be written in terms of the ${\bf x}$ and ${\bf w}$ gradients of this potential. Hence the final equations of motion can be written in a form independent of the particular $V_{\rm ext}$.

III.1 Constraints on the wave function

Here we make several assumptions about the physical system which have material effects on the values of the variational parameters. These are as follows:

We assume that each of the $N_{c}$ clouds are moving at sufficiently different velocities such that any integral of a quantity containing a factor like $e^{i\left(\bar{k}_{j\eta}-\bar{k}_{j^{\prime}\eta}\right)\bar{\eta}}$ where $j\neq j^{\prime}$ can be neglected. If the clouds move with sufficiently different velocities, these factors will be rapidly oscillating and their integrals can be neglected. 2. 2.

The amplitudes, $A_{j}(\bar{t})$, are all real. Here we assume that the phase of each cloud’s wave function is quadratic so that there is no phase is subsumed by the amplitudes. The quadratic phase enables the cloud centers to move and the cloud widths to expand and contract. 3. 3.

The number of atoms in each cloud is fixed. Clouds do not exchange atoms. This plus the normalization condition, fixes a relationship (derived below) between $A_{j}$ and the widths $\bar{w}_{j\eta}$ where $\eta=x,y,z$.

We can use these assumptions plus the normalization condition on the trial wave function to derive conditions that constrain the values of the $A_{j}$.

To find these conditions we require that the full trial wave function be normalized to unity:

[TABLE]

where we have used assumptions 1 and 2 above to reduce the normalization condition to a sum of Gaussian integrals which are easily evaluated. This last expression is the condition for the trial wave function to be normalized. However, our assumption that the number of atoms in each cloud is fixed adds a further restriction to the above expression. That is that each cloud is individually normalized. This gives, finally,

[TABLE]

These constraints together automatically satisfy Eq. (22) and enable the elimination of all of the $A_{j}$ in the final 3D Lagrangian. A similar condition on the norm of each cloud in the 1D case can be derived and the result is

[TABLE]

Using these conditions we can now derive the 3D and 1D Lagrangians.

III.2 The LVM model Lagrangians

The LVM Lagrangians for the particular choice of trial wave function given above can be derived by inserting Eq. (18) into Eq. (13) and performing the required integrations. The details are given in Appendix A. The result for the 3D Lagrangian is

[TABLE]

The variational potentials are defined as

[TABLE]

The derivation of the Lagrangian for the 1D case is very similar. The major difference is the absence of a rotating–frame term. Thus we will simply define the 1D variational potentials and present the 1D Lagrangian result.

The 1D external variational potential can be written as

[TABLE]

Where we have used constraints 1 and 3 to simplify the integrals. In order to apply this 1D model to a particular system, this potential must be calculated. Derivatives of $\bar{U}_{\rm ext}$ appear in the equations of motion.

The 1D interaction variational potential is defined by

[TABLE]

We can use these definitions to write the final form of the 1D Lagrangian as follows.

[TABLE]

Explicit expressions for the 1D and 3D interaction potentials (which are always the same) are derived in Appendix B. With the 3D and 1D Lagrangians in hand we can now derive the equations of motion for the variational parameters.

III.3 LVM model equations of motion

The 3D equations of motion for the Gaussian center coordinates, widths, and linear and quadratic phase parameters can obtained in two steps: (1) derive the E–L equation of motion for each variational parameter by using the 3D Lagrangian (Eq. (25)) in the standard E–L equation of motion (Eq. (6)) to obtain a first–order differential equation in time, (2) take the second time derivative of the equations containing $\dot{\bar{\eta}}_{j}$ and $\dot{\bar{w}}_{j\eta}$ where ($j=1,\dots,N_{c}$ and $\eta=x,y,z$) and use the other equations to eliminate the other variables from these equations. The equations for the centers and widths together form a closed system.

The LVM 3D equations of motion thus consist of a pair of second–order ordinary differential equation for the cloud centers and widths as well as expressions for the $\bar{\beta}_{j\eta}$ and the $\bar{\alpha}_{j\eta}$ in terms of the centers, widths and their first derivatives. The derivation of these equations can be found in Appendix C.

[TABLE]

In the above, $\bar{U}^{(3D)}\equiv\bar{U}^{(3D)}_{\rm ext}+\bar{U}^{(3D)}_{\rm int}$. The equations for the cloud centers and cloud widths (Eqs. (31a), (31b), (31c), and (31d)) form a closed set that contain only the $\bar{\eta}_{j}$, $\dot{\bar{\eta}}_{j}$, $\bar{w}_{j\eta}$, and $\dot{\bar{w}}_{j\eta}$. Once these quantities are obtained, all of the other variational parameters can be calculated.

The equations of motion for the 1D case are derived similarly. We will simply state the result here. They consist of a pair of second–order ordinary differential equations for the cloud centers and widths as well as expressions for the $\bar{\beta}_{j}$ and the $\bar{\alpha}_{j}$ in terms of the centers, widths and their time derivatives:

[TABLE]

Similar to the 3D case, $\bar{U}^{(1D)}\equiv\bar{U}^{(1D)}_{\rm ext}+\bar{U}^{(1D)}_{\rm int}$. The equations for the cloud centers and cloud widths (Eqs. (32a) and (32b)) form a closed set that contain only the $\bar{x}_{j}$, $\dot{\bar{x}}_{j}$,$\bar{w}_{j}$, and $\dot{\bar{w}}_{j}$. Once these quantities are obtained, all of the other variational parameters can be calculated.

These equations of motion for the variational parameters of the 1D and 3D trial wave functions form the central result of this paper. We note one more time that these equations hold for any external potential, $\bar{V}_{\rm ext}$. The 1D and 3D $N_{c}$–Gaussian–cloud wave function ansatz has several areas of application including atom interferometery involving BECs and BECs moving in waveguides.

IV AI Measurement of $G$ in Microgravity

As an example of the use of the LVM model, we describe how it could be applied to the design of a measurement of Newton’s universal gravitational constant, $G$, in a microgravity environment. Thus we consider the idealized 1D atom–interferometry sequence shown in Fig. 1 where two pieces of a BEC are separated, differentially pulled upon by a source mass (SM), recombined, and split again. In this section we describe this AI sequence in detail.

To extract the value of $G$ requires an interference pattern where only the effects of the SM are present. This can be obtained by repeating the above AI sequence many times with the SM present so that averaging over the interference patterns produced would wash out random effects of the environment. Systematic effects of the environment as well as the effect of the SM would still be present in the averaged pattern. To remove systematic effects of the environment the SM would then be taken away and the AI sequence repeated many more times and these SM–absent interference patterns would again be averaged. The difference between the averaged SM–present and SM-absent patterns would leave a pattern in which only the effects of the SM are present. The pattern due to the SM where environmental effects are neglected can be approximately simulated by the variational method described here.

The basic idea of the AI sequence is that the relative velocity developed between the two condensate pieces due to the differential acceleration produced by the gravity of the source mass. This relative velocity is imprinted on the final–state interference pattern since the velocity distribution is proportional to the condensate phase gradient.

The value of $G$ can then be extracted from the interference data. This sequence is imagined as the essential step in an AI measurement of $G$. It is assumed to be carried out in a microgravity environment such as NASA’s Cold Atom Laboratory aboard the International Space Station. A more sophisticated version of this sequence would be required for a precision measurement of $G$ under such conditions.

The basic sequence begins with a condensate formed at the center of a 1D harmonic trap with frequency $\omega_{T}$ with a SM present (Fig. 1(a)). The overall sequence consists of an “initial split” where this initial condensate is separated into two pieces which are then rejoined and a “final split” where the rejoined pieces are split again and the interference pattern is imaged. The initial split (IS) has three phases as we now describe.

The speed $v$ imparted to the condensate pieces is assumed in this work to come from applying identical counterpropagating laser beams to the condensate. An atom will receive a momentum kick by absorbing a photon from one beam and emitting it into the other beam. The momentum imparted to a condensate atom is $2\hbar k_{L}$ where $k_{L}=2\pi/\lambda_{L}$ where $\lambda_{L}$ is the wavelength of the laser. Laser pulses can be engineered so that this basic momentum–kick process happens $m$ times so that the velocity kick the atoms receive is $v=4\pi\hbar m/(M\lambda_{L})$ where $M$ is the mass of the condensate atom. In the examples given below we assume that $m=8$ so that, for 87Rb atoms, $v=9.14\times 10^{2}$ m/s.

In IS phase 1 the condensate is separated into two equal pieces by a sequence of optical–lattice pulses Edwards et al. (2010) the total effect of which is to change the velocity of one condensate piece by $+v$ and the other by $-v$ (Fig. 1(b)). The initially motionless condensate pieces thus fly apart with velocities $-v$ (blue circle to the left) and $+v$ (red circle to the right) and both come to a stop at the turning points of the harmonic potential after one quarter of the trap period. We call this time $t=T_{1/4}$ and is shown in Fig. 1(c).

In IS phase 2 the harmonic trap is turned off at time $t=T_{1/4}$ ideally leaving the two pieces motionless. The system is then allowed to evolve during the time interval $T_{1/4}\leq t\leq T_{1/4}+T_{\rm w}$. We call this evolution time, $T_{\rm w}$, the “wait period”. This is shown in Fig. 1(d). During this time the gravitational pull of the SM causes a differential acceleration of the two condensate pieces since one piece is further away from the SM than the other. This differential acceleration causes the two pieces to develop a relative velocity.

In IS phase 3 the trap is then turned back on which brings the condensate pieces back together during the time interval $T_{1/4}+T_{\rm w}\leq t\leq 2T_{1/4}+T_{\rm w}$ as shown in Fig. 1(e). Once they overlap again the blue (left) piece will be moving at velocity $+v$ plus a small difference due to the SM and the red (right) piece will be moving at $-v$ plus a small difference. The initial split takes place during the time interval $0\leq t\leq 2T_{1/4}+T_{\rm w}\equiv T_{\rm is}$ and shown in Figs. 1(a)–(e).

Once the two pieces are overlapped again the same set of optical–lattice pulses is applied as before (Fig. 1(f)). This causes the red condensate piece to split in two with a piece moving at velocity approximately $-2v$ and the other at zero velocity plus a small amount due to the pull of the SM. The blue half is also split into a piece moving at approximate velocity $+2v$ and the other approximately motionless except again for a small deviation. These two nearly motionless pieces will remain overlapped and will have the small relative velocity that was developed during the wait time when the trap was turned off. We assume that the wait time is much longer than a quarter period of the harmonic trap ($T_{\rm w}\gg T_{1/4}$).

Finally, the two fast condensate pieces will fly away leaving the nearly motionless pieces behind as shown in Fig. 1(g). We will refer to this sequence, shown in Figs. 1(f)–(g), as the “final split”. Imaging this middle cloud will leave an interference pattern due to the relative velocity of the two condensate pieces and this pattern will be due only to the differential gravitational pull of the SM during the wait time $T_{W}$.

V Applying the LVM model to the AI sequence

V.1 Comparison of the LVM model and GPE solutions

In this section we show how the LVM model can be used to analyze the atom interferometry sequence described above. We emphasize here that the LVM model is not sufficiently accurate to be used in the analysis of a precision measurement rather it is to be used in the preliminary design of the experiment. To show this, we simulated the evolution of the condensate wave function during the initial split of this AI sequence using the 1D time–dependent GPE and the 1D LVM model described earlier. The comparison is shown in Fig. 2. The details of how the simulation was conducted can be found in Appendix D.

This comparison shows that, while the Gaussian–approximate trial wave functions only qualitatively agree their GPE counterparts, the motion of the GPE and LVM wave packet centers and widths track very well. The variational wave function will only agree with the exact solution to the extent that variational wave function can be modified to fit the exact solution by varying its parameters.

Modeling performed for the final design of a precision AI measurement will require the numerical solution of the full 3D Gross–Pitaevskii equation. Since experimental conditions for precision measurements are typically extreme, numerical solution of the 3D GPE will require a substantial effort. This effort usually makes solving the 3D GPE too expensive for the purposes of preliminary AI design. We shall elaborate on this point below.

The purpose of the LVM model is to facilitate the preliminary design of a precision measurement. It can serve this purpose in several ways. First, it can provide a simple and physically intuitive picture of how the measured quantity can be extracted from the experimental data. In some cases an approximate analytical expression for the measured quantity in terms of experimental parameters can be derived. Second, the model can be used to estimate the sensitivity of the measurement. Finally, this expression can be used to speed up significantly the time spent on the preliminary design of an AI sequence by providing estimates of the values of these parameters thus constraining the size of the experimental parameter space. We illustrate these features of the LVM model by applying it to the AI sequence described above.

V.2 Approximating the interference pattern using the LVM model

The equations of motion for the case when an harmonic trap with frequency $\omega_{T}$ is present along with a source mass of mass $M_{SM}$ and located at $x_{SM}$ can be found by computing the derivatives found in Eqs. (32a) and (32b). The external varational potential for this case, given by Eq. (96), is derived in Appendix E and the interaction variational potential is given in Eq. (68) of Appendix B.

The resulting equations of motion have a striking similarity to equations of motion produced by Newton’s second law. They can be written as follows:

[TABLE]

The terms $X_{j}({\bf x},{\bf w})$ and $W_{j}({\bf x},{\bf w})$ represent “repulsion forces” due to the overlap of different clouds and are given by

[TABLE]

and

[TABLE]

where $j=1,2$.

The LVM model can provide an approximate expression for the interference pattern produced by the above AI sequence. This can be done in two steps: (1) use the 1D trial wave function with $N_{c}=2$ from Eq. (20) to derive the condensate density at the end of the sequence, and (2) use approximate analytical solutions to the equations of motion given by Eqs. (33a) and (33) to obtain expressions for the variational parameters $\bar{x}_{1}$, $\dot{\bar{x}}_{1}$, $\bar{x}_{2}$, $\dot{\bar{x}}_{2}$, $\bar{w}_{1}$, $\dot{\bar{w}}_{1}$, $\bar{w}_{2}$, and $\dot{\bar{w}}_{2}$ at the end of the sequence.

The condensate trial wave function at the end of the initial split ($t=T_{\rm is}$) for two clouds ($N_{c}=2$) is given by Eq. (20) and can be written as follows:

[TABLE]

where we have used Eq. (32d). All variational parameters above are evaluated at $t=T_{\rm is}$. Since $\bar{\beta}_{j}=\dot{\bar{w}}_{j}/(4\bar{w}_{j})$ (from Eq. (32c) the wave function is expressed entirely in terms of the $\bar{x}_{j}$, $\bar{w}_{j}$ and their derivatives. These quantities can be found by solving Eqs. (33a) and (33).

We can model the effect of the second splitting by multiplying $\psi(\bar{x},\bar{T}_{\rm is})$ by a factor which splits each exisiting cloud into two clouds: one which is boosted by $v$ and one boosted by $-v$. Thus, in scaled units, the final–split wave function at $t=T_{\rm is}$ becomes

[TABLE]

The final–split wave function consists of four clouds. The velocities of these clouds just after the second splitting are $\approx 0$ for clouds labeled (1,1) and (2,2), $\approx+2v$ for cloud (1,2) and $\approx-2v$ for cloud (2,1).

Clouds (1,2) and (2,1) fly rapidly away from the center and after a short time, $\tau$, the overlap of the fast clouds with those at the center can be neglected. At this time in the AI sequence the density of the condensate can be imaged to obtain an interference pattern. To obtain a simple approximation for the condensate density at this time we assume that the change in the variational parameters during the time $\tau$ can be neglected and we approximate their values at time $t=T_{\rm is}+\tau$ with their values at $t=T_{\rm is}$.

The condensate density at the trap center can be found from the squared modulus of the sum of the wave functions of the two nearly motionless clouds. The result for the condensate density near the trap center after the final split is (in SI units)

[TABLE]

where

[TABLE]

is the final condensate phase difference and

[TABLE]

The quantity $k_{f}$ is the spatial frequency near the center of the interference pattern. We note that these variational quantities are understood to be evaluated at $t=T_{\rm is}$. We can now find an approximate interference pattern by solving Eqs. (33a) and (33) for the variational parameters either numerically or by finding approximate analytical solutions.

The presence of a source mass has two major effects on the interference pattern: (1) it creates interference fringes and (2) the peak of the pattern is shifted away from the zero–mass position towards the source–mass. We will call the width of the central fringe, $\lambda_{f}$, and the shift of the peak due to the source mass is called, $x_{c}$. Later we will show how $G$ can be obtained from the experimental parameters and either of these two characteristics of the interference pattern.

The final condensate phase, $\phi(x)$, predicts an interference pattern that displays equally spaced fringes near the center of the pattern due to the linear term. The fringe spacing narrows away from the center due to the quadratic term. The gradient of the phase measures the local relative velocity of the overlapping clouds.

This relative velocity is caused by (1) the differential accelerations of the two clouds during the wait period and (2) the changing widths of the two clouds due to condensate atom–atom interactions. Interactions also cause the fringe spacing to narrow away from the center of the interference pattern. The $\beta$ parameters in the LVM model describe how the clouds expand and/or contract due to the competition between the repulsive interactions and the harmonic confinement. They also account for any shift, due to atom–atom interactions, in the value of $G$ inferred from the interference pattern.

V.3 Extracting $G$ from the inteference pattern

In modeling a precision measurement the value of $G$ would be inferred from repeated GPE simulations of the experiment using different values of $G$ until a satisfactory match between the measured and simulated interference patterns was obtained. This is a necessary (but very expensive, see below) procedure. Such repeated simulations with different $G$ values can also be performed using the LVM model even though its accuracy is not sufficient to make a final decision about any particular AI scheme.

The LVM model can also help in designing an experiment by providing an approximate expression for $G$ in terms of the experimental parameters and characteristics of the interference pattern. This can be done by approximating the solution of the LVM equations of motion for each step of the proposed AI sequence. We now show this by deriving expressions for $G$ for the AI sequence described in Section IV.

To find approximate solutions of Eqs. (33a) and (33) we will make three assumptions: (1) the gravity effect on the condensate is small compared to the harmonic confinement so that $\omega_{SM}\ll\omega_{T}$ (see Eq. (42)), (2) the overlap terms, $X_{j}$ and $W_{j}$, can be neglected, and (3) the condensate width at the end of IS phase 2 can be approximated by the product of the IS phase 2 width expansion rate and the wait time, i.e., $\bar{w}(\bar{T}_{2})\approx\dot{\bar{w}}(\bar{T}_{2})\bar{T}_{W}$.

One immediate consequence of this is that we can neglect $\omega_{SM}$ in the equations of motion for the $\bar{w}_{j}$ and so they are approximately unaffected by gravity during the entire AI sequence. With these approximations we can also infer that the time evolution of the widths of the two clouds are approximately the same, i.e., $\bar{w}_{1}(\bar{t})\approx\bar{w}_{2}(\bar{t})$ and also $\dot{\bar{w}}_{1}(\bar{t})\approx\dot{\bar{w}}_{2}(\bar{t})$. This also implies that $\bar{\beta}_{1}(\bar{t})\approx\bar{\beta}_{2}(\bar{t})$.

These approximations also enable us to find analytical solutions of Eq. (33a) for the cloud center positions in each of the three phases of the initial split (panels (c), (d), and (e) of Fig. 1, respectively). In phases 1 and 3 the solutions are $\sin(\bar{\omega}_{T}\bar{t})$ and $\cos(\bar{\omega}_{T}\bar{t})$ and, in phase 2, the solutions are $\sinh(\bar{\omega}_{SM}\bar{t})$ and $\cosh(\bar{\omega}_{SM}\bar{t})$. Using the initial conditions $\bar{x}_{1}(0)=\bar{x}_{2}(0)=0$, $\dot{\bar{x}}_{1}(0)=-\bar{v}$, and $\dot{\bar{x}}_{2}(0)=+\bar{v}$ we can solve equations of motion for the center coordinates in each IS phase using the final values of $\bar{x}_{j}$ and $\dot{\bar{x}}_{j}$ in the previous phase as the initial values in the next phase.

This straightforward procedure enables us to obtain approximate expressions for the center positions and velocities at $t=T_{is}$. The results are

[TABLE]

where $j=1,2$, $\epsilon_{1}=-1$, $\epsilon_{2}=+1$, $\phi_{SM}\equiv\sqrt{2}\omega_{SM}T_{W}$ and from Eq. (94)

[TABLE]

It is also possible to find an approximate value for $\bar{\beta}_{1}(\bar{T}_{is})$ (assumed to be equal to $\bar{\beta}_{2}(\bar{T}_{is})$) by approximating the width equation of motion (Eq. (33)) during phase 3 of the IS. An example of how the width behaves during this time is shown in Fig. 3 between the two rightmost vertical dotted lines. The cloud width is initially quite large and drops rapidly to a value close to its starting value.

Thus during the IS phase 3 period all of the interaction terms in Eq. (33) can be neglected and thus the equation of motion can be approximated as (for cloud 1)

[TABLE]

with initial conditions $\bar{w}_{1}(\bar{T}_{2})$ and $\dot{\bar{w}}_{1}(\bar{T}_{2})$ at $\bar{t}=\bar{T}_{2}$. This can be easily solved and, since the duration of phase 3 is one–quarter of the trap period, we can write $\bar{w}_{1}(\bar{T}_{is})$ and $\dot{\bar{w}}_{1}(\bar{T}_{is})$ in terms of $\bar{w}_{1}(\bar{T}_{2})$ and $\dot{\bar{w}}_{1}(\bar{T}_{2})$ as follows

[TABLE]

and

[TABLE]

where in the last equality we have used the approximation mentioned above that the cloud width at the end of IS phase 2 is approximately the width velocity at $t=T_{2}$ multiplied by the wait time, $T_{W}$. So, finally we have an approximate expression for the $\beta$ parameters at $t=T_{is}$:

[TABLE]

This expression enables us to derive an expression for $G$ solely in terms of the experimental parameters and characteristics of the measured interference pattern.

We can use the expressions for the $\bar{x}_{j}(\bar{T}_{is})$ and the $\dot{\bar{x}}_{j}(\bar{T}_{is})$ from Eqs. (V.3)) and the $\bar{\beta}_{j}(\bar{T}_{is})$ in Eq. (46) to approximate $k_{f}$ in Eq. (40). This equation can be written as (in SI units)

[TABLE]

The left–hand–side (LHS) of this equation contains only experimental parameters ($M$ and $v$) and a characteristic of the interference pattern ($\lambda_{f}$ is the width of the central fringe). The value of the LHS can thus be used to find the value of $\phi_{SM}=\sqrt{2}\omega_{SM}T_{W}$ from which $G$ can be obtained. Thus this expression can be used to estimate the value of $G$ from the experimental parameters and the interference pattern.

This is most easily seen if we assume that $\phi_{SM}\ll 1$ in which case the above expression for $G$ can be further approximated as $k_{f}\approx Mv\phi_{SM}^{2}/2\hbar$. So, in SI units, we can write

[TABLE]

where $M$ is the mass of a condensate atom and $\lambda_{f}$ is the width of the central fringe of the interference pattern.

This formula expresses $G$ in terms of the experimental parameters, $T_{W}$, $x_{SM}$, $M_{SM}$, and $v$, and, $\lambda_{f}$, the width of the central fringe of the interference pattern. This quantity is easily inferred from the pattern. Interference patterns for several sets of experimental parameters are shown in Fig. 4. In particular, the figure shows how the patterns change as the wait time, $T_{W}$, is increased. The figure also shows the values of $G$ given by the above formula for each pattern.

It is also possible to write an expression for $G$ in terms of the shift, $x_{c}$, of the interference pattern maximum away from the maximum of the pattern obtained when the source mass is absent. The result is

[TABLE]

Equations (48) and (49) display the ability of LVM model to show how $G$ can be obtained from the experimental parameters and the data. They can further be used to estimate the sensitivity of the measured $G$ value to each of the experimental parameters and to develop a preliminary error budget for the experiment.

V.4 Using the LVM model for preliminary AI design

One of the features of the LVM model is that different AI schemes can be rapidly evaluated and compared. The scheme presented above only splits the condensate into two fast pieces. Other possible AI schemes might, for example, split the condensate into three, five or more pieces (see Ref. Edwards et al. (2010)). Whether such ideas are useful can only be determined by modeling the particular AI sequence envisioned. In general this will require many iterations in order to come to a conclusion about the candidate AI sequence. We have also shown above that the model is also useful in shrinking the experimental parameter space.

If this modeling is conducted by numerical solution of the GPE, this will require many runs each of which is very expensive. For example, in the simulation shown in Fig. 2 the GPE run needed $10^{4}$ seconds to finish using a fast desktop computer while the LVM simulation required $<1$ second to finish. This simulation where $T_{W}=0.5$ s and $v=3.15\times 10^{-3}$ m/s was very different from the much more extreme conditions found in Fig. 4(c) where $T_{W}=20$ s and $v=9.14\times 10^{-2}$ m/s. The size of the grid box in the second simulation was 30 times larger and the evolution time was 40 times longer making the full simulation take 1200 times longer. Thus, if the same GPE simulator used for the first simulation were used to perform the second simulation, it would have taken the desktop computer more than $10^{7}$ seconds (or about 4 months) to finish this single simulation. The LVM simulation took less than three minutes. If dozens or hundreds of preliminary modeling runs were needed, it would clearly not be practical to use the GPE even if such modeling were done on a supercomputer. Furthermore, a realistic AI experiment for measuring $G$ would take place in three dimensions and would require 3D modeling. This would make simulation time for the GPE even longer but would not be significantly longer for the LVM model.

We believe that it is not necessary to use the GPE for preliminary design modeling of candidate AI sequences. This can be accomplished with the LVM and precision modeling can carried out using the GPE only for a small number of likely AI candidate schemes. A number of such likely 3D schemes will be the subject of a future publication.

VI Summary

In this paper we presented a variational technique for approximating the solution of the time–dependent Gross–Pitaevskii equation, for situations where the condensate is split into multiple parts as is common in atom interferometry processes. The technique presented can approximate the solution of either the one–dimensional GPE or the three–dimensional GPE in the laboratory or rotating frame.

Each part of the condensate is modeled as a Gaussian with time–dependent center, width, and linear and quadratic phase parameters that enable the centers and widths to change. Different parts of the condensate are allowed to overlap but exchange of atoms between different condensate parts is assumed to be negligible. The resulting equations of motion for the center coordinates, widths, and phase parameters is a system of ordinary differential equations in time that can be rapidly solved.

The power of this model is that this set of equations is cast in terms of derivatives of an arbitrary external potential and thus can be applied to many different AI schemes. This model should enable rapid prototyping and design of novel AI schemes for applications especially in microgravity environments.

Acknowledgements.

This material is based upon work supported by the U.S. National Science Foundation under grant numbers PHY–1004975, PHY–0758111, the Physics Frontier Center grant PHY–0822671 and by the National Institute of Standards and Technology.

Appendix A 3D LVM Lagrangian

Here we sketch the derivation of the 3D Lagrangian. This Lagrangian has five terms, defined in Eq. (13), and we will derive $\bar{L}^{(3D)}_{1}$, $\bar{L}^{(3D)}_{2}$ and $\bar{L}^{(3D)}_{5}$ explicitly. Further we will show that $\bar{L}^{(3D)}_{3}$ and $\bar{L}^{(3D)}_{4}$ only depend on the $\bm{x}$ and $\bm{w}$ variational parameters. This in turn motivates the introduction of the “variational” potentials, $\bar{U}_{\rm ext}(\bm{x},\bm{w})$ and $\bar{U}_{\rm int}(\bm{x},\bm{w})$, which later appear in the equations of motion for the Gaussian center coordinates and widths. The explicit expression for $\bar{U}_{\rm int}(\bm{x},\bm{w})$ is given in Appendix B. The form of $\bar{U}_{\rm ext}(\bm{x},\bm{w})$ for the case of a 1D harmonic trap plus the gravitational potential of an arbitrarily placed point mass will be given in Appendix E.

A.1 Derivation of $\bar{L}^{(3D)}_{1}$

The $\bar{L}_{1}$ term of the Lagrangian has the form

[TABLE]

The trial wave function and its time derivative are

[TABLE]

where the dot denotes partial differentiation with respect to $\bar{t}$. Inserting these into Eq. (50), neglecting the rapidly oscillating terms and then taking the imaginary part yields the following result

[TABLE]

where in the second line we have used Eqs. (23) to eliminate the $A_{j}^{2}(\bar{t})$ factors. The Gaussian integrals appearing in this last expression can easily be evaluated. The final result for $\bar{L}^{(3D)}_{1}$ is $\bar{L}_{1}$:

[TABLE]

Next we derive $\bar{L}^{(3D)}_{2}$.

A.2 Derivation of $\bar{L}^{(3D)}_{2}$

The expression for $\bar{L}^{(3D)}_{2}$ is given by

[TABLE]

If we differentiate $\Psi$ with respect to $\bar{\eta}$ ($\eta$ is $x$, $y$, or $z$), multiply the result with its complex conjugate, and neglect the rapidly oscillating terms (which integrate to zero), we obtain

[TABLE]

Inserting this into Eq. (54) gives

[TABLE]

Again we are left with Gaussian integrals that can be straightforwardly evaluated. The final result is

[TABLE]

In the above we have, again, used Eqs. (23) to eliminate the $A_{j}^{2}(\bar{t})$ factors.

A.3 3D Variational Potentials

Next we consider the terms $L^{(3D)}_{3}$ and $L^{(3D)}_{4}$ which are

[TABLE]

Note that, in terms of the trial wave function, these depend on $\left|\Psi({\bf\bar{r}},\bar{t})\right|^{2}$ or its square. If we form $\Psi^{\ast}\Psi$ with our trial wave function and neglect the rapidly oscillating terms, we obtain

[TABLE]

It is clear that $|\Psi|^{2}$ only depends on the centers, $\bm{x}$, and widths, $\bm{w}$, if rapidly oscillating terms can be neglected. This is evidently also the case for $|\Psi|^{4}$ which is the square of $|\Psi|^{2}$. Thus, clearly, we have $\bar{L}^{(3D)}_{3}=\bar{L}^{(3D)}_{3}(\bm{x},\bm{w})$ and $\bar{L}^{(3D)}_{4}=\bar{L}^{(3D)}_{4}(\bm{x},\bm{w})$, that is, these quantities only depend on the Gaussian center and width variational parameters.

Thus we introduce here the following variational potentials. First we define the “external” variational potential that accounts for the effects of $V_{\rm ext}$:

[TABLE]

We further introduce the “interaction” variational potential that accounts for atom–atom scattering interactions:

[TABLE]

Note that $\bar{U}^{(3D)}_{\rm ext}(\bm{x},\bm{w})$ depends on the particular form of $V_{\rm ext}({\bf r},t)$. The expression for $\bar{U}^{(3D)}_{\rm int}(\bm{x},\bm{w})$ is always the same. The derivation of $\bar{U}^{(3D)}_{\rm int}(\bm{x},\bm{w})$ is more complex than that of the other Lagrangian terms and is given in Appendix B. It is only necessary to know that these potentials depend on $\bm{x}$ and $\bm{w}$ in order to derive the equations of motion for the Gaussian centers and widths. The final equations of motion will be written in terms of the gradients of $\bar{U}^{(3D)}_{\rm ext}(\bm{x},\bm{w})$ and $\bar{U}^{(3D)}_{\rm int}(\bm{x},\bm{w})$.

A.4 Derivation of $\bar{L}^{(3D)}_{5}$

The term $\bar{L}^{(3D)}_{5}$ accounts for the possibility that the system is in a rotating frame. Here we will assume that there is rotation only about the $z$ axis. The $\bar{L}^{(3D)}_{5}$ term of the Lagrangian then has the form

[TABLE]

This term can be straightforwardly evaluated by inserting the expressions for the trial wave function and its space derivative into Eq. (63), neglecting the oscillating terms, evaluating the resulting Gaussian integrals, and eliminating the $A_{j}^{2}(\bar{t})$ factors. The result is

[TABLE]

A.5 The 3D Lagrangian

Combining Eqs. (53), (57), (26), (27), and (64) we have the full 3D Lagrangian

[TABLE]

We will use this form of the Lagrangian to derive the 3D equations of motion.

Appendix B Derivation of $\bar{U}_{int}^{(3D)}({\bf x},{\bf w})$ and $\bar{U}_{int}^{(1D)}({\bf x},{\bf w})$

Here we present the derivation of the expressions for $\bar{U}^{(3D)}_{\rm int}({\bf x},{\bf w})$ and $\bar{U}^{(1D)}_{\rm int}({\bf x},{\bf w})$. The expression for $\bar{U}^{(3D)}_{\rm int}({\bf x},{\bf w})$ is

[TABLE]

In order to perform this integral we must first calculate $|\Psi|^{4}$. We can write this quantity as the square of $|\Psi|^{2}$:

[TABLE]

where the have written the $|\Psi|^{2}$ appearing inside the large square brackets as a sum of two terms. The first term is a sum over terms that definitely do not oscillate in space ($j_{1}=j_{2}$) and the second term, where $j_{1}\neq j_{2}$, all of the terms definitely do oscillate in space.

When this expression is expanded by writing the terms of the overall square only the square of the first term and the square of the last term need be retained. We only need to integrate terms that definitely do not oscillate as the oscillatory terms will integrate to zero. Thus we can rewrite the above expression as follows:

[TABLE]

The second term above still has oscillating terms. We only want to keep the non–oscillating terms. The terms that don’t oscillate are those where $j_{1}^{\prime}=j_{2}$ and $j_{2}^{\prime}=j_{1}$ (since $j_{1}=j_{2}$ and $j_{1}^{\prime}=j_{2}^{\prime}$ are excluded already). Thus we can evaluate the primed sums keeping only those terms where $j_{1}^{\prime}=j_{2}$ and $j_{2}^{\prime}=j_{1}$:

[TABLE]

The two double sums appearing in the expression above are the same with one exception: the second double sum requires $j_{1}\neq j_{2}$ while the first double sum has no such restriction. We can write the first sum as two terms: a single sum where $j_{1}=j_{2}$ and a double sum where $j_{1}\neq j_{2}$. This second term will then be identical to the second term in the original expression. Carrying out this procedure writing $|\Psi|^{4}$ in terms of coordinates yields

[TABLE]

This form for $|\Psi|^{4}$ can now be inserted into Eq. (66) and then resulting Gaussian integrals can be straightforwardly evaluated. The final result is

[TABLE]

Using a similar derivation, we obtain the 1D version of $\bar{U}_{\rm int}$:

[TABLE]

The derivatives of these expressions with respect to the center and width variational parameters are computed straightforwardly.

Appendix C The LVM Equations of Motion

In this section we derive the LVM equations of motion for the variational parameters appearing in the 3D and 1D trial wave functions given above. The major steps in this undertaking are the following. First we write down the equation of motion for each variational parameter, $\bar{q}$, using the standard Euler–Lagrange (EL) equation (Eq. (6)). We will show that further manipulations of these equations can produce a closed set of second–order differential equations in time for the center position and width coordinates. Furthermore we will show that all of the other variational parameters can be expressed in terms of the position and width coordinates and their time derivatives. We begin with the 3D case.

C.1 3D LVM equations of motion

C.1.1 3D cloud center EOMs

To obtain the equations of motion for the cloud–center coordinates we need the EOMs associated with the $\bar{\alpha}_{j\eta}$ and the $\bar{\eta}_{j}$. We begin with the EOMs associated with the $\bar{\alpha}_{j\eta}$. The Euler–Lagrange equations for these are

[TABLE]

We can compute the derivatives of $\bar{L}$ using Eq. (25). The straightforward result is

[TABLE]

Next we obtain the EL equations associated with $\bar{x}_{j}$, $\bar{y}_{j}$, and $\bar{z}_{j}$. The $\bar{\eta}_{j}$ Euler–Lagrange EOMs can be written as

[TABLE]

where the second equality holds as there are no $\dot{\bar{\eta}}$ terms in the Lagrangian. Taking the necessary derivatives we obtain the EL equations associated with $\bar{j}$, $\bar{y}_{j}$, and $\bar{z}_{j}$, respectively

[TABLE]

We have written the above equations in a form for convenient use in deriving the final equations of motion. Equations (70), (71), and (72) can be used along with Eqs. (74), (75), and (76) to derive second–order differential equations that involve only the center and width coordinates.

We illustrate how this can be done by deriving the equation for $\ddot{\bar{x}}_{j}$. If we differentiate both sides of the equation for $\dot{\bar{x}}_{j}$ (Eq. (70)) with respect to time, the resulting second–order equation will contain both $\ddot{\bar{x}}_{j}$ and $\dot{\bar{x}}_{j}$. Using Eq. (70) again to eliminate $\dot{\bar{x}}_{j}$ from this second–order equation gives the following result.

[TABLE]

where we have obtained the second line by noting that the quantity in square brackets appearing in the first line is identical to the quantity in square brackets in Eq. (74). The second equation can be further simplified by noting that the quantity in the parenthesis appearing in the second term is also present in Eq. (70). The second–order equation of motion is thus written in the compact final form

[TABLE]

The corresponding equations for the other cloud–center coordinates can be derived similarly.

The resulting second–order equations of motion for the center coordinates of cloud $j$ are as follows.

[TABLE]

Note that these equations depend only on ${\bf x}$, ${\bf w}$ and their time derivatives. We now turn to the equations for the cloud widths.

C.1.2 3D cloud width EOMs

It is possible to derive a set of second–order equations of motion for the Gaussian cloud widths similar to that for the cloud centers. As will be seen below, the center equations and width equations form a closed system. All of the other variational parameters can be expressed in terms of the centers and widths and their time derivatives.

To obtain the width equations we start with the EL equation for the $\bar{\beta}_{j\eta}$. As before, we will only derive the equation of motion for $\bar{w}_{jx}$ to illustrate how the derivation is carried out. The EL equation for $\bar{\beta}_{jx}$ is

[TABLE]

Differentiating the 3D Lagrangian in Eq. (25) and inserting into the above equations gives

[TABLE]

If we replace $\dot{\bar{x}}_{j}$ on the left–hand–side with its expression given in Eq. (70) we obtain the following remarkably simple result

[TABLE]

This result holds for $\bar{y}$ and $\bar{z}$ as well and enables us to express the $\bar{\beta}_{j\eta}$ in terms of the widths.

To proceed we turn to the Euler–Lagrange equation for $\bar{w}_{jx}$ which reads

[TABLE]

where the last equality holds because no terms containing $\dot{\bar{w}}_{jx}$ appear in the Lagrangian given in equation Eq. (25). Again taking derivatives of $\bar{L}^{(3D)}$ we obtain

[TABLE]

Now note that, if we differentiate both sides of Eq. (82) with respect to time we obtain

[TABLE]

where the second equality comes from reusing Eq. (82) to replace $\dot{\bar{w}}_{jx}$ in the first line. Finally note that the right–hand–side of this last expression is identical to the left–hand–side of Eq. (84) and thus we obtain the second–order equation for $\bar{w}_{jx}$:

[TABLE]

The derivation of equations for $\bar{y}$ and $\bar{z}$ are similar. All three equations can be written in the following compact form:

[TABLE]

We can add these to EOMs for the cloud centers to get the full set of equations of motion in 3D. They consist of a pair of second–order ordinary differential equation for the cloud centers and widths as well as expressions for the $\bar{\beta}_{j\eta}$ and the $\bar{\alpha}_{j\eta}$ in terms of the centers, widths and their first derivatives:

[TABLE]

The equations for the cloud centers and cloud widths (Eqs. (88a), (88b), (88c), and (88d)) form a closed set that contain only the $\bar{\eta}_{j}$, $\dot{\bar{\eta}}_{j}$, $\bar{w}_{j\eta}$, and $\dot{\bar{w}}_{j\eta}$. Once these quantities are obtained, all of the other variational parameters can be calculated.

Appendix D Details of the GPE simulation

The conditions assumed in the simulation shown in Fig. 2 include a condensate of 10,000 87Rb atoms which were confined in a harmonic trap with a frequency of $\omega_{T}/2\pi=1$ Hz. The condensate pieces were given an initial velocity of $v=2\times 10^{-3}$ m/s and the wait time after reaching the trap turning points is $T=500$ ms. The source mass is absent in the simulation presented.

The GPE was solved numerically using the split–step, Crank–Nicolson algorithm Muruganandam and Adhikari (2009). Briefly, in this algorithm, the Hamiltonian is split as $\hat{H}=\hat{T}+\hat{V}$ where $\hat{T}$ is the kinetic energy and $\hat{V}$ is the potential plus the nonlinear interaction term. The condensate wave function is advanced from time $t$ to time $t+\delta t$ by multiplying $\psi(x,t)$ point by $e^{-iV(x,t)\delta t/\hbar}$ and then the result is multiplied by $e^{-i\hat{T}\delta t/\hbar}$. This second step is carried out using the Crank–Nicolson algorithm. The time–dependent variational equations of motion were solved by the Euler method.

The comparison of the GPE and the LVM simulations for the initial split phase of the AI sequence is shown in Fig. 2. The initial split itself divided into three parts. In the first part the condensate pieces fly out to the turning points in the presence of the harmonic potential whose frequency is $\omega_{T}/2\pi=1$ Hz and finally stop after a quarter period (250 ms) as shown in Figs. 1(b)–(c). Both the GPE and the LVM densities are plotted together on the same graph at four different times as seen in the left panel of Fig. 2. Each density consists of two symmetrically placed peaks and the left peak of each density is labeled with its time stamp. In the second part of the initial split sequence the two condensate pieces evolve at rest for a wait time of $T=500$ ms while the trap is off. This evolution is shown in the middle panel of Fig. 2. Finally the condensate pieces recombine after the trap is turned back on. This is shown in the right panel of Fig. 2.

Appendix E Derivation of $\bar{U}_{ext}^{(1D)}({\bf x},{\bf w})$ for an harmonic plus point–mass potential

In this appendix we derive the variational potentials for $\bar{U}_{ext}^{(1D)}({\bf x},{\bf w})$ for the case of an external potential consisting of an harmonic potential plus the gravatational potential of a point mass situated far from the condensate. The 1D version is used to obtain the equations of motion whose solutions are compared with the full 1D GPE simulation of the illustrative AI measurement of $G$ in a microgravity environment presented in Section IV.

The 1D external variational potential is given by Eq. (28)

[TABLE]

The first step is to find an approximate expression for $\bar{V}_{\rm ext}(\bar{x},\bar{t})$ for the combination of an harmonic potential plus a point mass.

We assume that the harmonic trap is centered at the origin of coordinates and that a point mass with mass $M_{SM}$ is located at ${x}_{SM}$. The exact potential can thus be written (in SI units) as

[TABLE]

We want to approximate $V_{\rm G}({\bf r})$ by assuming that the distance of any cloud to the origin is much smaller than the distance of the source mass to the origin. We have chosen the origin to be at the center of the harmonic potential confining the BEC.

First we consider only the gravitational part of the potential:

[TABLE]

we can approximate this exact expression by making a Taylor expansion about $x=0$ to second order in $x/x_{SM}$

[TABLE]

This expression is valid only for points $x$ such that $x_{SM}>|x|$ which we take to be the case since we are actually assuming $x_{SM}\gg|x|$. It will be convenient here to introduce the gravitational frequency

[TABLE]

Thus we can rewrite the approximate gravitational potential in terms of this quantity as follows.

[TABLE]

where, in the second line, we have expressed the gravitational potential in scaled units. Now we are ready to write down the full external potential. The full external potential in scaled units is thus

[TABLE]

Inserting the above expression into Eq. (LABEL:Uext_1d_appB) yields the final expression for $\bar{U}^{(1D)}_{\rm ext}$ for the case of an external harmonic trap plus the gravitational potential produced by a point mass far away from the condensate is given by

[TABLE]

We note that the last term is constant and never appears in the equations of motion.