# Approximate mean–field equations of motion for quasi–2D Bose–Einstein condensate systems

###### Abstract

We present a method for approximating the solution of the three–dimensional, time–dependent Gross–Pitaevskii equation (GPE) for Bose–Einstein condensate systems where the confinement in one dimension is much tighter than in the other two. This method employs a hybrid Lagrangian variational technique whose trial wave function is the product of a completely unspecified function of the coordinates in the plane of weak confinement and a gaussian in the strongly confined direction having a time–dependent width and quadratic phase. The hybrid Lagrangian variational method produces equations of motion that consist of (1) a two–dimensional, effective GPE whose nonlinear coefficient contains the width of the gaussian and (2) an equation of motion for the width that depends on the integral of the fourth power of the solution of the 2D effective GPE. We apply this method to the dynamics of Bose–Einstein condensates confined in ring–shaped potentials and compare the approximate solution to the numerical solution of the full 3D GPE.

###### pacs:

03.75.Gg,67.85.Hj,03.67.Dg## I Introduction

Recent advances in laser–control technology have enabled the laboratory realization of Bose–Einstein condensate (BEC) systems subjected to all–optical potentials which provide strong confinement in a horizontal plane and an arbitrary potential within this plane. These potentials can be produced by a combination of a horizontal light sheet combined with a rapidly moving red– or blue–detuned vertical laser that “paints” an arbitrary time–averaged optical dipole potential in the horizontal plane Henderson et al. (2009). Horizontal light sheets can also be combined with vertically propagating beams in specialized laser modes, such as Laguerre–Gauss modes, to produce other types of novel potentials Ramanathan et al. (2011). In addition to providing strong vertical confinement and counteracting the effect of gravity, the light sheet provides stabilization against dynamic excitations of the condensate Onofrio et al. (2000) as well as thermal phase fluctuations Dettmer et al. (2001).

The ability to create and probe quasi–2D BECs in arbitray 2D potentials is motivated by several areas of current ultra–cold atom research. For example, condensates in toroidal and ring lattices can be studied. Stable states of multiple vortices and persistent currents can be created and studied by stirring the condensate Brand and Reinhardt (2001); Ryu et al. (2007); Ramanathan et al. (2011). There are proposals for creating ring lattices and for studying non–equilibrium phase transitions within this geometry Amico et al. (2005); Dziarmaga et al. (2008). Toroidal geometries are well–suited for studying topological defects that may appear during a rapid cooling process that produces a condensate Anglin and Zurek (1999); Weiler et al. (2008). There is some indication that stirring within a ring–lattice geometry can produce a coherent superposition of states with different circulation which can lead to a reduction in the threshold of the Mott–Insulator phase transition Rey et al. (2007). These systems also offer an excellent finite–sized testbed for systems of ultra–cold atoms that mimic condensed–matter systems Lewenstein et al. (2007).

Quasi–2D BECs may also provide a convenient platform for studying systems of ultra–cold neutral atoms that are analogs of electronic materials, devices, and circuits Seaman et al. (2007). Such systems are called “atomtronic” because strongly interacting Bose gases in a lattice potential are analogous to “electronic” systems of electrons moving in the periodic lattice potential of a crystalline solid. The ability to produce arbitrary potentials in the plane of the quasi–2D condensate may enable the controlled study of novel atomtronic systems. In particular it may be possible to produce circuit–like potentials within the plane.

The behavior of many of the above–mentioned ultra–cold bosonic systems can be described using mean–field theory. In this case, the governing equation is the time–dependent Gross–Pitaevskii equation (TDGPE) Gross (1961); Pitaevskii (1961). This is a partial differential equation in three space variables and one time variable whose solution represents the wave function of the single–particle orbital that all of the condensate atoms occupy. Experiments conducted on these systems typically involve releasing the condensate for imaging. In this case solution of the 3D TDGPE becomes a challenging numerical problem due to the volume that must be accounted for in simulating the experiment.

