Tools for designing atom interferometers in a microgravity environment
Elizabeth Ashwood, Ed Wesley Wells, Doga Murat Kurkcuoglu, Robert Colson Sapp, Charles W Clark, Mark Edwards
aa r X i v : . [ c ond - m a t . qu a n t - g a s ] M a r Tools for designing atom interferometers in a microgravity environment
Elizabeth Ashwood, Ed Wesley Wells, Doga Murat Kurkcuoglu, Robert Colson Sapp, Charles W. Clark, and Mark Edwards
1, 2 Department of Physics, Georgia Southern University, Statesboro, GA 30460–8031 USA Joint Quantum Institute, National Institute of Standards and Technologyand the University of Maryland, Gaithersburg, MD 20899, USA (Dated: March 12, 2019)We present a variational model suitable for rapid preliminary design of atom interferometers in amicrogravity environment. The model approximates the solution of the 3D rotating–frame Gross–Pitaevskii equation (GPE) as the sum of N c Gaussian clouds. Each Gaussian cloud is assumedto have time–dependent center positions, widths, and linear and quadratic phase parameters. Weapplied the Lagrangian Variational Method (LVM) with this trial wave function to derive equationsof motion for these parameters that can be adapted to any external potential . We also present a 1Dversion of this variational model. As an example we apply the model to a 1D atom interferometryscheme for measuring Newton’s gravitational constant, G , in a microgravity environment. We showhow the LVM model can (1) constrain the experimental parameter space size, (2) show how the valueof 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 errorbudget. PACS numbers: 03.75.Gg,67.85.Hj,03.67.Dg
I. INTRODUCTION
Atom interferometers (AIs) have been used in manyapplications [1, 2]. In particular this is the case forgravitational measurements, such as the determinationof Earth’s gravitational acceleration, g , [3] and its gradi-ents [4, 5]. Atom interferometers have also been proposedfor use in Sagnac gyroscopes [6, 7], in testing general rel-ativity [8–10], in searching for dark energy [11], and mea-suring gravity waves [12, 13]. The Newtonian constant ofgravitation, G , has also been measured with AIs [14, 15].In a typical atom interferometer involving a gaseousBose–Einstein condensate (BEC), the basic sequenceconsists of splitting the condensate into two or moreclouds that then experience different environments whileseparated. Some differing property of these environmentsgenerates a phase difference between the clouds whilethey are separated. The phase difference produces an in-terference pattern when the clouds are reunited and thispattern can be used to infer the value of this property.An emerging area of AI applications is to systems thatare freely falling in a gravitational field. Such systems of-fer the prospect of much longer interrogation times thanare possible in fixed terrestrial laboratories, where atomsare subject to constant vertical acceleration. Such “mi-crogravity” environments have been created in drop tow-ers [16], sounding rockets [17], and recently aboard theInternational Space Station (ISS).The National Aeronautics and Space Administra-tion (NASA) has deployed its Cold Atom Labora-tory (CAL) [18] to the ISS enabling future atom–interferometry experiments to be performed on Bose–Einstein condensates there. This will provide a user fa-cility in which atom–interferometry experiments can beconducted in a microgravity environment. The basic AI sequence of splitting, reuniting, and split-ting the BEC again to produce an interference patterncan be viewed as a sequence of sudden splits of the con-densate alternating with periods where the condensateclouds evolve freely according to their local environments.In general, the performance of AIs of this type dependson several factors including the number of interferometerpathways, the area of the interferometer, and the inter-rogation time (the time the clouds spend in the differentenvironments).AI performance can be enhanced both by increasingany or all of these quantities. Specific performance en-hancements can also be achieved by designing differentAI sequences as was done recently with neutron inter-ferometers [19]. Differential AI experiments, where thesame AI sequence is carried out on identical condensatesat the same time, can also be used in making precisionmeasurements.Enhancing AI performance by increasing interrogationtime is especially possible in a microgravity environment.The extended interrogation times available in micrograv-ity are also valuable for precision measurement applica-tions. Microgravity environments in general and the CALin particular are important platforms for making preci-sion measurements. To further this end, designing newAI sequences that take advantage of microgravity wouldbe highly desirable.Designing new AI sequences involving BECs suitablefor precision measurements conducted in microgravityenvironments presents a demanding numerical challenge.Detailed modeling of condensate behavior for AI se-quences having long interrogation times and/or multiplepathways requires large amounts of computer time andstorage assuming that the condensate obeys the Gross-Pitaevskii equation (GPE). Many different possible se-quences would also need to be simulated and the resultsevaluated. Thus a method for rapid simulation of a givenAI sequence would be an extremely useful tool for AI de-sign. Such a tool would also be helpful in determininghow the quantity measured by the AI could be extractedfrom the resulting interference patterns.This paper presents a tool suitable for preliminary de-sign of candidate AI sequences. The tool can providerapid approximate modeling of condensate behavior be-tween the splits of an AI sequence where this evolution isassumed to be governed by the GPE. This model is basedon the Lagrangian Variational Method (LVM) where theassumed trial wave function approximates the GPE so-lution as a set of N c separate but possibly overlappingGaussian clouds. We emphasize that the final choice ofan AI sequence would have to be validated by simulatingit with the GPE.In section II we describe the general Lagrangian Vari-ational Method, introduce a set of scaled units suitablefor computation, and present the trial wave functions forboth the one-dimensional and three-dimensional versionsof the method. In section III we derive the equationsof motion for the variational parameters and show thatthese can be cast in terms of the space and width gradi-ents of a “variational potential” which is just the expec-tation value of the external potential over the trial wavefunction.In section IV we describe an example of an AI schemethat could be used to measure G in microgravity. In sec-tion V we apply the LVM model to this AI scheme. Wecompare a GPE simulation of this scheme with the LVMsimulation. We also use the model to obtain an expres-sion for the LVM interference pattern, show how G canbe extracted from this pattern, and how the model canfacilitate preliminary design of an AI sequence. Finally,we present a summary of the work in section VI. II. THE LAGRANGIAN VARIATIONALMETHOD AND TRIAL WAVE FUNCTIONA. The general Lagrangian variational method
The Lagrangian Variational Method provides approx-imations to the solutions of the time–dependent Gross–Pitaevskii equation [20–24]. In three dimensions and inthe rotating frame, this equation has the form [25]: i ¯ h ∂ Φ ∂t = − ¯ h M ∇ Φ + V ext ( r , t )Φ + g D N | Φ | Φ+ i ¯ h Ω · ( r × ∇ ) Φ , (1)where Φ( r , t ) is the condensate wave function, M is themass of a condensate atom, N is the number of atoms inthe condensate, g D = 4 π ¯ h a s /M measures the strengthof the atom–atom scattering with a s being the scatteringlength, V ext ( r , t ) is the potential exerted on a condensateatom by external fields, and Ω is the angular velocity ofthe rotating frame. The LVM is based on the fact that the GPE can bederived as the Euler–Lagrange equation of motion pro-duced by the following Lagrangian density: L [Φ ∗ ] = ¯ h Im { Φ ∗ Φ t } + ¯ h M ∇ Φ ∗ · ∇ Φ+ V ext ( r , t )Φ ∗ Φ + g D N (Φ) (Φ ∗ ) + Φ (cid:16) Ω · ˆ L (cid:17) Φ ∗ , (2)where ˆ L is the single–particle angular momentum oper-ator.This Lagrangian density along with the followingEuler–Lagrange equation of motion produces the GPE: X η = x,y,z,t ∂∂η (cid:18) ∂ L ∂ Φ ∗ η (cid:19) − ∂ L ∂ Φ ∗ = 0 , where Φ η ≡ ∂ Φ ∂η . (3)The Lagrangian Variational Method consists of devis-ing a trial wave function,Φ trial ( r , t ) = Φ trial ( q ( t ) , . . . , q n ( t ); r ) (4)where the { q i ( t ) } , i = 1 , . . . , n are variational parameters that only depend on the time, t . The equations of motionof these variational parameters are derived by computingthe ordinary Lagrangian: L ( q ( t ) , . . . , q n ( t )) = Z d r L [ (cid:0) Φ trial (cid:1) ∗ ] (5)and then using the standard Euler–Lagrange equations, ddt (cid:18) ∂L∂ ˙ q k (cid:19) − ∂L∂q k = 0 , k = 1 , . . . , n (6)to produce an equation of motion associated with eachvariational parameter. B. Scaled units
We can simplify the equations produced by the abovemethod by introducing a set of units appropriate to theproblem and a set of scaled variables (both independentand dependent). The scaled variables are defined by firstchoosing a length unit , L , and then defining energy andtime units in terms of L : E ≡ ¯ h M L and T ≡ ¯ hE = 2 M L ¯ h . (7)We then introduce scaled variables which are generallydenoted by barred quantities. These consist of scaledspace and time coordinates:¯ x ≡ xL ¯ y ≡ yL ¯ z ≡ zL and ¯ t ≡ tT . (8)We also introduce the scaled condensate wave functionfor the solution of the GPE:Φ( r , t ) = L − / Ψ( ¯r , ¯ t ) . (9)We can express the original GPE in terms of scaled quan-tities and this can be done for the Lagrangian density andits associated Euler–Lagrange equation as well.In terms of scaled quantities, the rotating–frame GPEbecomes: i ∂ Ψ ∂ ¯ t = − ¯ ∇ Ψ + ¯ V ext ( ¯r , ¯ t )Ψ + ¯ g D N | Ψ | Ψ+ i ¯Ω z (cid:18) ¯ x ∂ Ψ ∂ ¯ y − ¯ y ∂ Ψ ∂ ¯ x (cid:19) (10)where we are considering the form of the rotating–frameGPE for the special case of rotation around the z axis.In the above we have g D ≡ ¯ g D E L , V ext ( ¯r , ¯ t ) =¯ V ext ( r , t ) E and Ω z = ¯Ω z /T . The scaled Lagrangiandensity for this version of the GPE becomes¯ L (3 D ) [Ψ ∗ ] = Im { Ψ ∗ Ψ ¯ t } + ¯ ∇ Ψ ∗ · ¯ ∇ Ψ + ¯ V ext ( ¯r , ¯ t )ΨΨ ∗ + ¯ gN (Ψ ∗ ) (Ψ) + i ¯Ω z Ψ (cid:0) ¯ y Ψ ∗ x − ¯ x Ψ ∗ y (cid:1) ≡ ¯ L (3 D )1 [Ψ ∗ ] + ¯ L (3 D )2 [Ψ ∗ ] + ¯ L (3 D )3 [Ψ ∗ ]+ ¯ L (3 D )4 [Ψ ∗ ] + ¯ L (3 D )5 [Ψ ∗ ] (11)and the scaled Euler–Lagrange equation is given by X η = x,y,z,t ∂∂ ¯ η (cid:18) ∂ ¯ L (3 D ) ∂ Ψ ∗ ¯ η (cid:19) − ∂ ¯ L (3 D ) ∂ Ψ ∗ = 0 . (12)It is straightforward to insert Eq. (11) into Eq. (12) andobtain Eq. (10). The 3D Lagrangian is computed byintegrating over all 3D space:¯ L (3 D ) [ x , w , α , β ] = Z d ¯ r ¯ L (3 D ) [Ψ ∗ ]= X k =1 Z d ¯ r ¯ L (3 D ) k ≡ ¯ L (3 D )1 + ¯ L (3 D )2 + ¯ L (3 D )3 + ¯ L (3 D )4 + ¯ L (3 D )5 . (13)In the above have denoted the dependence of ¯ L on the dif-ferent types of variational parameters that will appear inthe trial wave function defined in the next section. Theseare: the cloud center coordinates by x ≡ (¯ x , . . . , ¯ z N c ),the cloud widths denoted by w ≡ ( ¯ w x , . . . , ¯ w N c z ), thelinear phase coefficients by α ≡ ( α x , . . . , α N c z ) and thequadratic phase coefficients by β ≡ ( β x , . . . , β N c z ).The LVM formulation is straightforward to apply inone–dimension where we denote the 1D GPE solutionby φ ( x, t ). For scaled units we will write φ ( x, t ) ≡ L − / ψ (¯ x, ¯ t ). In terms of scaled quantities, the 1D GPEbecomes: i ∂ψ∂ ¯ t = − ∂ ψ∂ ¯ x + ¯ V ext (¯ x, ¯ t ) ψ + ¯ g D N | ψ | ψ. (14) where g D ≡ ¯ g D E L and V ext (¯ x, ¯ t ) = ¯ V ext ( x, t ) /E .The 1D scaled Lagrangian density becomes¯ L (1 D ) [ ψ ∗ ] = Im { ψ ∗ ψ ¯ t } + ψ ∗ ¯ x ψ + ¯ V ext (¯ x, ¯ t ) ψψ ∗ + ¯ g D N ( ψ ∗ ) ( ψ ) (15)and the 1D scaled Euler–Lagrange equation is given by ∂∂ ¯ x (cid:18) ∂ ¯ L (1 D ) ∂ψ ∗ ¯ x (cid:19) + ∂∂ ¯ t (cid:18) ∂ ¯ L (1 D ) ∂ψ ∗ ¯ t (cid:19) − ∂ ¯ L (1 D ) ∂ψ ∗ = 0 . (16)We can calculate the Lagrangian associated with thistrial wave function by integrating ¯ L (1 D ) over all 1D space:¯ L (1 D ) [ x , w , α , β ] = Z ∞−∞ d ¯ x ¯ L (1 D ) [ ψ ∗ ] . (17)This Lagrangian will depend only on the variational pa-rameters contained in the trial wave function. The 1Dequations of motion are found with the standard Euler–Lagrange equations. Next we will present the 1D and 3DLVM trial wave functions. C. N c
3D LVM trial wave function
In the 3D, N c –gaussian model we take the trial wavefunction to be a sum of N c three–dimensional Gaussianclouds. The j th cloud has its own initial momentum, ¯k j and set of variational parameters. These parameters arethe cartesian coordinates of the cloud center: ¯ x j , ¯ y j , and¯ z j ; the widths along the x , y , and z directions: ¯ w jx , ¯ w jy ,and ¯ w jz ; the linear phase coefficients along the x , y , and z directions: ¯ α jx , ¯ α jy , and ¯ α jz ; and the quadratic phasecoefficients along the x , y , and z directions: ¯ β jx , ¯ β jy ,and ¯ β jz . The j th cloud also has its own normalizationcoefficient, A j , which will be eliminated by fixing thenumber of atoms in each cloud.The mathematical form (in scaled units) of the trialwave function is the following:Ψ( ¯r , ¯ t ) = 1 √ N c N c X j =1 A j (¯ t ) e f j ( ¯r , ¯ t )+ i ¯k j · ¯r (18)where f j ( ¯r , ¯ t ) = X η = x,y,z − (¯ η − ¯ η j ) w jη + i ¯ α jη ¯ η + i ¯ β jη ¯ η ! . (19)This trial wave function will be used below in derivingthe 3D Lagrangian function from which the 3D equationsof motion will be obtained. D. N c
1D LVM trial wave function
In the 1D, N c –gaussian–cloud model we take the trialwave function to be a sum of N c one–dimensional Gaus-sian clouds. The j th cloud has its own initial momentum,¯ k j and set of variational parameters. These parametersare the cartesian coordinate of the cloud center, ¯ x j , thecloud width, ¯ w j , the linear phase coefficient, ¯ α j , and thequadratic phase coefficient, ¯ β j . The j th cloud also hasits own normalization coefficient, A j , which will be elim-inated by fixing the number of atoms in each cloud.The mathematical form (in scaled units) of the trialwave function is the following: ψ (¯ x, ¯ t ) = 1 √ N c N c X j =1 A j (¯ t ) e f j (¯ x, ¯ t )+ i ¯ k j ¯ x (20)where f j (¯ x, ¯ t ) = − (¯ x − ¯ x j (¯ t )) w j (¯ t ) + i ¯ α j (¯ t )¯ x + i ¯ β j (¯ t )¯ x (21) III. THE LVM EQUATIONS OF MOTION
In this section we derive the 3D and 1D LVM equationsof motion. We first describe our assumed constraintsplaced on the trial wave function, then we sketch thederivation of the 3D and 1D Lagrangians for the giventrial wave functions and finally sketch the derivation ofthe final LVM equations of motion. Most of the de-tails are relegated to appendices. The terms in theseLagrangians containing the external potential ¯ V ext onlydepend on the cloud centers, x , and widths, w , no matterthe form of the potential. This is due solely to our choiceof trial wave function. Thus the 1D and 3D Lagrangianscan be written in terms of a “variational” external po-tential, U ext ( x , w ). This enables the terms in the finalequations of motion that account for the external poten-tial to be written in terms of the x and w gradients ofthis potential. Hence the final equations of motion canbe written in a form independent of the particular V ext . A. Constraints on the wave function
Here we make several assumptions about the physicalsystem which have material effects on the values of thevariational parameters. These are as follows:1. We assume that each of the N c clouds are mov-ing at sufficiently different velocities such thatany integral of a quantity containing a factor like e i ( ¯ k jη − ¯ k j ′ η ) ¯ η where j = j ′ can be neglected. If theclouds move with sufficiently different velocities,these factors will be rapidly oscillating and theirintegrals can be neglected.2. The amplitudes, A j (¯ t ), are all real. Here we as-sume that the phase of each cloud’s wave functionis quadratic so that there is no phase is subsumedby the amplitudes. The quadratic phase enablesthe cloud centers to move and the cloud widths toexpand and contract. 3. The number of atoms in each cloud is fixed. Cloudsdo not exchange atoms. This plus the normaliza-tion condition, fixes a relationship (derived below)between A j and the widths ¯ w jη where η = x, y, z .We can use these assumptions plus the normalizationcondition on the trial wave function to derive conditionsthat constrain the values of the A j .To find these conditions we require that the full trialwave function be normalized to unity:1 = Z d ¯ r | Ψ( ¯r , ¯ t ) | = 1 N c N c X j =1 A j (¯ t ) (cid:16) π / ¯ w jx (¯ t ) ¯ w jy (¯ t ) ¯ w jz (¯ t ) (cid:17) , (22)where we have used assumptions 1 and 2 above to re-duce the normalization condition to a sum of Gaussianintegrals which are easily evaluated. This last expressionis the condition for the trial wave function to be nor-malized. However, our assumption that the number ofatoms in each cloud is fixed adds a further restriction tothe above expression. That is that each cloud is individ-ually normalized. This gives, finally, A j (¯ t ) π / ¯ w jx (¯ t ) ¯ w jy (¯ t ) ¯ w jz (¯ t ) = 1 , j = 1 , . . . , N c . (23)These constraints together automatically satisfy Eq. (22)and enable the elimination of all of the A j in the final 3DLagrangian. A similar condition on the norm of eachcloud in the 1D case can be derived and the result is A j (¯ t ) π / ¯ w j (¯ t ) = 1 , j = 1 , . . . , N c . (24)Using these conditions we can now derive the 3D and 1DLagrangians. B. The LVM model Lagrangians
The LVM Lagrangians for the particular choice of trialwave function given above can be derived by insertingEq. (18) into Eq. (13) and performing the required inte-grations. The details are given in Appendix A. The resultfor the 3D Lagrangian is¯ L (3 D ) = 1 N c N c X j =1 X η = x,y,z " ˙¯ α jη ′ ¯ η j + ˙¯ β jη (cid:18) ¯ η j + 12 ¯ w jη (cid:19) + 12 ¯ w jη + 2 ¯ β jη ¯ w jη + (cid:16) β jη ¯ η j + ¯ α jη + ¯ k jη (cid:17) + ¯Ω z N c N c X j =1 " ¯ y j (cid:0) ¯ α jx + ¯ k jx + 2 ¯ β jx ¯ x j (cid:1) − ¯ x j (cid:0) ¯ α jy + ¯ k jy + 2 ¯ β jy ¯ y j (cid:1) + 12 N c (cid:16) ¯ U (3 D )ext ( x , w )+ ¯ U (3 D )int ( x , w ) (cid:17) . (25)The variational potentials are defined as¯ U (3 D )ext ( x , w ) ≡ N c Z d ¯ r ¯ V ext ( ¯r , ¯ t ) | Ψ( ¯r , ¯ t ) | (26)¯ U (3 D )int ( x , w ) ≡ N c (cid:18)
12 ¯ gN (cid:19) Z d ¯ r | Ψ( ¯r , ¯ t ) | . (27)The derivation of the Lagrangian for the 1D case isvery similar. The major difference is the absence of arotating–frame term. Thus we will simply define the 1Dvariational potentials and present the 1D Lagrangian re-sult.The 1D external variational potential can be writtenas¯ U (1 D )ext ( x , w ) ≡ N c ¯ L (1 D )3 ( x , w )= 2 N c Z + ∞−∞ d ¯ x ¯ V ext (¯ x, ¯ t ) ψ ∗ (¯ x, ¯ t ) ψ (¯ x, ¯ t )= N c X j =1 (cid:18) π / ¯ w j (cid:19) Z + ∞−∞ d ¯ x × exp ( − (¯ x − ¯ x j ) ¯ w j ) ¯ V ext (¯ x, ¯ t ) (28)Where we have used constraints 1 and 3 to simplify theintegrals. In order to apply this 1D model to a particularsystem, this potential must be calculated. Derivatives of¯ U ext appear in the equations of motion.The 1D interaction variational potential is defined by¯ U (1 D )int ( x , w ) ≡ N c ¯ L (1 D )4 ( x , w )= (2 N c ) 12 ¯ g D N Z ∞−∞ d ¯ x | ψ | (29)We can use these definitions to write the final form ofthe 1D Lagrangian as follows.¯ L (1 D ) = 1 N c N c X j =1 ˙¯ α j ¯ x j + ˙¯ β j (cid:18) ¯ x j + 12 ¯ w j (cid:19) + 12 ¯ w j + 2 ¯ β j ¯ w j + (cid:16) β j ¯ x j + ¯ α j + ¯ k j (cid:17) ! + 12 N c (cid:16) ¯ U (1 D )ext ( x , w ) + ¯ U (1 D )int ( x , w ) (cid:17) (30)Explicit expressions for the 1D and 3D interaction poten-tials (which are always the same) are derived in AppendixB. With the 3D and 1D Lagrangians in hand we can nowderive the equations of motion for the variational param-eters. C. LVM model equations of motion
The 3D equations of motion for the Gaussian center co-ordinates, widths, and linear and quadratic phase param-eters can obtained in two steps: (1) derive the E–L equa-tion of motion for each variational parameter by using the 3D Lagrangian (Eq. (25)) in the standard E–L equa-tion of motion (Eq. (6)) to obtain a first–order differentialequation in time, (2) take the second time derivative ofthe equations containing ˙¯ η j and ˙¯ w jη where ( j = 1 , . . . , N c and η = x, y, z ) and use the other equations to eliminatethe other variables from these equations. The equationsfor the centers and widths together form a closed system.The LVM 3D equations of motion thus consist of apair of second–order ordinary differential equation for thecloud centers and widths as well as expressions for the¯ β jη and the ¯ α jη in terms of the centers, widths and theirfirst derivatives. The derivation of these equations canbe found in Appendix C.¨¯ x j = 2 ¯Ω z ˙¯ y j + ¯Ω z ¯ x j − ∂ ¯ U (3 D ) ∂ ¯ x j , (31a)¨¯ y j = − z ˙¯ x j + ¯Ω z ¯ y j − ∂ ¯ U (3 D ) ∂ ¯ y j , (31b)¨¯ z j = − ∂ ¯ U (3 D ) ∂ ¯ z j , (31c)¨¯ w jη = 4¯ w jη − ∂ ¯ U (3 D ) ∂ ¯ w jη , (31d)¯ β jη = ˙¯ w jη w jη , (31e)¯ α jx = ( ˙¯ x j − ¯Ω z ¯ y j ) − β jx ¯ x j − ¯ k jx , (31f)¯ α jy = ( ˙¯ y j + ¯Ω z ¯ x j ) − β jy ¯ y j − ¯ k jy , (31g)¯ α jz = ˙¯ z j − β jx ¯ x j − ¯ k jx , (31h) η = x, y, z j = 1 , . . . , N c In the above, ¯ U (3 D ) ≡ ¯ U (3 D )ext + ¯ U (3 D )int . The equations forthe cloud centers and cloud widths (Eqs. (31a), (31b),(31c), and (31d)) form a closed set that contain only the¯ η j , ˙¯ η j , ¯ w jη , and ˙¯ w jη . Once these quantities are obtained,all of the other variational parameters can be calculated.The equations of motion for the 1D case are derivedsimilarly. We will simply state the result here. They con-sist of a pair of second–order ordinary differential equa-tions for the cloud centers and widths as well as expres-sions for the ¯ β j and the ¯ α j in terms of the centers, widthsand their time derivatives:¨¯ x j = − ∂ ¯ U (1 D ) ∂ ¯ x j , (32a)¨¯ w j = 4¯ w j − ∂ ¯ U (1 D ) ∂ ¯ w j , (32b)¯ β j = ˙¯ w j w j , (32c)¯ α j = ˙¯ x j − β j ¯ x j − ¯ k j , (32d) j = 1 , . . . , N c Similar to the 3D case, ¯ U (1 D ) ≡ ¯ U (1 D )ext + ¯ U (1 D )int . Theequations for the cloud centers and cloud widths (Eqs.(32a) and (32b)) form a closed set that contain only the (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (a) (b) (c)(f) (g)(d) (e) SMSMSM SMSM SMSM(source mass) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)
FIG. 1: (color online) 1D Atom interferometry scheme for measuring G . (a) A BEC is formed in the presence of an harmonictrap plus source mass (SM). (b)–(c) Initial Split Phase 1: the condensate is split by laser light and the two condensate piecesallowed to fly apart with initial speeds ± v until they stop at which point the harmonic trap is turned off. (d) Initial SplitPhase 2: the two condensate pieces now experience differential gravitational accelerations due to the SM during a wait time, T W . (e) Initial Split Phase 3: the harmonic trap is turned on pushing the condensate pieces back together. (f)–(g) Final Split:once they overlap they are split again creating four clouds: two that fly apart with approximate speeds ± v and two that aremotionless except for the small relative velocity developed during the wait time. ¯ x j , ˙¯ x j , ¯ w j , and ˙¯ w j . Once these quantities are obtained,all of the other variational parameters can be calculated.These equations of motion for the variational param-eters of the 1D and 3D trial wave functions form thecentral result of this paper. We note one more time thatthese equations hold for any external potential, ¯ V ext . The1D and 3D N c –Gaussian–cloud wave function ansatz hasseveral areas of application including atom interferome-tery involving BECs and BECs moving in waveguides. IV. AI MEASUREMENT OF G INMICROGRAVITY
As an example of the use of the LVM model, we de-scribe how it could be applied to the design of a measure-ment of Newton’s universal gravitational constant, G , ina microgravity environment. Thus we consider the ide-alized 1D atom–interferometry sequence shown in Fig. 1where two pieces of a BEC are separated, differentiallypulled upon by a source mass (SM), recombined, andsplit again. In this section we describe this AI sequencein detail.To extract the value of G requires an interference pat-tern where only the effects of the SM are present. Thiscan be obtained by repeating the above AI sequencemany times with the SM present so that averaging overthe interference patterns produced would wash out ran-dom effects of the environment. Systematic effects of theenvironment as well as the effect of the SM would still bepresent in the averaged pattern. To remove systematiceffects of the environment the SM would then be takenaway and the AI sequence repeated many more times andthese SM–absent interference patterns would again be av-eraged. The difference between the averaged SM–presentand SM-absent patterns would leave a pattern in which only the effects of the SM are present. The pattern dueto the SM where environmental effects are neglected canbe approximately simulated by the variational methoddescribed here.The basic idea of the AI sequence is that the relativevelocity developed between the two condensate pieces dueto the differential acceleration produced by the gravityof the source mass. This relative velocity is imprintedon the final–state interference pattern since the velocitydistribution is proportional to the condensate phase gra-dient.The value of G can then be extracted from the inter-ference data. This sequence is imagined as the essentialstep in an AI measurement of G . It is assumed to be car-ried out in a microgravity environment such as NASA’sCold Atom Laboratory aboard the International SpaceStation. A more sophisticated version of this sequencewould be required for a precision measurement of G un-der such conditions.The basic sequence begins with a condensate formedat the center of a 1D harmonic trap with frequency ω T with a SM present (Fig. 1(a)). The overall sequence con-sists of an “initial split” where this initial condensate isseparated into two pieces which are then rejoined and a“final split” where the rejoined pieces are split again andthe interference pattern is imaged. The initial split (IS)has three phases as we now describe.The speed v imparted to the condensate pieces is as-sumed in this work to come from applying identical coun-terpropagating laser beams to the condensate. An atomwill receive a momentum kick by absorbing a photonfrom one beam and emitting it into the other beam.The momentum imparted to a condensate atom is 2¯ hk L where k L = 2 π/λ L where λ L is the wavelength of thelaser. Laser pulses can be engineered so that this basicmomentum–kick process happens m times so that the ve- BE C den s i t y ( a r b . un i t s ) x ( µ m)BEC halves fly apart, trap ON LVMGPE x ( µ m)BEC halves are stationary, trap OFF LVMGPE x ( µ m)BEC halves fly back together, trap ON LVMGPE
FIG. 2: (color online) Comparison of the evolution of the condensate density during the “initial split” phase of the AI scheme(see Fig. 1 panels (a)–(f)). The left panel shows the LVM (red, solid lines) and GPE (blue, dashed lines) results at four differenttimes as the condensate pieces are flying apart and the trap is on. Note that each curve consists of two symmetically placedpeaks and only the left peak is given a time label. The middle panel shows their evolution at the beginning and end of theperiod when the trap is off. The right panel shows the evolution after the trap is turned back on and the condensate piecesre–overlap. The conditions for this case are N = 10 Rb atoms, the trap potential frequency is ω T / π = 1 Hz, the wait timeis T W = 500 ms., the initial velocity is v = 3 . × − m/s and there is no source mass present. locity kick the atoms receive is v = 4 π ¯ hm/ ( M λ L ) where M is the mass of the condensate atom. In the examplesgiven below we assume that m = 8 so that, for Rbatoms, v = 9 . × m/s.In IS phase 1 the condensate is separated into two equalpieces by a sequence of optical–lattice pulses [26] the totaleffect of which is to change the velocity of one conden-sate piece by + v and the other by − v (Fig. 1(b)). Theinitially motionless condensate pieces thus fly apart withvelocities − v (blue circle to the left) and + v (red cir-cle to the right) and both come to a stop at the turningpoints of the harmonic potential after one quarter of thetrap period. We call this time t = T / and is shown inFig. 1(c).In IS phase 2 the harmonic trap is turned off at time t = T / ideally leaving the two pieces motionless. Thesystem is then allowed to evolve during the time interval T / ≤ t ≤ T / + T w . We call this evolution time, T w ,the “wait period”. This is shown in Fig. 1(d). Duringthis time the gravitational pull of the SM causes a differ-ential acceleration of the two condensate pieces since onepiece is further away from the SM than the other. Thisdifferential acceleration causes the two pieces to developa relative velocity.In IS phase 3 the trap is then turned back on whichbrings the condensate pieces back together during thetime interval T / + T w ≤ t ≤ T / + T w as shown inFig. 1(e). Once they overlap again the blue (left) piecewill be moving at velocity + v plus a small difference dueto the SM and the red (right) piece will be moving at − v plus a small difference. The initial split takes placeduring the time interval 0 ≤ t ≤ T / + T w ≡ T is andshown in Figs. 1(a)–(e).Once the two pieces are overlapped again the same setof optical–lattice pulses is applied as before (Fig. 1(f)).This causes the red condensate piece to split in two witha piece moving at velocity approximately − v and theother at zero velocity plus a small amount due to the pull of the SM. The blue half is also split into a piece movingat approximate velocity +2 v and the other approximatelymotionless except again for a small deviation. These twonearly motionless pieces will remain overlapped and willhave the small relative velocity that was developed duringthe wait time when the trap was turned off. We assumethat the wait time is much longer than a quarter periodof the harmonic trap ( T w ≫ T / ).Finally, the two fast condensate pieces will fly awayleaving the nearly motionless pieces behind as shownin Fig. 1(g). We will refer to this sequence, shown inFigs. 1(f)–(g), as the “final split”. Imaging this middlecloud will leave an interference pattern due to the rela-tive velocity of the two condensate pieces and this patternwill be due only to the differential gravitational pull ofthe SM during the wait time T W . V. APPLYING THE LVM MODEL TO THE AISEQUENCEA. Comparison of the LVM model and GPEsolutions
In this section we show how the LVM model can beused to analyze the atom interferometry sequence de-scribed above. We emphasize here that the LVM modelis not sufficiently accurate to be used in the analysis ofa precision measurement rather it is to be used in thepreliminary design of the experiment. To show this, wesimulated the evolution of the condensate wave functionduring the initial split of this AI sequence using the 1Dtime–dependent GPE and the 1D LVM model describedearlier. The comparison is shown in Fig. 2. The detailsof how the simulation was conducted can be found inAppendix D.This comparison shows that, while the Gaussian–approximate trial wave functions only qualitatively agreetheir GPE counterparts, the motion of the GPE and LVMwave packet centers and widths track very well. The vari-ational wave function will only agree with the exact solu-tion to the extent that variational wave function can bemodified to fit the exact solution by varying its parame-ters.Modeling performed for the final design of a precisionAI measurement will require the numerical solution ofthe full 3D Gross–Pitaevskii equation. Since experimen-tal conditions for precision measurements are typicallyextreme, numerical solution of the 3D GPE will requirea substantial effort. This effort usually makes solving the3D GPE too expensive for the purposes of preliminary AIdesign. We shall elaborate on this point below.The purpose of the LVM model is to facilitate the pre-liminary design of a precision measurement. It can servethis purpose in several ways. First, it can provide a sim-ple and physically intuitive picture of how the measuredquantity can be extracted from the experimental data.In some cases an approximate analytical expression forthe measured quantity in terms of experimental param-eters can be derived. Second, the model can be usedto estimate the sensitivity of the measurement. Finally,this expression can be used to speed up significantly thetime spent on the preliminary design of an AI sequenceby providing estimates of the values of these parametersthus constraining the size of the experimental parameterspace. We illustrate these features of the LVM model byapplying it to the AI sequence described above.
B. Approximating the interference pattern usingthe LVM model
The equations of motion for the case when an harmonictrap with frequency ω T is present along with a sourcemass of mass M SM and located at x SM can be found bycomputing the derivatives found in Eqs. (32a) and (32b).The external varational potential for this case, given byEq. (E8), is derived in Appendix E and the interactionvariational potential is given in Eq. (B3) of Appendix B.The resulting equations of motion have a striking sim-ilarity to equations of motion produced by Newton’s sec-ond law. They can be written as follows:¨¯ x j + (cid:0) ¯ ω T − ω SM (cid:1) ¯ x j = ¯ ω SM ¯ x SM − X j ( x , w ) , (33a)¨¯ w j + (cid:0) ¯ ω T − ω SM (cid:1) ¯ w j = 4¯ w j + ¯ gN (2 π ) / ¯ w j − W j ( x , w ) j = 1 , . (33b)The terms X j ( x , w ) and W j ( x , w ) represent “repulsionforces” due to the overlap of different clouds and are given by X j ( x , w ) = (cid:18) gNπ / (cid:19) (34) × X j = j (¯ x j − ¯ x j ) exp (cid:26) − ( ¯ x j − ¯ x j ) ¯ w j + ¯ w j (cid:27)(cid:0) ¯ w j + ¯ w j (cid:1) / and W j ( x , w ) = (cid:18) gNπ / (cid:19) X j = j ¯ w j exp (cid:26) − ( ¯ x j − ¯ x j ) ¯ w j + ¯ w j (cid:27)(cid:0) ¯ w j + ¯ w j (cid:1) / × h x j − ¯ x j ) − (cid:0) ¯ w j + ¯ w j (cid:1)i , (35)where j = 1 , N c = 2 from Eq. (20) to de-rive the condensate density at the end of the sequence,and (2) use approximate analytical solutions to the equa-tions of motion given by Eqs. (33a) and (33b) to obtainexpressions for the variational parameters ¯ x , ˙¯ x , ¯ x , ˙¯ x ,¯ w , ˙¯ w , ¯ w , and ˙¯ w at the end of the sequence.The condensate trial wave function at the end of theinitial split ( t = T is ) for two clouds ( N c = 2) is given byEq. (20) and can be written as follows: ψ is (¯ x, ¯ T is ) = (2 π / ¯ w ) − / e − (¯ x − ¯ x w + i (( 12 ˙¯ x − β ¯ x )¯ x + ¯ β ¯ x ) + (2 π / ¯ w ) − / e − (¯ x − ¯ x w + i (( 12 ˙¯ x − β ¯ x )¯ x + ¯ β ¯ x ) ≡ √ (cid:0) ψ is , (¯ x, ¯ T is ) + ψ is , (¯ x, ¯ T is ) (cid:1) (36)where we have used Eq. (32d). All variational parame-ters above are evaluated at t = T is . Since ¯ β j = ˙¯ w j / (4 ¯ w j )(from Eq. (32c) the wave function is expressed entirely interms of the ¯ x j , ¯ w j and their derivatives. These quanti-ties can be found by solving Eqs. (33a) and (33b).We can model the effect of the second splitting by mul-tiplying ψ (¯ x, ¯ T is ) by a factor which splits each exisitingcloud into two clouds: one which is boosted by v andone boosted by − v . Thus, in scaled units, the final–splitwave function at t = T is becomes ψ fs (¯ x, ¯ T is ) = 12 (cid:16) e + i ¯ v ¯ x/ + e − i ¯ v ¯ x/ (cid:17) ( ψ is , + ψ is , )= 12 e − i ¯ v ¯ x/ ψ is , + 12 e + i ¯ v ¯ x/ ψ is , + 12 e − i ¯ v ¯ x/ ψ is , + 12 e + i ¯ v ¯ x/ ψ is , ≡ ψ , + ψ , + ψ , + ψ , . (37)The final–split wave function consists of four clouds. Thevelocities of these clouds just after the second splittingare ≈ ≈ +2 v forcloud (1,2) and ≈ − v for cloud (2,1).Clouds (1,2) and (2,1) fly rapidly away from the centerand after a short time, τ , the overlap of the fast cloudswith those at the center can be neglected. At this timein the AI sequence the density of the condensate can beimaged to obtain an interference pattern. To obtain asimple approximation for the condensate density at thistime we assume that the change in the variational param-eters during the time τ can be neglected and we approx-imate their values at time t = T is + τ with their valuesat t = T is .The condensate density at the trap center can be foundfrom the squared modulus of the sum of the wave func-tions of the two nearly motionless clouds. The result forthe condensate density near the trap center after the finalsplit is (in SI units) ρ ( x, T is ) = | ψ , + ψ , | = ( π w ) − e − ( x − x ) /w (38)+ ( π w ) − e − ( x − x ) /w + ( π w ) − ( π w ) − × e − ( x − x ) / (2 w ) − ( x − x ) / (2 w ) cos ( φ ( x )) , where φ ( x ) = k f x + ( β − β ) x (39)is the final condensate phase difference and k f = M ¯ h ( ˙ x − ˙ x ) + v + 2 β x − β x . (40)The quantity k f is the spatial frequency near the center ofthe interference pattern. We note that these variationalquantities are understood to be evaluated at t = T is . Wecan now find an approximate interference pattern by solv-ing Eqs. (33a) and (33b) for the variational parameterseither numerically or by finding approximate analyticalsolutions.The presence of a source mass has two major effects onthe interference pattern: (1) it creates interference fringesand (2) the peak of the pattern is shifted away from thezero–mass position towards the source–mass. We willcall the width of the central fringe, λ f , and the shift ofthe peak due to the source mass is called, x c . Later wewill show how G can be obtained from the experimentalparameters and either of these two characteristics of theinterference pattern.The final condensate phase, φ ( x ), predicts an inter-ference pattern that displays equally spaced fringes nearthe center of the pattern due to the linear term. Thefringe spacing narrows away from the center due to thequadratic term. The gradient of the phase measures thelocal relative velocity of the overlapping clouds.This relative velocity is caused by (1) the differentialaccelerations of the two clouds during the wait period and t=T t=T is c l oud w i d t h ( mm ) time (s) FIG. 3: (color online) Cloud width versus time, w ( t ), dur-ing the IS from the numerical solution of the full system ofequations, (33a) and (33b). The conditions are M SM = 100kg, x SM = 0 . v = 9 . × − m/s, ω T / π = 1 Hz,and T W = 10 seconds. The vertical dotted lines indicate theend–times of the three phases of the IS. (2) the changing widths of the two clouds due to conden-sate atom–atom interactions. Interactions also cause thefringe spacing to narrow away from the center of the in-terference pattern. The β parameters in the LVM modeldescribe how the clouds expand and/or contract due tothe competition between the repulsive interactions andthe harmonic confinement. They also account for anyshift, due to atom–atom interactions, in the value of G inferred from the interference pattern. C. Extracting G from the inteference pattern In modeling a precision measurement the value of G would be inferred from repeated GPE simulations of theexperiment using different values of G until a satisfactorymatch between the measured and simulated interferencepatterns was obtained. This is a necessary (but very ex-pensive, see below) procedure. Such repeated simulationswith different G values can also be performed using theLVM model even though its accuracy is not sufficient tomake a final decision about any particular AI scheme.The LVM model can also help in designing an exper-iment by providing an approximate expression for G interms of the experimental parameters and characteristicsof the interference pattern. This can be done by approx-imating the solution of the LVM equations of motion foreach step of the proposed AI sequence. We now showthis by deriving expressions for G for the AI sequencedescribed in Section IV.To find approximate solutions of Eqs. (33a) and (33b)we will make three assumptions: (1) the gravity effecton the condensate is small compared to the harmonic0 (a) BE C den s i t y ( a r b . un i t s ) x ( µ m)T W = 5 s, G = 6.13x10 -11 m kg -1 s -2 λ f = 163.0 µ m 0 0.02 0.04 0.06 0.08 0.1-300 -200 -100 0 100 200 300 (b) x ( µ m)T W = 10 s, G = 6.47x10 -11 m kg -1 s -2 λ f = 38.6 µ m 0 0.02 0.04 0.06 0.08 0.1-300 -200 -100 0 100 200 300 (c) x ( µ m)T W = 20 s, G = 6.59x10 -11 m kg -1 s -2 λ f = 9.48 µ m FIG. 4: (color online) Interference patterns for three different wait times showing how the width of the central fringe decreasesas T W increases. (a) T W = 5 s, λ f = 163 . µ m (b) T W = 10 s, λ f = 38 . µ m and (c) T W = 20 s, λ f = 9 . µ m. Except forthe wait times the conditions are the same as in Fig. 3: M SM = 100 kg, x SM = 0 . v = 9 . × − m/s, ω T / π = 1 Hz.The resulting values for G are (a) G = 6 . × − m kg − s − , (b) G = 6 . × − m kg − s − , and (c) G = 6 . × − m kg − s − . confinement so that ω SM ≪ ω T (see Eq. (42)), (2) theoverlap terms, X j and W j , can be neglected, and (3)the condensate width at the end of IS phase 2 can beapproximated by the product of the IS phase 2 widthexpansion rate and the wait time, i.e., ¯ w ( ¯ T ) ≈ ˙¯ w ( ¯ T ) ¯ T W .One immediate consequence of this is that we can ne-glect ω SM in the equations of motion for the ¯ w j andso they are approximately unaffected by gravity dur-ing the entire AI sequence. With these approxima-tions we can also infer that the time evolution of thewidths of the two clouds are approximately the same,i.e., ¯ w (¯ t ) ≈ ¯ w (¯ t ) and also ˙¯ w (¯ t ) ≈ ˙¯ w (¯ t ). This alsoimplies that ¯ β (¯ t ) ≈ ¯ β (¯ t ).These approximations also enable us to find analyticalsolutions of Eq. (33a) for the cloud center positions ineach of the three phases of the initial split (panels (c),(d), and (e) of Fig. 1, respectively). In phases 1 and 3the solutions are sin(¯ ω T ¯ t ) and cos(¯ ω T ¯ t ) and, in phase 2,the solutions are sinh(¯ ω SM ¯ t ) and cosh(¯ ω SM ¯ t ). Using theinitial conditions ¯ x (0) = ¯ x (0) = 0, ˙¯ x (0) = − ¯ v , and˙¯ x (0) = +¯ v we can solve equations of motion for thecenter coordinates in each IS phase using the final valuesof ¯ x j and ˙¯ x j in the previous phase as the initial values inthe next phase.This straightforward procedure enables us to obtainapproximate expressions for the center positions and ve-locities at t = T is . The results are¯ x j ( ¯ T is ) = √ ω SM ¯ ω T (cid:18) ¯ x SM + ǫ j ¯ v ¯ ω T (cid:19) sinh( φ SM ) (41)˙¯ x j ( ¯ T is ) = ¯ ω T (cid:20) ¯ x SM − (cid:18) ¯ x SM + ǫ j ¯ v ¯ ω T (cid:19) cosh( φ SM ) (cid:21) , where j = 1 , ǫ = − ǫ = +1, φ SM ≡ √ ω SM T W and from Eq. (E6) ω SM = GM SM | x SM | ! / . (42) It is also possible to find an approximate value for¯ β ( ¯ T is ) (assumed to be equal to ¯ β ( ¯ T is )) by approxi-mating the width equation of motion (Eq. (33b)) duringphase 3 of the IS. An example of how the width behavesduring this time is shown in Fig. 3 between the two right-most vertical dotted lines. The cloud width is initiallyquite large and drops rapidly to a value close to its start-ing value.Thus during the IS phase 3 period all of the interactionterms in Eq. (33b) can be neglected and thus the equationof motion can be approximated as (for cloud 1)¨¯ w + ¯ ω T ¯ w ≈ , (43)with initial conditions ¯ w ( ¯ T ) and ˙¯ w ( ¯ T ) at ¯ t = ¯ T . Thiscan be easily solved and, since the duration of phase 3 isone–quarter of the trap period, we can write ¯ w ( ¯ T is ) and˙¯ w ( ¯ T is ) in terms of ¯ w ( ¯ T ) and ˙¯ w ( ¯ T ) as follows¯ w ( ¯ T is ) = ˙¯ w ( ¯ T )¯ ω T (44)and ˙¯ w ( ¯ T is ) = − ¯ w T ¯ w ( ¯ T ) ≈ − ¯ w T ˙¯ w ( ¯ T ) ¯ T W , (45)where in the last equality we have used the approxima-tion mentioned above that the cloud width at the end ofIS phase 2 is approximately the width velocity at t = T multiplied by the wait time, T W . So, finally we have anapproximate expression for the β parameters at t = T is :¯ β ( ¯ T is ) ≈ ¯ β ( ¯ T is ) = ˙¯ w ( ¯ T is )4 ¯ w ( ¯ T is ) ≈ −
14 ¯ ω T ¯ T W . (46)This expression enables us to derive an expression for G solely in terms of the experimental parameters andcharacteristics of the measured interference pattern.We can use the expressions for the ¯ x j ( ¯ T is ) and the˙¯ x j ( ¯ T is ) from Eqs. (42)) and the ¯ β j ( ¯ T is ) in Eq. (46) to ap-proximate k f in Eq. (40). This equation can be written1as (in SI units)¯ hk f M v = π ¯ hM vλ f = 1+ φ SM sinh( φ SM ) − cosh( φ SM ) . (47)The left–hand–side (LHS) of this equation contains onlyexperimental parameters ( M and v ) and a characteristicof the interference pattern ( λ f is the width of the centralfringe). The value of the LHS can thus be used to findthe value of φ SM = √ ω SM T W from which G can beobtained. Thus this expression can be used to estimatethe value of G from the experimental parameters and theinterference pattern.This is most easily seen if we assume that φ SM ≪ G can be furtherapproximated as k f ≈ M vφ SM / h . So, in SI units, wecan write G ≈ (cid:18) ¯ hx SM M vM SM T W (cid:19) k f = (cid:18) ¯ hx SM M vM SM T W (cid:19) (cid:18) πλ f (cid:19) (48)where M is the mass of a condensate atom and λ f is thewidth of the central fringe of the interference pattern.This formula expresses G in terms of the experimentalparameters, T W , x SM , M SM , and v , and, λ f , the widthof the central fringe of the interference pattern. Thisquantity is easily inferred from the pattern. Interferencepatterns for several sets of experimental parameters areshown in Fig. 4. In particular, the figure shows how thepatterns change as the wait time, T W , is increased. Thefigure also shows the values of G given by the above for-mula for each pattern.It is also possible to write an expression for G in termsof the shift, x c , of the interference pattern maximumaway from the maximum of the pattern obtained whenthe source mass is absent. The result is x c ≈ ( x + x ) x c = x SM ω T T W φ SM sinh( φ SM ) x c ≈ x SM ω T T W φ SM , φ SM ≪ G ≈ (cid:18) x SM ω T T W M SM (cid:19) x c . (49)Equations (48) and (49) display the ability of LVM modelto show how G can be obtained from the experimentalparameters and the data. They can further be used toestimate the sensitivity of the measured G value to eachof the experimental parameters and to develop a prelim-inary error budget for the experiment. D. Using the LVM model for preliminary AI design
One of the features of the LVM model is that differentAI schemes can be rapidly evaluated and compared. The scheme presented above only splits the condensate intotwo fast pieces. Other possible AI schemes might, for ex-ample, split the condensate into three, five or more pieces(see Ref. [26]). Whether such ideas are useful can only bedetermined by modeling the particular AI sequence en-visioned. In general this will require many iterations inorder to come to a conclusion about the candidate AI se-quence. We have also shown above that the model is alsouseful in shrinking the experimental parameter space.If this modeling is conducted by numerical solution ofthe GPE, this will require many runs each of which isvery expensive. For example, in the simulation shown inFig. 2 the GPE run needed 10 seconds to finish using afast desktop computer while the LVM simulation required < T W = 0 . v = 3 . × − m/s was very different from themuch more extreme conditions found in Fig. 4(c) where T W = 20 s and v = 9 . × − m/s. The size of the gridbox in the second simulation was 30 times larger and theevolution time was 40 times longer making the full sim-ulation take 1200 times longer. Thus, if the same GPEsimulator used for the first simulation were used to per-form the second simulation, it would have taken the desk-top computer more than 10 seconds (or about 4 months)to finish this single simulation. The LVM simulation tookless than three minutes. If dozens or hundreds of prelim-inary modeling runs were needed, it would clearly notbe practical to use the GPE even if such modeling weredone on a supercomputer. Furthermore, a realistic AIexperiment for measuring G would take place in threedimensions and would require 3D modeling. This wouldmake simulation time for the GPE even longer but wouldnot 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 precisionmodeling can carried out using the GPE only for a smallnumber of likely AI candidate schemes. A number ofsuch likely 3D schemes will be the subject of a futurepublication. VI. SUMMARY
In this paper we presented a variational techniquefor approximating the solution of the time–dependentGross–Pitaevskii equation, for situations where the con-densate is split into multiple parts as is common in atominterferometry processes. The technique presented canapproximate the solution of either the one–dimensionalGPE or the three–dimensional GPE in the laboratory orrotating frame.Each part of the condensate is modeled as a Gaus-sian with time–dependent center, width, and linear andquadratic phase parameters that enable the centers andwidths to change. Different parts of the condensate areallowed to overlap but exchange of atoms between dif-ferent condensate parts is assumed to be negligible. The2resulting equations of motion for the center coordinates,widths, and phase parameters is a system of ordinarydifferential equations in time that can be rapidly solved.The power of this model is that this set of equationsis cast in terms of derivatives of an arbitrary externalpotential and thus can be applied to many different AIschemes. This model should enable rapid prototypingand design of novel AI schemes for applications especiallyin microgravity environments.
Acknowledgments
This material is based upon work supported by theU.S. National Science Foundation under grant numbersPHY–1004975, PHY–0758111, the Physics Frontier Cen-ter grant PHY–0822671 and by the National Institute ofStandards and Technology.3 [1] G. Tino and M. Kasevich,
Atom Interferometry , vol.188 (IOS Press, Amsterdam, 2014), 1st ed., ISBN9781614994473.[2] A. D. Cronin, J. Schmiedmayer, and D. E. Pritchard,Rev. Mod. Phys. , 1051 (2009).[3] A. Peters, K. Y. Chung, and S. Chu, Nature , 849(1999).[4] J. M. McGuirk, G. T. Foster, J. B. Fixler, M. J. Snadden,and M. A. Kasevich, Phys. Rev. A , 033608 (2002).[5] F. Sorrentino, Y.-H. Lien, G. Rosi, L. Cacciapuoti,M. Prevedelli, and G. M. Tino, New Journal of Physics , 095009 (2010).[6] T. Gustavson, P. Bouyer, and M. Kasevich, ClassicalQuantum Gravity , 2385 (2000).[7] B. Canuel, F. Leduc, D. Holleville, A. Gauguet, J. Fils,A. Virdis, A. Clairon, N. Dimarcq, C. J. Bord´e, A. Lan-dragin, et al., Phys. Rev. Lett. , 010402 (2006).[8] S. Fray, C. A. Diez, T. W. H¨ansch, and M. Weitz, Phys.Rev. Lett. , 240404 (2004).[9] S. Dimopoulos, P. W. Graham, J. M. Hogan, and M. A.Kasevich, Phys. Rev. Lett. , 111102 (2007).[10] M. G. Tarallo, T. Mazzoni, N. Poli, D. V. Sutyrin,X. Zhang, and G. M. Tino, Phys. Rev. Lett. , 023005(2014).[11] C. Burrage, E. J. Copeland, and E. Hinds, Journal ofCosmology and Astroparticle Physics , 042 (2015).[12] S. Dimopoulos, P. W. Graham, J. M. Hogan, M. A. Kase-vich, and S. Rajendran, Phys. Rev. D , 122002 (2008).[13] P. W. Graham, J. M. Hogan, M. A. Kasevich, and S. Ra-jendran, Phys. Rev. Lett. , 171102 (2013).[14] G. Rosi, F. Sorrentino, L. Cacciapuoti, M. Prevedelli,and G. M. Tino, Nature , 518 (2014).[15] G. Lamporesi, A. Bertoldi, L. Cacciapuoti, M. Prevedelli,and G. M. Tino, Phys. Rev. Lett. , 050801 (2008).[16] T. van Zoest, N. Gaaloul, Y. Singh, H. Ahlers, W. Herr,S. T. Seidel, W. Ertmer, E. Rasel, M. Eckart, E. Kajari,et al., Science , 1540 (2010).[17] S. T. Seidel, D. Becker, M. D. Lachmann, W. Herr, E. M.Rasel, and Quantus Collaboration, in APS Division ofAtomic, Molecular and Optical Physics Meeting Abstracts (2016).[18] D. Aveline, E. Elliot, J. Williams, and R. Thompson, in
APS Division of Atomic, Molecular and Optical PhysicsMeeting Abstracts (2017).[19] D. A. Pushin, M. Arif, and D. G. Cory, Phys. Rev. A ,053635 (2009).[20] V. M. Perez-Garcia, H. Michinel, J. I. Cirac, M. Lewen-stein, and P. Zoller, Phys. Rev. A , 1424 (1997).[21] M. Edwards, M. Krygier, H. Seddiqi, B. Benton, andC. W. Clark, Phys. Rev. E , 056710 (2012), URL https://link.aps.org/doi/10.1103/PhysRevE.86.056710 .[22] E. Infeld, M. Matuszewski, and M. Trippen-bach, Journal of Physics B: Atomic, Molecu-lar and Optical Physics , L113 (2006), URL http://stacks.iop.org/0953-4075/39/i=6/a=L02 .[23] L. Salasnich, A. Parola, and L. Reatto,Phys. Rev. A , 043614 (2002), URL https://link.aps.org/doi/10.1103/PhysRevA.65.043614 .[24] L. Salasnich and B. A. Malomed,Phys. Rev. A , 053620 (2009), URL https://link.aps.org/doi/10.1103/PhysRevA.79.053620 . [25] L. Pitaevskii and S. Stringari, Bose–Einstein Condensa-tion (Oxford University Press, 2003).[26] M. Edwards, B. Benton, J. Heward, and C. W. Clark,Phys. Rev. A , 063613 (2010).[27] P. Muruganandam and S. Adhikari, Computer PhysicsCommunications , 1888 (2009), ISSN 0010-4655. Appendix A: 3D LVM Lagrangian
Here we sketch the derivation of the 3D Lagrangian.This Lagrangian has five terms, defined in Eq. (13), andwe will derive ¯ L (3 D )1 , ¯ L (3 D )2 and ¯ L (3 D )5 explicitly. Fur-ther we will show that ¯ L (3 D )3 and ¯ L (3 D )4 only dependon the x and w variational parameters. This in turnmotivates the introduction of the “variational” poten-tials, ¯ U ext ( x , w ) and ¯ U int ( x , w ), which later appear inthe equations of motion for the Gaussian center coordi-nates and widths. The explicit expression for ¯ U int ( x , w )is given in Appendix B. The form of ¯ U ext ( x , w ) for thecase of a 1D harmonic trap plus the gravitational poten-tial of an arbitrarily placed point mass will be given inAppendix E.
1. Derivation of ¯ L (3 D )1 The ¯ L term of the Lagrangian has the form¯ L (3 D )1 = Z d ¯ r Im { Ψ ∗ ( ¯r , ¯ t )Ψ ¯ t ( ¯r , ¯ t ) } . (A1)The trial wave function and its time derivative areΨ( ¯r , ¯ t ) = 1 √ N c N c X j =1 A j (¯ t ) e f j ( ¯r , ¯ t )+ ¯k j · ¯r , Ψ ¯ t ( ¯r , ¯ t ) = 1 √ N c N c X j =1 (cid:16) ˙ A j a (¯ t ) + A j (¯ t ) ˙ f j ( ¯r , ¯ t ) (cid:17) × e f j ( ¯r , ¯ t )+ ¯k j · ¯r (A2)where the dot denotes partial differentiation with re-spect to ¯ t . Inserting these into Eq. (A1), neglecting therapidly oscillating terms and then taking the imaginarypart yields the following result¯ L (3 D )1 = 1 N c N c X j =1 A j (¯ t ) X η = x,y,z Z d ¯ r (cid:16) ˙¯ α jη ¯ η + ˙¯ β jη ¯ η (cid:17) × Y η ′ = x,y,z exp ( − (cid:0) ¯ η ′ − ¯ η ′ j (cid:1) ¯ w jη ′ ) = 1 N c N c X j =1 X η = x,y,z (cid:18) π / ¯ w jη (cid:19) (A3) × Z d ¯ r (cid:16) ˙¯ α jη ¯ η + ˙¯ β jη ¯ η (cid:17) exp ( − (¯ η − ¯ η j ) ¯ w jη ) , where in the second line we have used Eqs. (23) to elimi-nate the A j (¯ t ) factors. The Gaussian integrals appearingin this last expression can easily be evaluated. The finalresult for ¯ L (3 D )1 is ¯ L :¯ L (3 D )1 = 1 N c N c X j =1 X η = x,y,z (cid:18) ˙¯ α jη ¯ η j + ˙¯ β jη (cid:18) ¯ η j + 12 ¯ w jη (cid:19)(cid:19) (A4) Next we derive ¯ L (3 D )2 .
2. Derivation of ¯ L (3 D )2 The expression for ¯ L (3 D )2 is given by¯ L (3 D )2 ≡ X η = x,y,z Z d ¯ r Ψ ∗ ¯ η Ψ ¯ η . (A5)If we differentiate Ψ with respect to ¯ η ( η is x , y , or z ), multiply the result with its complex conjugate, andneglect the rapidly oscillating terms (which integrate tozero), we obtainΨ ∗ ¯ η Ψ ¯ η ≈ N c N c X j =1 A j (¯ t ) exp − X η ′ = x,y,z (cid:0) η ′ − η ′ j (cid:1) ¯ w jη ′ × (cid:0) ¯ α jη + ¯ k jη + 2 ¯ β jη η (cid:1) + ( η − η j ) ¯ w jη ! . (A6)Inserting this into Eq. (A5) gives¯ L (3 D )2 = 1 N c N c X j =1 η = x,y,z A j (¯ t ) Z d ¯ r (cid:0) ¯ α jη + ¯ k jη + 2 ¯ β jη η (cid:1) + ( η − η j ) ¯ w jη ! exp − X η ′ = x,y,z (cid:0) η ′ − η ′ j (cid:1) ¯ w jη ′ (A7)Again we are left with Gaussian integrals that can bestraightforwardly evaluated. The final result is¯ L (3 D )2 = 1 N c N c X j =1 X η = x,y,z
12 ¯ w jη + 2 ¯ β jη ¯ w jη + (cid:16) ¯ α jη + ¯ k jη + 2 ¯ β jη ¯ η j (cid:17) ! . (A8)In the above we have, again, used Eqs. (23) to eliminatethe A j (¯ t ) factors.
3. 3D Variational Potentials
Next we consider the terms L (3 D )3 and L (3 D )4 which are¯ L (3 D )3 ≡ Z d ¯ r ¯ V ext ( ¯r , ¯ t ) | Ψ( ¯r , ¯ t ) | (A9)¯ L (3 D )4 ≡
12 ¯ gN Z d ¯ r | Ψ( ¯r , ¯ t ) | . (A10)Note that, in terms of the trial wave function, these de-pend on | Ψ( ¯r , ¯ t ) | or its square. If we form Ψ ∗ Ψ with5our trial wave function and neglect the rapidly oscillat-ing terms, we obtain | Ψ | ≈ N c N c X j =1 A j (¯ t ) exp ( − X η = x,y,z (¯ η − ¯ η j ) ¯ w j ) . (A11)It is clear that | Ψ | only depends on the centers, x , andwidths, w , if rapidly oscillating terms can be neglected.This is evidently also the case for | Ψ | which is the squareof | Ψ | . Thus, clearly, we have ¯ L (3 D )3 = ¯ L (3 D )3 ( x , w ) and¯ L (3 D )4 = ¯ L (3 D )4 ( x , w ), that is, these quantities only de-pend on the Gaussian center and width variational pa-rameters.Thus we introduce here the following variational poten-tials. First we define the “external” variational potentialthat accounts for the effects of V ext :¯ U (3 D )ext ( x , w ) ≡ N c ¯ L (3 D )3 ( x , w ) . (A12)We further introduce the “interaction” variational poten-tial that accounts for atom–atom scattering interactions:¯ U (3 D )int ( x , w ) ≡ N c ¯ L (3 D )4 ( x , w ) . (A13)Note that ¯ U (3 D )ext ( x , w ) depends on the particular form of V ext ( r , t ). The expression for ¯ U (3 D )int ( x , w ) is always thesame. The derivation of ¯ U (3 D )int ( x , w ) is more complexthan that of the other Lagrangian terms and is givenin Appendix B. It is only necessary to know that thesepotentials depend on x and w in order to derive theequations of motion for the Gaussian centers and widths.The final equations of motion will be written in terms ofthe gradients of ¯ U (3 D )ext ( x , w ) and ¯ U (3 D )int ( x , w ).
4. Derivation of ¯ L (3 D )5 The term ¯ L (3 D )5 accounts for the possibility that thesystem is in a rotating frame. Here we will assume thatthere is rotation only about the z axis. The ¯ L (3 D )5 termof the Lagrangian then has the form¯ L (3 D )5 = i ¯Ω z Z d ¯ r Ψ (cid:0) ¯ y Ψ ∗ x − ¯ x Ψ ∗ y (cid:1) . (A14)This term can be straightforwardly evaluated by insert-ing the expressions for the trial wave function and itsspace derivative into Eq. (A14), neglecting the oscillat-ing terms, evaluating the resulting Gaussian integrals,and eliminating the A j (¯ t ) factors. The result is¯ L (3 D )5 = ¯Ω z N c N c X j =1 " ¯ y j (cid:0) ¯ α jx + ¯ k jx + 2 ¯ β jx ¯ x j (cid:1) − ¯ x j (cid:0) ¯ α jy + ¯ k jy + 2 ¯ β jy ¯ y j (cid:1) (cid:21) (A15)
5. The 3D Lagrangian
Combining Eqs. (A4), (A8), (26), (27), and (A15) wehave the full 3D Lagrangian¯ L (3 D ) = 1 N c N c X j =1 X η = x,y,z " ˙¯ α jη ′ ¯ η j + ˙¯ β jη (cid:18) ¯ η j + 12 ¯ w jη (cid:19) + 12 ¯ w jη + 2 ¯ β jη ¯ w jη + (cid:16) β jη ¯ η j + ¯ α jη + ¯ k jη (cid:17) + ¯Ω z N c N c X j =1 " ¯ y j (cid:0) ¯ α jx + ¯ k jx + 2 ¯ β jx ¯ x j (cid:1) − ¯ x j (cid:0) ¯ α jy + ¯ k jy + 2 ¯ β jy ¯ y j (cid:1) + 12 N c (cid:16) ¯ U (3 D )ext ( x , w )+ ¯ U (3 D )int ( x , w ) (cid:17) (A16)We will use this form of the Lagrangian to derive the 3Dequations of motion. Appendix B: Derivation of ¯ U (3 D ) int ( x , w ) and ¯ U (1 D ) int ( x , w ) Here we present the derivation of the expressionsfor ¯ U (3 D )int ( x , w ) and ¯ U (1 D )int ( x , w ). The expression for¯ U (3 D )int ( x , w ) is¯ U (3 D )int ( x , w ) ≡ N c (cid:18)
12 ¯ gN Z d ¯ r | Ψ | (cid:19) (B1)In order to perform this integral we must first calculate | Ψ | . We can write this quantity as the square of | Ψ | : | Ψ | = 1 N c " N c X j =1 A j e P η = x,y,z Re { f j η } + N c X j ,j =1 j = j A j A j × exp ( X η = x,y,z (cid:16) f ∗ j η + f j η + i (cid:0) ¯ k j η − ¯ k j η (cid:1) ¯ η (cid:17)) where the have written the | Ψ | appearing inside the largesquare brackets as a sum of two terms. The first term isa sum over terms that definitely do not oscillate in space( j = j ) and the second term, where j = j , all of theterms definitely do oscillate in space.When this expression is expanded by writing the termsof the overall square only the square of the first term andthe square of the last term need be retained. We onlyneed to integrate terms that definitely do not oscillate asthe oscillatory terms will integrate to zero. Thus we can6rewrite the above expression as follows: | Ψ | ≈ N c " N c X j ,j =1 A j A j × exp ( X η = x,y,z (cid:16) f ∗ j η + f j η + f ∗ j η + f j η (cid:17)) + X j ,j =1 j = j X j ′ ,j ′ =1 j ′ = j ′ A j A j A j ′ A j ′ × exp ( X η = x,y,z (cid:16) f ∗ j η + f j η + f ∗ j ′ η + f j ′ η + i (cid:0) ¯ k j η − ¯ k j η + ¯ k j ′ η − ¯ k j ′ η (cid:1) ¯ η (cid:17)) The second term above still has oscillating terms. Weonly want to keep the non–oscillating terms. The termsthat don’t oscillate are those where j ′ = j and j ′ = j (since j = j and j ′ = j ′ are excluded already). Thus wecan evaluate the primed sums keeping only those termswhere j ′ = j and j ′ = j : | Ψ | ≈ N c " N c X j ,j =1 A j A j × exp ( X η = x,y,z (cid:16) f ∗ j η + f j η + f ∗ j η + f j η (cid:17)) + X j ,j =1 j = j A j A j × exp ( X η = x,y,z (cid:16) f ∗ j η + f j η + f ∗ j η + f j η (cid:17)) The two double sums appearing in the expression aboveare the same with one exception: the second double sumrequires j = j while the first double sum has no suchrestriction. We can write the first sum as two terms:a single sum where j = j and a double sum where j = j . This second term will then be identical to thesecond term in the original expression. Carrying out thisprocedure writing | Ψ | in terms of coordinates yields | Ψ | ≈ N c " N c X j =1 A j exp ( − X η = x,y,z (cid:16) η − ¯ η j ) ¯ w j η (cid:17)) + 2 X j ,j =1 j = j A j A j × exp ( − X η = x,y,z (cid:16) (¯ η − ¯ η j ) ¯ w j η + (¯ η − ¯ η j ) ¯ w j η (cid:17)) This form for | Ψ | can now be inserted into Eq. (B1)and then resulting Gaussian integrals can be straightfor- wardly evaluated. The final result is¯ U (3 D )int = ¯ gN (2 π ) / N c " N c X j =1 (cid:18) w j x ¯ w j y ¯ w j z (cid:19) + 2 / N c X j ,j =1 j = j Y η = x,y,z exp (cid:26) − ( ¯ η j − ¯ η j ) ¯ w j η + ¯ w j η (cid:27)(cid:0) ¯ w j η + ¯ w j η (cid:1) / (B2)Using a similar derivation, we obtain the 1D version of¯ U int : ¯ U (1 D )int = (¯ gN )(2 π ) / N c " N c X j =1 w j + 2 / N c X j ,j =1 j = j exp (cid:26) − ( ¯ x j − ¯ x j ) ¯ w j + ¯ w j (cid:27)(cid:0) ¯ w j + ¯ w j (cid:1) / (B3)The derivatives of these expressions with respect to thecenter and width variational parameters are computedstraightforwardly. Appendix C: The LVM Equations of Motion
In this section we derive the LVM equations of motionfor the variational parameters appearing in the 3D and1D trial wave functions given above. The major steps inthis undertaking are the following. First we write downthe equation of motion for each variational parameter, ¯ q ,using the standard Euler–Lagrange (EL) equation (Eq.(6)). We will show that further manipulations of theseequations can produce a closed set of second–order dif-ferential equations in time for the center position andwidth coordinates. Furthermore we will show that allof the other variational parameters can be expressed interms of the position and width coordinates and theirtime derivatives. We begin with the 3D case.
1. 3D LVM equations of motion a. 3D cloud center EOMs
To obtain the equations of motion for the cloud–centercoordinates we need the EOMs associated with the ¯ α jη and the ¯ η j . We begin with the EOMs associated withthe ¯ α jη . The Euler–Lagrange equations for these are dd ¯ t (cid:18) ∂ ¯ L (3 D ) ∂ ˙¯ α jη (cid:19) − ∂ ¯ L (3 D ) ∂ ¯ α jη = 0 , j = 1 , . . . , N c , η = x, y, z. (C1)7We can compute the derivatives of ¯ L using Eq. (25). Thestraightforward result is˙ x j = 2 (cid:0) β jx ¯ x j + ¯ α jx + ¯ k jx (cid:1) + ¯Ω z ¯ y j , (C2)˙ y j = 2 (cid:0) β jy ¯ y j + ¯ α jy + ¯ k jy (cid:1) − ¯Ω z ¯ x j , (C3)˙ z j = 2 (cid:0) β jz ¯ z j + ¯ α jz + ¯ k jz (cid:1) , j = 1 , . . . , N c . (C4)Next we obtain the EL equations associated with ¯ x j ,¯ y j , and ¯ z j . The ¯ η j Euler–Lagrange EOMs can be writtenas ∂ ¯ L∂ ¯ η j = dd ¯ t (cid:18) ∂ ¯ L∂ ˙¯ η j (cid:19) ,∂ ¯ L∂ ¯ η j = 0 , j = 1 , . . . , N c η = x, y, z (C5)where the second equality holds as there are no ˙¯ η termsin the Lagrangian. Taking the necessary derivatives weobtain the EL equations associated with ¯ j , ¯ y j , and ¯ z j ,respectively0 = h ˙ α jx + 2 ˙¯ β jx ¯ x j + 2 (cid:0) β jx ¯ x j + ¯ α jx + ¯ k jx (cid:1) (cid:0) β jx (cid:1) + 2 ¯Ω z ¯ β jx ¯ y j i + 12 ∂ ¯ U (3 D )ext ∂ ¯ x j + ∂ ¯ U (3 D )int ∂ ¯ x j ! − ¯Ω z (cid:0) β jy ¯ y j + ¯ α jy + ¯ k jy (cid:1) (C6)0 = h ˙ α jy + 2 ˙¯ β jy ¯ y j + 2 (cid:0) β jy ¯ y j + ¯ α jy + ¯ k jy (cid:1) (cid:0) β jy (cid:1) − z ¯ β jy ¯ x j i + 12 ∂ ¯ U (3 D )ext ∂ ¯ y j + ∂ ¯ U (3 D )int ∂ ¯ y j ! + ¯Ω z (cid:0) β jx ¯ x j + ¯ α jx + ¯ k jx (cid:1) (C7)0 = h ˙ α jz + 2 ˙¯ β jz ¯ z j + 2 (cid:0) β jz ¯ z j + ¯ α jz + ¯ k jz (cid:1) (cid:0) β jz (cid:1) i + 12 ∂ ¯ U (3 D )ext ∂ ¯ z j + ∂ ¯ U (3 D )int ∂ ¯ z j ! j = 1 , . . . , N c (C8)We have written the above equations in a form for conve-nient use in deriving the final equations of motion. Equa-tions (C2), (C3), and (C4) can be used along with Eqs.(C6), (C7), and (C8) to derive second–order differentialequations that involve only the center and width coordi-nates.We illustrate how this can be done by deriving theequation for ¨¯ x j . If we differentiate both sides of the equa-tion for ˙¯ x j (Eq. (C2)) with respect to time, the resultingsecond–order equation will contain both ¨¯ x j and ˙¯ x j . Us-ing Eq. (C2) again to eliminate ˙¯ x j from this second–orderequation gives the following result.¨¯ x j = 2 h ˙¯ α jx + 2 ˙¯ β jx ¯ x j + 2 (cid:0) β jx ¯ x j + ¯ α jx + ¯ k jx (cid:1) (cid:0) β jx (cid:1) + 2 ¯Ω z ¯ β jx ¯ y j i + ¯Ω z ˙¯ y j = − ∂ ¯ U (3 D )ext ∂ ¯ x j + ∂ ¯ U (3 D )int ∂ ¯ x j ! + 2 ¯Ω z (cid:0) β jy ¯ y j + ¯ α jy + ¯ k jy (cid:1) + ¯Ω z ˙¯ y j (C9) where we have obtained the second line by noting thatthe quantity in square brackets appearing in the first lineis identical to the quantity in square brackets in Eq. (C6).The second equation can be further simplified by notingthat the quantity in the parenthesis appearing in the sec-ond term is also present in Eq. (C2). The second–orderequation of motion is thus written in the compact finalform ¨¯ x j = 2 ¯Ω z ˙¯ y j + ¯Ω z ¯ x j − ∂ ¯ U (3 D ) ∂ ¯ x j , where¯ U (3 D ) ( x , w ) ≡ ¯ U (3 D )ext ( x , w ) + ¯ U (3 D )int ( x , w ) . (C10)The corresponding equations for the other cloud–centercoordinates can be derived similarly.The resulting second–order equations of motion for thecenter coordinates of cloud j are as follows.¨¯ x j = 2 ¯Ω z ˙¯ y j + ¯Ω z ¯ x j − ∂ ¯ U (3 D ) ∂ ¯ x j , (C11a)¨¯ y j = − z ˙¯ x j + ¯Ω z ¯ y j − ∂ ¯ U (3 D ) ∂ ¯ y j , (C11b)¨¯ z j = − ∂ ¯ U (3 D ) ∂ ¯ z j , j = 1 , . . . , N c . (C11c)Note that these equations depend only on x , w and theirtime derivatives. We now turn to the equations for thecloud widths. b. 3D cloud width EOMs It is possible to derive a set of second–order equationsof motion for the Gaussian cloud widths similar to thatfor the cloud centers. As will be seen below, the centerequations and width equations form a closed system. Allof the other variational parameters can be expressed interms of the centers and widths and their time deriva-tives.To obtain the width equations we start with the ELequation for the ¯ β jη . As before, we will only derive theequation of motion for ¯ w jx to illustrate how the deriva-tion is carried out. The EL equation for ¯ β jx is dd ¯ t ∂ ¯ L∂ ˙¯ β jx ! − ∂ ¯ L∂ ¯ β jx = 0 , j = 1 , . . . , N c . (C12)Differentiating the 3D Lagrangian in Eq. (25) and insert-ing into the above equations gives2¯ x j ˙¯ x j + ¯ w jx ˙¯ w jx = (4¯ x j ) (cid:0) β jx ¯ x j + ¯ α jx + ¯ k jx (cid:1) + 2 ¯Ω z ¯ y j ¯ x j + 4 ¯ β jx ¯ w jx . (C13)If we replace ˙¯ x j on the left–hand–side with its expressiongiven in Eq. (C2) we obtain the following remarkablysimple result ˙¯ w jx = 4 ¯ β jx ¯ w jx . (C14)8This result holds for ¯ y and ¯ z as well and enables us toexpress the ¯ β jη in terms of the widths.To proceed we turn to the Euler–Lagrange equationfor ¯ w jx which reads dd ¯ t (cid:18) ∂ ¯ L∂ ˙¯ w jx (cid:19) − ∂ ¯ L∂ ¯ w jx = 0 = ∂ ¯ L∂ ¯ w jx , (C15)where the last equality holds because no terms containing˙¯ w jx appear in the Lagrangian given in equation Eq. (25).Again taking derivatives of ¯ L (3 D ) we obtain4 ˙¯ β jx ¯ w jx + 16 ¯ β jx ¯ w jx = 4¯ w jx − ∂ ¯ U (3 D ) ∂ ¯ w jx . (C16)Now note that, if we differentiate both sides of Eq. (C14)with respect to time we obtain¨¯ w jx = 4 ˙¯ β jx ¯ w jx + 4 ¯ β jx ˙¯ w jx = 4 ˙¯ β jx ¯ w jx + 16 ¯ β jx ¯ w jx , (C17)where the second equality comes from reusing Eq. (C14)to replace ˙¯ w jx in the first line. Finally note that theright–hand–side of this last expression is identical tothe left–hand–side of Eq. (C16) and thus we obtain thesecond–order equation for ¯ w jx :¨¯ w jx = 4¯ w jx − ∂ ¯ U (3 D ) ∂ ¯ w jx . (C18)The derivation of equations for ¯ y and ¯ z are similar. Allthree equations can be written in the following compactform:¨¯ w jη = 4¯ w jη − ∂ ¯ U (3 D ) ∂ ¯ w jη , η = x, y, z j = 1 , . . . , N c . (C19)We can add these to EOMs for the cloud centers to getthe full set of equations of motion in 3D. They consist ofa pair of second–order ordinary differential equation forthe cloud centers and widths as well as expressions forthe ¯ β jη and the ¯ α jη in terms of the centers, widths andtheir first derivatives:¨¯ x j = 2 ¯Ω z ˙¯ y j + ¯Ω z ¯ x j − ∂ ¯ U (3 D ) ∂ ¯ x j , (C20a)¨¯ y j = − z ˙¯ x j + ¯Ω z ¯ y j − ∂ ¯ U (3 D ) ∂ ¯ y j , (C20b)¨¯ z j = − ∂ ¯ U (3 D ) ∂ ¯ z j , (C20c)¨¯ w jη = 4¯ w jη − ∂ ¯ U (3 D ) ∂ ¯ w jη , (C20d)¯ β jη = ˙¯ w jη w jη , (C20e)¯ α jx = ( ˙¯ x j − ¯Ω z ¯ y j ) − β jx ¯ x j − ¯ k jx , (C20f)¯ α jy = ( ˙¯ y j + ¯Ω z ¯ x j ) − β jy ¯ y j − ¯ k jy , (C20g)¯ α jz = ˙¯ z j − β jx ¯ x j − ¯ k jx , (C20h) η = x, y, z j = 1 , . . . , N c The equations for the cloud centers and cloud widths(Eqs. (C20a), (C20b), (C20c), and (C20d)) form a closedset that contain only the ¯ η j , ˙¯ η j , ¯ w jη , and ˙¯ w jη . Oncethese quantities are obtained, all of the other variationalparameters can be calculated. Appendix D: Details of the GPE simulation
The conditions assumed in the simulation shown inFig. 2 include a condensate of 10,000 Rb atoms whichwere confined in a harmonic trap with a frequency of ω T / π = 1 Hz. The condensate pieces were given aninitial velocity of v = 2 × − m/s and the wait timeafter reaching the trap turning points is T = 500 ms. Thesource mass is absent in the simulation presented.The GPE was solved numerically using the split–step,Crank–Nicolson algorithm [27]. Briefly, in this algorithm,the Hamiltonian is split as ˆ H = ˆ T + ˆ V where ˆ T is thekinetic energy and ˆ V is the potential plus the nonlinearinteraction term. The condensate wave function is ad-vanced from time t to time t + δt by multiplying ψ ( x, t )point by e − iV ( x,t ) δt/ ¯ h and then the result is multipliedby e − i ˆ T δt/ ¯ h . This second step is carried out using theCrank–Nicolson algorithm. The time–dependent vari-ational equations of motion were solved by the Eulermethod.The comparison of the GPE and the LVM simulationsfor the initial split phase of the AI sequence is shown inFig. 2. The initial split itself divided into three parts. Inthe first part the condensate pieces fly out to the turningpoints in the presence of the harmonic potential whosefrequency is ω T / π = 1 Hz and finally stop after a quar-ter period (250 ms) as shown in Figs. 1(b)–(c). Both theGPE and the LVM densities are plotted together on thesame graph at four different times as seen in the leftpanel of Fig. 2. Each density consists of two symmetri-cally placed peaks and the left peak of each density islabeled with its time stamp. In the second part of theinitial split sequence the two condensate pieces evolve atrest for a wait time of T = 500 ms while the trap isoff. This evolution is shown in the middle panel of Fig. 2.Finally the condensate pieces recombine after the trap isturned back on. This is shown in the right panel of Fig. 2. Appendix E: Derivation of ¯ U (1 D ) ext ( x , w ) for anharmonic plus point–mass potential In this appendix we derive the variational potentialsfor ¯ U (1 D ) ext ( x , w ) for the case of an external potential con-sisting of an harmonic potential plus the gravatationalpotential of a point mass situated far from the conden-sate. The 1D version is used to obtain the equations ofmotion whose solutions are compared with the full 1DGPE simulation of the illustrative AI measurement of G in a microgravity environment presented in Section IV.9The 1D external variational potential is given by Eq.(28)¯ U (1 D )ext = N c X j =1 (cid:18) π / ¯ w j (cid:19) Z + ∞−∞ d ¯ x e − (¯ x − ¯ x j ) / ¯ w j ¯ V ext (¯ x, ¯ t ) . (E1)The first step is to find an approximate expression for¯ V ext (¯ x, ¯ t ) for the combination of an harmonic potentialplus a point mass.We assume that the harmonic trap is centered at theorigin of coordinates and that a point mass with mass M SM is located at x SM . The exact potential can thusbe written (in SI units) as V ext ( x ) = 12 M ω T,x x − GM M SM | x SM − x | ≡ V H ( x ) + V G ( x ) . (E2)We want to approximate V G ( r ) by assuming that thedistance of any cloud to the origin is much smaller thanthe distance of the source mass to the origin. We havechosen the origin to be at the center of the harmonicpotential confining the BEC.First we consider only the gravitational part of thepotential: V G ( x ) = − GM M SM | x SM − x | . (E3)we can approximate this exact expression by making aTaylor expansion about x = 0 to second order in x/x SM V G ( x ) ≈ − GM M SM | x SM | (cid:0) x SM + x SM x + x (cid:1) (E4)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 ≫ | x | . It will be convenient hereto introduce the gravitational frequency ω SM ≡ (cid:18) GM SM | x SM | (cid:19) / . (E5)Thus we can rewrite the approximate gravitational po-tential in terms of this quantity as follows. V G ( r ) ≈ − M ω SM (cid:0) x SM + x SM x + x (cid:1) ¯ V G ( r ) ≈ − ¯ ω SM (cid:0) ¯ x SM + ¯ x SM ¯ x + ¯ x (cid:1) , (E6)where, in the second line, we have expressed the gravi-tational potential in scaled units. Now we are ready towrite down the full external potential. The full externalpotential in scaled units is thus¯ V ext (¯ x ) = (cid:0) ¯ ω T − ω SM (cid:1) ¯ x − ¯ ω SM ¯ x SM ¯ x − ¯ ω SM ¯ x SM (E7)Inserting the above expression into Eq. (E1) yields thefinal expression for ¯ U (1 D )ext for the case of an external har-monic trap plus the gravitational potential produced bya point mass far away from the condensate is given by¯ U (1 D )ext ( x , w ) = N c X j =1 (cid:16) (cid:0) ¯ ω T − ω SM (cid:1) (cid:0) ¯ x j + ¯ w j (cid:1) − ¯ ω SM ¯ x SM ¯ x j (cid:17) − N c ¯ ω SM ¯ x SM ..