Collisions of anisotropic two-dimensional bright solitons in dipolar Bose-Einstein condensates
Rüdiger Eichler, Damir Zajec, Patrick Köberle, Jörg Main, Günter Wunner
CCollisions of anisotropic two-dimensional bright solitons in dipolar Bose-Einsteincondensates
R¨udiger Eichler, ∗ Damir Zajec, † Patrick K¨oberle, J¨org Main, and G¨unter Wunner
Institut f¨ur Theoretische Physik 1, Universit¨at Stuttgart, 70550 Stuttgart, Germany (Dated: August 20, 2018)We investigate the coherent collision of anisotropic quasi-two-dimensional bright solitons in dipolarBose-Einstein condensates. Our analysis is based on the extended Gross-Pitaevskii equation, and weuse the split-operator method for the grid calculations and the time-dependent variational principlewith an ansatz of coupled Gaussian functions to calculate the time evolution of the ground state. Wecompare the results of both approaches for collisions where initially the solitons are in the repellingside-by-side configuration and move towards each other with a specific momentum. We change therelative phases of the condensates, and introduce a total angular momentum by shifting the solitonsin opposite direction along the polarization axis. Our calculations show that collisions result inbreathing-mode-like excitations of the solitons.
PACS numbers: 03.75.-b, 05.45.-a, 67.85.-d, 34.50.-s
I. INTRODUCTION
Bose-Einstein condensates (BECs) of magnetic atomshave attracted much attention since their experimentalrealization with Cr atoms [1]. Recently, the creation ofcondensates of
Dy [2] and
Er [3] atoms with muchlarger magnetic moments than Cr have also been re-ported. Furthermore, there has been fast progress to-wards the condensation of polar molecules with electricdipole moments [4], where the dipole-dipole interaction(DDI) is even more dominant. A review of the physics ofdipolar bosonic quantum gases has recently been given byLahaye et al. [5]. The features of the DDI being a non-local long-ranged and anisotropic interaction give rise toa variety of new effects. One example is the creation ofsolitary waves, where in analogy to nonlinear optics theeffects of dispersion and nonlinearity may cancel eachother. This leads to a condensate with a shape constantin time. The experimental realization of one-dimensionalsolitons in self-attractive BECs of Li atoms has beenreported [6]. Tikhonenkov et al. have theoretically pre-dicted 2D solitons [7] and K¨oberle et al. have proposed arealistic experimental setup for the creation of a 2D soli-ton [8]. An exciting aspect of multidimensional solitonsis their anisotropic nature, based on the in-plane polar-ization of the dipoles of such solitons. 2D solitons havealready been studied using a variational ansatz with asingle Gaussian and with coupled Gaussian functions [9].Adhikari et al. have recently investigated axially sym-metric and vortex solitons on a one-dimensional opticallattice [10]. Note that in contrast to systems with har-monic traps, where the density distribution in the trapdirection is an approximate Gaussian, systems in an op-tical lattice will have an exponential density distribution. ∗ [email protected] † [email protected](Both authors contributed equally to this work) The collision of axially symmetric bright 2D solitonshas been studied by Pedri et al. [11] and Adhikari etal. [10]. Pedri et al. investigated a system with dipolesaligned parallel to the harmonic trap, while Adhikari etal. used an optical lattice instead. In both cases, thesign of the DDI has to be inverted by fast rotation of theorientation of the dipoles [12]. The resulting interactionenergy becomes U d ( R ) = − α (3 cos ϑ − / R , where ϑ is the angle between the polarization axis and R = r − r (cid:48) .The factor α can continuously be changed from − / et al. [13] have investigated the collision ofone-dimensional bright and vortex solitons. The inves-tigations in [7–9] concentrated on the creation and thestability of 2D solitons with respect to small perturba-tions. However, one important property of solitons isthat their shape is constant in time even when they aremoving. Therefore, the collision of two solitons is an ad-equate scenario for the investigation of soliton dynamicsfar beyond small excitations. The influence of the non-linear contact interaction and the DDI are of particularinterest in such calculations.As mentioned above, the creation of a BEC of magneticatoms has been realized with a variety of species. Ourresults are valid for all dipolar systems, but we will addthe corresponding values for a system with 20 000 Cr-atoms per soliton in parentheses.At sufficiently low temperatures, the dynamics of aBose-Einstein condensate can be described by the ex-tended Gross-Pitaevskii equation (GPE) which in atomicunits and with particle-number scaling [14] reads a r X i v : . [ c ond - m a t . qu a n t - g a s ] O c t H ( t )Ψ( r , t ) = ( − ∆ + V har + V sc + V d ) Ψ( r , t )=i ∂ t Ψ( r , t ) , (1)with V har = γ y y , V sc = 8 πa | Ψ( r , t ) | ,V d = (cid:90) d r (cid:48) − ϑ | r − r (cid:48) | | Ψ( r (cid:48) , t ) | . Here a is the scattering length and Ψ designates themean-field wave function. The dipoles are polarized alongthe z -axis, so that ϑ is the angle between the z -axis andthe vector r − r (cid:48) . We choose the y -direction as the axis ofconfinement perpendicular to the polarization axis where γ y = 20 000 (420 Hz), while the condensate is free in x -and z -direction. All simulations deal with condensates oflow densities, and only a small period of time in which thetwo condensates merge to one transient condensate withhigher density. This means that we do not need to take athree-body-loss term [8] into account, as the resulting ab-sorption images ( | ψ | integrated along the y -axis) wouldonly be slightly affected. We checked this assumption forthe calculation of the collision without difference in phaseand without angular momentum which up to the time of t = 0 .
06 ( t = 0 .
001 corresponds to 15 ms) resulted onlyin a loss of about 5 .
5% of the particles.As has been shown in [9], solitons only exist in a cer-tain range of values of the scattering length, which canbe tuned by the use of Feshbach resonances [15]. For toolarge values, the condensate will disperse, while too smallvalues lead to the collapse of the condensate. In the fol-lowing the scattering length is chosen to be 0 .
14 (12 . a B ,where a B is the Bohr radius) . II. NUMERICAL APPROACH
The main theoretical task for the grid calculations ishow to apply the time evolution operator U = e − iHt ona state | ψ (cid:105) . For this, one splits U symmetrically by usingthe Baker-Campell-Hausdorff formula [16] U (∆ t ) = e − iH ∆ t = e − i ( T + V )∆ t ≈ e − i T ∆ t e − iV ∆ t e − i T ∆ t , (2)where V = V har + V sc + V d . One projects the action ofthe approximated time evolution operator on the basis ofthe position operator and makes use of the possibility toinsert (cid:82) d ν | ν (cid:105) (cid:104) ν | = 1: ψ ( r , t + ∆ t ) = (cid:104) r | U (∆ t ) | ψ (cid:105) = (cid:90) d p (cid:48) d r (cid:48) d p (cid:104) r | e − i p ∆ t | p (cid:48) (cid:105)(cid:104) p (cid:48) | e − iV ( r )∆ t | r (cid:48) (cid:105) (cid:104) r (cid:48) | e − i p ∆ t | p (cid:105) (cid:104) p | ψ (cid:105) = 1 √ π (cid:90) d p (cid:48) d r (cid:48) d pe i rp (cid:48) e − i p (cid:48) ∆ t e − i p (cid:48) r (cid:48) e − iV ( r (cid:48) )∆ t e i r (cid:48) p e − i p ∆ t ˜ ψ ( p ) . (3)The structure of (3) suggests the following algorithmicprocedure: • Fourier transform of ψ ( r ) in order to obtain ˜ ψ ( p ) • Multiply by e − i p ∆ t • Inverse Fourier transform to real space • Multiply by e − iV ( r )∆ t • Fourier transform to momentum space • Multiply by e − i p ∆ t • Inverse Fourier transform to real spaceThe potential V consists of the harmonic potential, thescattering potential and the DDI potential. The scatter-ing potential and the DDI potential have to be calculatedat each time step. The latter can be evaluated by meansof the convolution theorem, which results in two moreFourier transforms:Φ dd ( r ) = 4 π F − (cid:26)(cid:18) k z k − (cid:19) F {| ψ ( r ) | } (cid:27) . (4)Here k and k z denote the momentum and the momen-tum in z -direction, respectively. Altogether, we have toperform six Fourier transforms for each time step. Notethat the first and last Fourier transforms described inthe algorithmic procedure of the time evolution are onlynecessary if one is interested in physical quantities whoseevaluation requires the wave function in real space.For the simulations, the spatial domain was discretizedwith up to 512 × ×
512 grid points. Since this schemeis numerically very demanding, it has been implementedfor graphics processing units (GPUs) using CUDA, en-abling a very high degree of parallelization. Using theTesla C2070 improves the performance of our algorithmby a factor of about 80 for double precision in compar-ison to the corresponding C algorithm using the wellknown FFTW library for computing the discrete Fouriertransform on a IBM System x3400 with a Quad-Core In-tel Xeon Processor E5430 (2.66GHz 12MB L2 1333MHz80w) and 4 x 4GB PC2-5300 CL5 ECC DDR2 ChipkillLow Power FBDIMM 667MHz.To investigate the coherent collision of solitons we haveapplied the following procedure. The first step is thecomputation of the ground state of one condensate usingthe split-operator method with imaginary time evolution( t = − iτ ). Afterwards we double the size of the gridin the x -direction and place two solitons in the repellingside-by-side configuration. The distance between the con-densates is chosen such that they do not feel the mutualdipole-dipole interaction. To introduce momentum in thesystem, we multiply the left hand-side of the wave func-tion by a plane wave e ikx (for the soliton moving to theright) and the right hand-side by e − ikx (for the solitonmoving to the left), respectively. III. TIME-DEPENDENT VARIATIONALANSATZ
Variational calculations using coupled Gaussian wavepackets (GWPs) have shown to be a full-fledged alterna-tive to numerical grid calculations for the calculation ofground states of dipolar BECs [9, 17]. The applicabilityof such ansatzes to dynamical simulations is a challengingtask. The decisive extension of the previous work [9, 17]is that additional translational and rotational degrees offreedom are included in the ansatz with coupled GWPs todescribe the dynamics of the condensate wave function.For the convenience of the reader we shortly review thetime-dependent variational principle (TDVP) in this sec-tion, and subsequently apply it to the ansatz of coupledGWPs. We make use of the TDVP in the formulation ofMcLachlan [18] where φ is varied such that I = || i φ − H Ψ( t ) || = min , (5)and set φ ≡ ˙Ψ afterwards. The wave function Ψ is con-sidered to be parametrized by the variational parametersΨ = Ψ( z ( t )). The minimization of the quantity I inEq. (5) leads to (cid:28) ∂ Ψ ∂ z (cid:12)(cid:12)(cid:12)(cid:12) i ˙Ψ − H Ψ (cid:29) = 0 , (6)which can be written in the short form K ˙ z = − i h , (7)with the positive definite Hermitian matrix K . We use alinear superposition of N Gaussian wave packets (GWPs)Ψ = N (cid:88) k =1 e − (cid:16) ( x T − q k ) T A k ( x − q k ) − i ( p k ) T ( x − q k ) + γ k (cid:17) ≡ N (cid:88) k =1 g k , (8)as an ansatz for the wave function in Eq. (5). In gen-eral, A k are 3 × p k and q k arethree-dimensional real vectors (representing momentumand center of the GWP) and γ k are complex numbers(where the real part stands for the amplitude and theimaginary part for the phase of the GWP, respectively).In this work we will make use of the strong confinementin one direction perpendicular to the dipole axis and omitthe translational and rotational degrees of freedom in y -direction A kyσ = A kσy = 0 , p ky = 0 , q ky = 0 , (9)with σ = x, z . Inserting the ansatz Eq. (8) in Eq. (6),sorting the result by powers of x and identifying theseterms with the coefficients of a time-dependent effectiveharmonic potential V k eff = v k + v k x + x V k x , (10)yields the equations of motion (EOM) for the variationalparameters˙ A k = − (cid:0) A k (cid:1) + i V k , (11a)˙ p k = − Re v k − A k (cid:0) ˙ q k − p k (cid:1) − V k q k , (11b)˙ q k = 2 p k + 12 (cid:0) Re A k (cid:1) − (cid:0) Im v k + 2 Im V k q k (cid:1) , (11c)˙ γ k = 2i Tr A k − i q k V k q k + 4 p k A k q k + i (cid:0) p k (cid:1) − i q k ˙ p k − i p k ˙ q k − q k A k ˙ q k + i v k . (11d)If we write Eq. (7) explicitly for GWPs, the set of linearequations for ˙ z can be rewritten to one for the vector v containing the coefficients of V k eff K v = r , (12)for details, see [19, 20]. With the transformation givenin Appendix A the EOM can now be integrated with astandard algorithm such as Runge-Kutta, where Eq. (12)has to be solved at every time step. The right-hand sidevector r with the components r l = N (cid:88) k =1 (cid:10) g l (cid:12)(cid:12) x mα x nβ V ( x ) (cid:12)(cid:12) g k (cid:11) , (13)where l = 1 , . . . , N ; α, β = 1 , . . . ,
3; 0 ≤ n + m ≤ (cid:104) Ψ | V d | Ψ (cid:105) can onlybe calculated analytically for GWPs centered in the ori-gin without the additional translational degrees of free-dom introduced in the ansatz (8). The analytical andnumerical treatment of the dipolar integral is presentedin Appendix B.The procedure for the calculations is as follows: Atfirst the equations of motion (11) for one soliton are in-tegrated in imaginary time, with the wave function be-ing normalized after every time step. Afterwards everyGWP of the wave function is copied and the resultingtwo solitons are positioned in the same way as given inSec. II. Then for each GWP a corresponding momentum p k = ± p kx e x is added. For this starting configuration theEOM are finally integrated in real time. IV. RESULTS
In Fig. 1 three grid calculations of colliding solitonswithout angular momenta prepared in the way givenabove are shown. For no phase difference constructiveinterference occurs and the condensates merge and splitup in two solitons again. Note that the condensates afterthe split-up ( t = 0 . t = 0 .
001 corresponds to 15 ms)have a larger spatial distribution than before ( t = 0 . φ = π/ φ = 0, resulting inonly slightly larger condensates at t = 0 . φ = π in phase we can seedestructive interference (column on the right), the soli-tons effectively repel each other. The transfer of kineticenergy into internal energy is once again smaller, corre-sponding to a just slightly larger condensate at t = 0 . x -directioncan be understood if one considers that a difference inphase of φ = 0 and φ = π yields a wave function whichis an eigenfunction of the parity operator, in the senseof Π ± Ψ( r , t ) = ± Ψ( − r , t ). A difference of φ = π/ t = 0 . FIG. 1. (Color online) Absorption images ( | ψ | integratedalong the y -axis) of grid calculations for the collision of soli-tons. The value of the momentum for each soliton is k = 10(velocity v = 127 µ m / s) and the field of view is 1 . × . µ m × µ m). All absorption images have been normal-ized to the maximum value. Left column: Absorption imagesfor a collision without difference in phase. Middle column:Absorption images for a collision with a difference of φ = π/ φ = π in phase. condensates with smaller extension at t = 0 .
049 thantheir corresponding condensates in the simulation pre-sented above. A difference of φ = π/ t = 0 . φ = π (Fig. 2c) shows the condensates effec-tively repelling each other, but in this case introducing (a) (b) (c) FIG. 2. (Color online) Absorption images for grid calculations and the variational ansatz of the simulation of two collidingsolitons with angular momentum. All absorption images have been normalized to the maximum value. The parameters are thesame as given in Fig. 1. The columns (a), (b) and (c) show calculations for a difference of φ = 0, φ = π/ φ = π in phase,where the left column is the result of the grid calculations and the column on the right hand-side presents the results of thevariational ansatz. For all three calculations six GWPs (three for each soliton) were used. The variational calculation is able toreproduce the transient ring-like structure during the collision for a difference of φ = π in phase and yields the correct resultsfor the configuration at the end of all three calculations. angular momentum leads to a transient ring-like struc-ture. The extension of the condensates after the collisionis much larger compared to the case with no differencein phase, which means that the amount of transferredkinetic energy in internal energy is larger than in the for-mer case. The case with angular momentum is suitedbest to show how the transfer of kinetic energy affects thespatial distribution of the condensate. In Fig. 3 we showthe kinetic energy as a function of time for the collisionswith angular momentum. Comparing the curves in Fig. 3with the absorption images in Fig. 2, it is obvious that alarger transfer of kinetic energy implies a larger conden-sate at t = 0 . | Ψ( x , t ) | increases and a small one diminishesthe coupling. This can be seen best in Fig. 1 (right col-umn) where the destructive interference for the calcula- E k i n t var. φ = 0var. φ = π /2var. φ = π num. φ = 0num. φ = π /2num. φ = π FIG. 3. (Color online) Kinetic energy as a function of timefor the collisions with angular momentum (Fig. 2). The dotsshow the numerical results, the lines show the results obtainedby the variational calculations. The kinetic energy increaseswhile the condensates merge. After the split up, the conden-sates have a lower kinetic energy than before, indicative of atransfer from kinetic to internal energy, thus resulting in exci-tation of the condensates. The inset shows the kinetic energyobtained by the variational calculations for large timescales.The oscillatory behavior indicates the excitation of the soli-tons. Note the larger kinetic energy obtained by the gridcalculations shortly after the split up. This is due to the fi-nite grid size, which manifests itself in oscillations of the wavefunction’s amplitude for large times. ∆ σ t ∆ x ∆ z ∆ x ∆ z FIG. 4. (Color online) Variance ∆ σ of a single soliton (rightsoliton in the absorption images in the upper panel) as a func-tion of time. The collision occurs at t ≈ .
03. The thick solidline and the double-dashed line show the ∆ x and ∆ z varianceof three GWPs, respectively. The thin solid line and the nor-mal dashed line show the ∆ x and ∆ z variance of the dominantGWP g . tion with phase difference φ = π leads to | Ψ( , t ) | = 0.For the corresponding calculation with nonzero angularmomentum (Fig. 2c) we find | Ψ( , t ) | = 0, too. How-ever, the ring-like structure increases the overlap duringthe collision.In Fig. 4 the variance ∆ σ = (cid:10) σ (cid:11) − (cid:104) σ (cid:105) with σ = x, z is plotted as a function of time. The variance has beencalculated for the three GWPs representing the solitonson the left-hand side in the starting configuration and forthe GWP which has the largest amplitude after the col-lision process. This dominant GWP g shows oscillatorybehavior while the other GWPs with much smaller am-plitudes describe particles leaving the soliton. This effectcan hardly be seen in the absorption images in the upperpanel of Fig. 4. However, the absorption images showthat a soliton still exists, although this would be difficultto see in an actual experiment due to the very long timescale.We have also performed simulations with smaller andlarger momenta of the solitons. The former case leadsto one merged condensate which does not split up againafter the collision but shows oscillatory behavior. Thisis very similar to the collision presented in [11]. In thelatter case the wavelength of the interference pattern issmaller and becomes more pronounced. Note that gridcalculations with high momenta are problematic, becausethe condensates quickly reach the edge of the grid. Anapproach with a variational ansatz is better suited toanalyze these scenarios. V. CONCLUSION
We have studied the collisions of anisotropic two-dimensional bright solitons in dipolar Bose-Einstein con-densates both with a fully-numerical ansatz and a time-dependent variational principle with coupled Gaussians.The calculations presented show that the collision pro-cess leads to an energy transfer from kinetic energy to“inner” energy of the solitons which leads to excited soli-tons with larger extent. The absorption images show verygood qualitative agreement of the results gained by thetwo different methods.The advantages of the grid calculations are the sim-plicity of the numerical scheme (although the implemen-tation for the massively parallel computation requiressome effort), the freedom in describing all different shapesof wave functions, and the numerical stability of themethod. The advantages of the variational calculationsare the much smaller numerical effort, enabling one to runlong calculations on standard PCs, the independence offinite grid size, and the small amount of parameters tobe saved.Both methods can be used to simulate the time-dependent GPE, supporting each other mutually. Onefurther application would be the inclusion of additionalexternal potentials such as optical lattices and the com-parison of the methods in such scenarios. Our resultsshould stimulate experimental efforts to study the colli-sions of 2D anisotropic solitons.
VI. ACKNOWLEDGEMENTS
We thank Boris Malomed for valuable discussions.This work was supported by Deutsche Forschungsgemein-schaft. R.E. is grateful for support from the Landes-graduiertenf¨orderung of the Land Baden-W¨urttemberg.
Appendix A: Transformation to CB -variables The direct numerical integration of Eq. (11) leads tonumerical difficulties [21]. These can be dealt with by theintroduction of two auxiliary matrices B and C . With A = BC − the equations of motion for the width matri-ces can be written as˙ A k = − (cid:0) A k (cid:1) + i V k , (A1)˙ A k = ˙ B k (cid:0) C k (cid:1) − − B k (cid:0) C k (cid:1) − ˙ C k , (A2)where C and B are 3 × k we obtain from these equations B − ˙ BC − − C − ˙ C = − B − A + i B − V (A3) ⇒ B − ˙ B − C − ˙ CC = − C − B + i B − V C . (A4)By comparison we yield the equations of motion for C and B ˙ B k = i V k C k , (A5)˙ C k = 4i B k . (A6)The reduction (9) can be done for those matrices, too.Note however, that the matrices B and C do not preservethe same symmetry as the matrices A which are complexsymmetric. Therefore, all five complex entries in B and C have to be integrated. Appendix B: Solution of the dipolar integral
The calculation of the dipolar integrals needed in theTDVP (cid:104) Ψ | α n β m V d | Ψ (cid:105) with α, β = x, y, z and 0 ≤ n + m ≤ n = m = 0.The other integrals are calculated analogously. We start from the six-dimensional non-local integral (cid:104) Ψ | V d | Ψ (cid:105) = (cid:88) l,k,j,i (cid:90) (cid:90) d r d r (cid:48) g l ∗ ( r ) g j ∗ ( r (cid:48) ) g i ( r (cid:48) ) g k ( r ) × (cid:32) − z − z (cid:48) ) | r − r (cid:48) | (cid:33) | r − r (cid:48) | . (B1)By the use of the convolution theorem of Fourier analysiswe can evaluate one of the three-dimensional integralsdirectly, while the inverse Fourier transform (cid:104) Ψ | V d | Ψ (cid:105) = 16 π (cid:88) l,k,j,i I kl I ij × (cid:90) d k exp (cid:26) − k T ¯ A klij k + 12 i (cid:0) ¯ p klij (cid:1) T k (cid:27) × (cid:18) k z k − (cid:19) , (B2)remains to be done. Here I kl denotes the overlap integralof the Gaussian functions k and l , and we have used theabbreviations¯ A klij = ( A kl ) − + ( A ij ) − , (B3)¯ p klij = ( A ij ) − p ij − ( A kl ) − p kl , (B4)with A kl = A k + A l ∗ and p kl = p k + p l ∗ and analo-gously for i and j . The integral (B2) can be split intwo parts, one leading to a shift in the scattering length(this is the short-range part of the DDI) and a secondpart (cid:104) Ψ | V d , eff | Ψ (cid:105) = (cid:80) l,k,j,i I kl I ij J klij . After a principalcomponent analysis of the exponential in Eq. (B2) theanalytical integration in k y -direction is possible when wemake use of Eq. (9). The remaining result reads J klij = 14 π ∞ (cid:90) d ρ w (cid:18) i (cid:113) ¯ A klijy ρ (cid:19) ρ e − ( ¯ A klijx + ¯ A klijz ) ρ × (cid:88) ± x c ± x s (cid:90) − d x ( ± c x c ± c x s ) √ − x × e − ( ¯ A klijx − ¯ A klijz ) ρ x + i2 ( ± ¯ p klijx ρx c ± ¯ p klijz ρx s ) , (B5)with x c = (cid:112) (1 + x ) / , x s = (cid:112) (1 − x ) /
2, the coeffi-cients c , c of the rotation matrix from the principalcomponent analysis and the Faddeeva function w ( z ) =e − z erfc( − i z ). The numerical evaluation of this integralcan efficiently be performed by a Taylor expansion of theFaddeeva function for which the single terms can be ob-tained by a recursion formula and using a Chebyshevquadrature for the x -integration. To improve the resultwe apply a Pad´e-approximation to the Taylor series.The numerical integration of the dipolar integrals is thecrucial part in this method. Dependent on the number ofGaussian functions N there is a total number of C num =( N + N − N ) / C elliptic = N ( N + 1) / [1] A. Griesmaier, J. Werner, S. Hensler, J. Stuhler,and T. Pfau, Phys. Rev. Lett. , 160401 (2005);Q. Beaufils, R. Chicireanu, T. Zanon, B. Laburthe-Tolra,E. Mar´echal, L. Vernac, J.-C. Keller, and O. Gorceix,Phys. Rev. A , 061601(R) (2008).[2] M. Lu, S. H. Youn, and B. L. Lev, Phys. Rev. Lett. ,063001 (2010); M. Lu, N. Q. Burdick, S. H. Youn, andB. L. Lev, Phys. Rev. Lett. , 190401 (2011).[3] K. Aikawa, A. Frisch, M. Mark, S. Baier, A. Rietzler,R. Grimm, and F. Ferlaino, Phys. Rev. Lett. , 210401(2012).[4] K.-K. Ni, S. Ospelkaus, M. H. G. de Miranda, A. Pe’er,B. Neyenhuis, J. J. Zirbel, S. Kotochigova, P. S. Julienne,D. S. Jin, and J. Ye, Science , 231 (2008).[5] T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, andT. Pfau, Rep. Progr. Phys. , 126401 (2009).[6] L. Khaykovich, F. Schreck, G. Ferrari, T. Bourdel, J. Cu-bizolles, L. D. Carr, and Y. Castin, Science , 1290(2002); K. E. Strecker, G. B. Partridge, A. G. Truscott,and R. G. Hulet, Nature , 150 (2002).[7] I. Tikhonenkov, B. A. Malomed, and A. Vardi, Phys.Rev. Lett. , 090406 (2008).[8] P. K¨oberle, D. Zajec, G. Wunner, and B. A. Malomed,Physical Review A , 023630 (2012). [9] R. Eichler, J. Main, and G. Wunner, Phys. Rev. A ,053604 (2011).[10] S. K. Adhikari and P. Muruganandam, J. Phys. B: At.Mol. Opt. Phys. , 045301 (2012).[11] P. Pedri and L. Santos, Phys. Rev. Lett. , 200404(2005).[12] S. Giovanazzi, A. G¨orlitz, and T. Pfau, Phys. Rev. Lett. , 130401 (2002).[13] L. E. Young, P. Muruganandam, and S. K. Adhikari, J.Phys. B: At. Mol. Opt. Phys. , 101001 (2011).[14] P. K¨oberle, H. Cartarius, T. Fabˇciˇc, J. Main, andG. Wunner, New Journal of Physics , 023017 (2009).[15] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev.Mod. Phys. , 1225 (2010).[16] M. D. Feit, J. A. Fleck, Jr., and A. Steiger, J. Comp.Phys. , 412 (1982).[17] S. Rau, J. Main, P. K¨oberle, and G. Wunner, Phys. Rev.A , 031605(R) (2010).[18] A. D. McLachlan, Mol. Phys. , 39 (1964).[19] S. Rau, J. Main, and G. Wunner, Phys. Rev. A ,023610 (2010).[20] S. Rau, J. Main, H. Cartarius, P. K¨oberle, and G. Wun-ner, Phys. Rev. A , 023611 (2010).[21] E. J. Heller, J. Chem. Phys.65