In this paper, we present a variational approximation to the solution of the TDGPE for these quasi–2D systems which produces equations of motion whose numerical solution can be obtained 100 to 1000 times faster than solving the full 3D TDGPE. This approximation is based on a variant of the standard Lagrangian Variational Method (LVM) Perez-Garcia et al. (1997) in which some of the variational parameters are functions of the space coordinates and time while others are only functions of time. The work presented here applies this “hybrid” version of the LVM to a quasi–2D systems of bosonic atoms. The hybrid LVM was previously applied to a quasi–1D system where only one dimension was weakly confined compared to the two other dimensions Edwards et al. (2005).

This paper is organized as follows. In Section II we describe the LVM and its hybrid form and derive the approximate equations of motion. We also derive the equations that provide the proper variational stationary solution in which a condensate is trapped in a confining potential. Section III presents a comparison of the solution of the hybrid LVM equations of motion with the numerical solution of the 3D TDGPE for a BEC confined in a ring–shaped potential. The parameters for this example system were taken from an actual experiment. Section IV presents a summary of the work.

## Ii Two–dimensional hybrid LVM equations of motion

The condensate wave function of a BEC that is strongly confined in one dimension ( direction) relative to the confinement in the other two dimensions ( plane) can often be approximated as the product of a function of and only with a gaussian function of only. In the mean–field approximation, the actual behavior of the condensate wave function is governed by the 3D, time–dependent, Gross–Pitaevskii equation. However, it is possible to find equations of motion from which the approximate product wave function can be constructed at each moment of time using the Lagrangian Variational Method. We briefly describe this method next.

### ii.1 The Lagrangian variational method

The Lagrangian Variational Method provides approximate solutions to the 3D TDGPE in the form of equations of motion for time–dependent parameters that appear in an assumed trial wave function. Thus, in the standard LVM, the exact solution of the 3D TDGPE requiring the solution of a partial differential equation in three space and one time variable is traded for the solution of ordinary differential equations in time for the variational parameters of a trial wave function of fixed functional form.

The 3D TDGPE can be written as

(1) |

where is the mass of a condensate atom, is the interaction strength of low–energy binary scattering events with being the –wave scattering length, is the number of atoms in the condensate, and is the external potential.

The TDGPE is itself a variational equation–of–motion and is derived from the following Lagrangian density:

(2) | |||||

where and . The associated Euler–Lagrange equation that produces the TDGPE with the above Lagrangian density is given by

(3) |

The LVM is an approximation method that produces an equation of motion for the time–dependent variational parameters, , appearing in a given trial wave function . The equations of motion for these parameters are obtained by inserting the trial wave function into the LVM Lagrangian density, integrating this over the spatial variables:

(4) |

and applying the usual Euler–Lagrange equations:

(5) |

This is the standard Lagrangian Variational Method Perez-Garcia et al. (1997).

The LVM can be regarded as having two limits in terms of the chosen trial wave function. The first limit consists of choosing a trial wave function where the variational “parameter” is . This choice enables the variational solution to vary in any possible way. As noted above, when Eq. (3) is applied to the Lagrangian density to derive the equation of motion, it turns out to be the full TDGPE. In the other limit, the trial wave function is chosen to have a fixed functional form of the spatial coordinates where the time dependence resides entirely within a set of variational parameters, , so that the Lagrangian depends only on these parameters, . The shape of this trial wave function can only be varied by changing the values of the . The equations of motion for the are ordinary differential equations in time and are obtained from the usual Euler–Lagrange equations, Eqs. (5). It is also possible to choose a “hybrid” trial wave function that plots a course midway between these two limits. We describe this approach now.

### ii.2 The hybrid LVM

The “Hybrid Lagrangian Varational Method” (HLVM) is an LVM in which the trial wave function consists of a completely unspecified function of some of the spatial coordinates, , multiplied by a fixed function of the rest of the coordinates that also contains some time–dependent variational parameters, . The HLVM is expected to apply to systems where there is tight confinement in one or two dimensions. The coordinates appearing in are those for which the confinement is weak while the trial wave function is assumed to be gaussian in the coordinates of tight confinment. The HLVM for tight confinement in two dimensions has been studied earlier Edwards et al. (2005). Here we study the case where there is tight confinement in one dimension only.

Coupled equations of motion can be derived from a “hybrid” Lagrangian which is constructed by integrating the Lagrangian density in Eq. (2) over the space coordinate(s) of the tightly confined direction(s). The resulting hybrid Lagrangian can be used to derive coupled equations of motion for both and the set of variational parameters .

Before proceeding with the derivation of these equations of motion, we will first introduce scaled variables and rewrite the LVM equations in terms of these variables. Scaled units are referenced to a chosen unit of length, denoted by , and scaled spatial coordinates are given by

(6) |

Energy and time units are defined in terms of enabling the definition of a scaled time:

(7) |

Hereafter barred symbols will denote quantities expressed in their appropriate scaled units. It will also be convenient to express the solution of the 3D TDGPE in terms of scaled units:

(8) |

In terms of these variables the TDGPE becomes

(9) |

where . In scaled units, the Lagrangian density takes the form:

(10) | |||||

and the scaled Euler–Lagrange equation becomes:

(11) |

Now we turn to the description of the hybrid Lagrangian Variational Method.

In deriving the HLVM equations of motion we will assume that the trapping potential can be written (at least approximately) as the sum of a part that depends only on the loosely confined coordinates (here and ) and a part that is harmonic in the tightly bound direction. Under this assumption we can write the potential as:

(12) |

This form of the potential applies in many realistic experimental cases such as the painted potentials mentioned earlier.

The trial wave function for the HLVM equations of motion is written as follows:

(13) |

Here the trial wave function is a product of a completely unspecified function with a gaussian function having a time–dependent width and quadratic phase coefficient . These are the variational parameters that will appear in the HLVM equations of motion.

The first step in the hybrid LVM consists of constructing a “hybrid” Lagrangian by integrating only over the spatial coordinate along which the system is strongly confined:

(14) |

The resulting hybrid Lagrangian is given by

In the above we have eliminated the variational parameter using the normalization constraint:

and by requiring that the separate parts of the product wave function to be separately normalized to unity:

(16) |

The second step in the HLVM is to apply the Euler–Lagrange equations of motion to to obtain the equations of motion. The equation for is a modified version of Eq. (11):

(17) |

and the Euler–Lagrange equations for and are the usual ones:

(18) |

Applying Eq. (17) yields the following equation for :

(19) | |||||

where

(20) |

This seemingly complicated function of can be transformed away by defining

(21) |

Inserting this into the equation of motion for yields an effective 2D Gross–Pitaevskii–like equation for :

Applying the Euler–Lagrange equation for gives the following result:

(23) |

We can obtain a simplified equation of motion by integrating both sides of the above over all and :

It is easy to show that the integral on the left is zero by using the equation of motion for . The integral on the right is unity by normalization and so we obtain the following relationship between and , :

(25) |

Thus, if and are known, is determined.

Applying the Euler–Lagrange equation for gives

(26) |

Integrating this equation over all on both sides as before we obtain

(27) |

where

(28) |

Note that we have used Eq. (21) to replace with .

It is possible to eliminate from the above equation by differentiating both sides of Eq. (25) with respect to time. We obtain

(29) |

where the second equality results from using Eq. (25) to replace with . Now we see that the right–hand–side of the above equation is identical to the first two terms on the left–hand–side of Eq. (27). Thus we can rewrite this equation as follows:

(30) |

This is the final equation of motion for .

### ii.3 The HLVM equations of motion and the variational initial state

The full set of HLVM equations of motion consist of a 2D effective GP–like equation for :

and an equation for :

(32) |

These two equations form a closed system from which , , and can be obtained. From these, the value of and can be calculated:

(33) |

Using these quantities, the full value of the variational trial wave function can be constructed:

(34) | |||||

where

(35) |

Note that Eqs. (LABEL:hlvm_phi_eq) and (32) are coupled. The nonlinear term in the 2D GPE for contains the gaussian width, , while the equation for contains the factor which is the integral of the fourth power of .

The final element required for this method to be used as a means to find an approximation to the solution of the 3D TDGPE is a set of initial conditions. We present one possibility here based on the physics of Bose–Einstein condensate systems.

In a typical BEC experiment a condensate is formed in an atom trap. If no further changes in the condensate’s environment occur, the condensate wave function should then, in principle, only acquire an overall time–dependent phase as it evolves in time. In the HLVM this situation should therefore be represented by the stationary solution of the above equations. We denote this stationary solution as

(36) |

This solution satifies the following time–independent equations.

(37) | |||||

and

(38) |

where

(39) |

The factor in the equation for is the chemical potential of the initial condensate. These equations for the stationary variational solution must be solved self–consistently.

## Iii Comparison with 3D GPE

In this section we illustrate the ability of the HLVM equations of motion to approximate the exact solution of the 3D TDGPE by comparing the two solutions for a case of current interest. The system we will consider is that of a Bose–Einstein condensate of Na atoms confined in a ring–shaped potential under the same conditions as in a recent experiment Ramanathan et al. (2011) conducted at NIST. To simplify the analysis we will only compute the profile of the condensate density integrated along the vertical direction for each point in the plane of weak confinement. This quantity predicted by the TDGPE will be compared with that predicted by the HLVM equations of motion. This is the quantity that can be compared with experiment.

In the NIST experiment a vertical Laguerre–Gauss (LG) laser beam (LG) was intersected with a horizontal light sheet. The shape of the vertical LG beam was approximately a hollow cylinder with thick walls so that its intersection with the horizontal light sheet created a ring–shaped region of maximum light intensity. Tuning the frequency of the beams to the red of the lowest electronic transition created a potential that caused the atoms to seek the maximum intensity.

In this comparison, we simulate an experiment in which a condensate is created in this ring potential, optionally stirred, and then probed. We simulate two types of probes: (1) direct release of the condensate by turning off all trapping potentials after stirring, and (2) release of the condensate after the ramp down of the Laguerre–Gauss potential. The stirring, which adds units of angular momentum to the condensate, is simulated by phase imprint. That is, the initial condensate wave function is multiplied by where is the azimuthal angle around the vertical axis and is an integer.

In each case we will compare what would be the measured density profile, as predicted by the 3D GPE and by the HLVM, for different times during the ramp down or expansion where the value of is fixed, and also for a fixed final time–of–flight for a range of different values. For maximum clarity, we present the two density profiles as a plot of the density along a line that cuts through the center of the ring. Since all of the density profiles are cylindrically symmetric, these plots will convey all of the available density information.

In these simulations the trap potential is modeled as the sum of a Laguerre–Gauss optical potential Wright et al. (2000) plus a vertical gaussian due to the light sheet. This potential can be written (in scaled units) as:

where the factor is included so that becomes the depth of the potential due to the Laguerre–Gaussian beam and is the radial position of its minimum. A time–dependent, dimensionless turn–on function, , is inserted to simulate the ramp down of the Laguerre–Gauss potential. The factor is the depth of the potential due to the light sheet and the –dependent gaussian light–sheet potential has been approximated by an harmonic oscillator.

Both the 3D TDGPE and the 2D GPE part of the HLVM equations of motion were solved using the split–step, Crank-Nicolson method. The 3D TDGPE was solved on a grid in which there were 400 points along and and 200 points along . The 2D GPE part of the HLVM equations of motion was solved on a grid of 800 points along both and . The codes that were used to solve these equations were extensively modified versions of codes publicly available in the literature Muruganandam and Adhikari (2009). The initial condensate wave function for the 3D TDGPE was obtained by solving it in imaginary time. Initial conditions for the HLVM equation of motion were obtained by solving equations (37) and (38) self consistently as follows. First, a value for was chosen, the associated was then found by integrating Eq. (LABEL:hlvm_phi_eq) in imaginary time, next the value of was calculated which was then used to compute the value of in Eq. (38). The value of was incremented and the process was repeated to compute a new value of . This process was continued until a root of was found. The value of for this case is the self–consistent gaussian width of the stationary solution of the HLVM equations of motion.

In the direct–release process simulated, the number of condensate atoms was 750,000 atoms and the scattering length of Na was taken to be 53 bohr. The minimum of the LG potential was set at m. The depth of the LG potential was taken to be nK which is equivalent (via ) to a radial harmonic frequency of Hz. The frequency of the harmonic light sheet potential was taken as Hz and the light–sheet depth was nK although this last quantity makes no difference in the shape of the initial–state density.

Figure 1 displays a comparison of the integrated column density of a released ring BEC predicted by the 3D TDGPE with that predicted by the HLVM equations of motion at six different times–of–flight (TOF) after release beginning with Fig. (1a) showing the moment of release. The condensate has been stirred so that it is released having one unit of angular momentum. We note that there is good agreement with quantitative differences occuring in the heights of individual peaks and in the position of the peaks at later times. The comparisons are typical of a variational solution in that they are the “best fit” to the exact solution for the given trial wave function.

A comparison of the GPE and LVM results for a ring BEC directly released from the trap and allowed to expand for a fixed TOF for different initial angular momenta is exhibited in Fig. 2. The figure displays comparisons for values ranging from 0 to 5. Again the agreement is good although there are some quantitative differences as to positions of the individual peaks. It is clear that there is qualitative and almost–quantitative agreement between the 3D TDGPE and the HLVM equations of motion for these cases.

We next compare the results of the TDGPE and HLVM for ring–BEC evolution while the LG potential is ramped down from its initial value. This differs from the previous comparison in that the confining light–sheet potential remains unchanged during the ramp down. In the simulated ramp–down process the number of condensate atoms was 500,000 atoms. The LG potential depth was ramped linearly down from its initial value of nK (the same as previously) to 20% of this value over a timespan of 50 ms. The light sheet potential was the same as in the direct–release simulations. A phase imprint was applied to the condensate to simulate one unit of angular momentum added by stirring.

Figure 3 shows the comparison starting at ms and for every 10 ms thereafter until the rampdown is complete. The two solutions again exhibit the type of agreement that is usual for variational approximations in that the variational solution is a “best fit” to the numerical solultion. With that caveat, the agreement here is quite good during the entire ramp down.

## Iv Summary

In this paper we derived equations of motion whose solution approximates the solution of the 3D TDGPE applied to a quasi–2D Bose–Einstein condensate. The equations were derived using a hybrid Lagrangian Variational Method. Similar equations of motion were derived earlier Edwards et al. (2005) for quasi–1D BEC systems. The main advantage of solving these equations is that numerical solution of the HLVM equation can be performed 100 to 1000 times faster than solving the 3D TDGPE. In the comparison simulation presented in Section III, solving the realtime 3D TDGPE required more than 24 hours of CPU time while solving the HLVM equation took about 10 minutes on a commodity desktop PC. The resulting speedup here is roughly a factor of 150.

This advantage enables rapid simulation of many different possible quasi–2D systems. It should be noted that when the HLVM was applied to a quasi–1D system describing soliton splitting in optical fibers Infeld et al. (2006), it was found that occasionally the HLVM equations of motion did not provide any advantage over the regular LVM technique. In the work cited, the authors recommended that any important results that come out of HLVM simulations be confirmed by simulations using the full equations. We did not find any case where the HLVM equations predicted behavior that was qualitatively different from the 3D TDGPE. However, we agree with the recommendation of the authors of Ref. Infeld et al. (2006).

This caveat notwithstanding, the HLVM equations derived in this paper enable rapid study of different systems of current experimental interest. In particular, they should be useful in simulating time–dependent behavior of quasi–2D atomtronic systems where mean–field theory applies. We expect this approximation to become a useful tool in studying future quasi–2D Bose–Einstein condensate systems.

###### 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. The authors acknowledge helpful discussions from Kevin Wright, Gretchen Campbell, and Noel Murray.## References

- Henderson et al. (2009) K. Henderson, C. Ryu, C. MacCormick, and M. G. Boshier, New Journal of Physics 11, 043030 (2009).
- Ramanathan et al. (2011) A. Ramanathan, K. C. Wright, S. R. Muniz, M. Zelan, W. T. Hill, C. J. Lobb, K. Helmerson, W. D. Phillips, and G. K. Campbell, Phys. Rev. Lett. 106, 130401 (2011).
- Onofrio et al. (2000) R. Onofrio, D. S. Durfee, C. Raman, M. Köhl, C. E. Kuklewicz, and W. Ketterle, Phys. Rev. Lett. 84, 810 (2000).
- Dettmer et al. (2001) S. Dettmer, D. Hellweg, P. Ryytty, J. J. Arlt, W. Ertmer, K. Sengstock, D. S. Petrov, G. V. Shlyapnikov, H. Kreutzmann, L. Santos, et al., Phys. Rev. Lett. 87, 160406 (2001).
- Brand and Reinhardt (2001) J. Brand and W. P. Reinhardt, J. Phys. B: At. Mol. Opt. Phys. 34, L113 (2001).
- Ryu et al. (2007) C. Ryu, M. F. Andersen, P. Cladé, V. Natarajan, K. Helmerson, and W. D. Phillips, Phys. Rev. Lett. 99, 260401 (2007).
- Amico et al. (2005) L. Amico, A. Osterloh, and F. Cataliotti, Phys. Rev. Lett. 95, 063201 (2005), URL http://link.aps.org/doi/10.1103/PhysRevLett.95.063201.
- Dziarmaga et al. (2008) J. Dziarmaga, J. Meisner, and W. H. Zurek, Phys. Rev. Lett. 101, 115701 (2008), URL http://link.aps.org/doi/10.1103/PhysRevLett.101.115701.
- Anglin and Zurek (1999) J. R. Anglin and W. H. Zurek, Phys. Rev. Lett. 83, 1707 (1999), URL http://link.aps.org/doi/10.1103/PhysRevLett.83.1707.
- Weiler et al. (2008) C. N. Weiler, T. W. Neely, D. R. Scherer, A. S. Bradley, M. J. Davis, and B. P. Anderson, Nature 455, 948 (2008).
- Rey et al. (2007) A. M. Rey, K. Burnett, I. I. Satija, and C. W. Clark, Phys. Rev. A 75, 063616 (2007), URL http://link.aps.org/doi/10.1103/PhysRevA.75.063616.
- Lewenstein et al. (2007) M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A. Sen, and U. Sen, Advances in Physics 56, 243 (2007).
- Seaman et al. (2007) B. T. Seaman, M. Krämer, D. Z. Anderson, and M. J. Holland, Phys. Rev. A 75, 023615 (2007), URL http://link.aps.org/doi/10.1103/PhysRevA.75.023615.
- Gross (1961) E. Gross, Il Nuovo Cimento 20, 454 (1961).
- Pitaevskii (1961) L. Pitaevskii, Soviet Physics JETP 13, 451 (1961).
- Perez-Garcia et al. (1997) V. M. Perez-Garcia, H. Michinel, J. I. Cirac, M. Lewenstein, and P. Zoller, Phys. Rev. A 56, 1424 (1997).
- Edwards et al. (2005) M. Edwards, L. M. DeBeer, M. Demenikov, J. Galbreath, T. J. Mahaney, B. Nelsen, and C. W. Clark, Journal of Physics B: Atomic, Molecular and Optical Physics 38, 363 (2005).
- Wright et al. (2000) E. M. Wright, J. Arlt, and K. Dholakia, Phys. Rev. A 63, 013608 (2000).
- Muruganandam and Adhikari (2009) P. Muruganandam and S. Adhikari, Computer Physics Communications 180, 1888 (2009).
- Infeld et al. (2006) E. Infeld, M. Matuszewski, and M. Trippenbach, Journal of Physics B: Atomic, Molecular and Optical Physics 39, L113 (2006).