From Feshbach-Resonance Managed Bose-Einstein Condensates to Anisotropic Universes: Applications of the Ermakov-Pinney equation with Time-Dependent Nonlinearity G. Herring1, P.G. Kevrekidis1, F. Williams1, T. Christodoulakis2 and D.J. Frantzeskakis2 1 Department of Mathematics and Statistics, University of Massachusetts, Amherst MA 01003-4515, USA 2 Department of Physics, University of Athens, Panepistimiopolis, Zografos, Athens 15784, Greece 7 0 Inthisworkwerevisitthetopicoftwo-dimensionalBose-Einsteincondensatesundertheinfluence 0 oftime-dependentmagneticconfinementandtime-dependentscatteringlength. Amomentapproach 2 reducestheexaminationofmomentsofthewavefunction(inparticular,ofitswidth)toanErmakov- n Pinney(EP)ordinarydifferentialequation(ODE).Weusethewell-knownstructureofthesolutions a of this nonlinear ODE to “engineer” trapping and interatomic interaction conditions that lead to J condensates dispersing, breathing or even collapsing. The advantage of the approach is that it is 0 fullytractableanalytically,inexcellentagreementwithournumericalobservations. Asanaside,we 3 also discuss how similar time-dependent EP equations may arise in the description of anisotropic scalar field cosmologies. ] r PACSnumbers: e h t o I. INTRODUCTION . t a m The recent progress in experimental and theoretical studies of Bose-Einstein condensates (BECs) of dilute atomic - gases [1] has been tremendous after their experimental realization [2]. This has also led to an explosion of interest in d the theme ofnonlinearmatter-wavessuchasdark[3], bright[4]andgap[5]solitons. Two-dimensional(2D)nonlinear n excitations of BECs, such as vortices [6] and vortex lattices [7], were also realized experimentally, while a multitude o of other coherent nonlinear structures were also theoretically predicted; these include, e.g., Faraday waves [8], ring c [ dark solitons and vortex necklaces [9], stable solitons and localized vortices in attractive BECs trapped in periodic potentials [10], matter-wave gap vortices [11], 2D bright solitons in dipolar BECs [12], and so on. 1 From the theoretical standpoint, the dynamics of such higher-dimensional structures is, generally, difficult to be v treatedanalytically,especiallysointime-dependentsettings. This,inturn,raisesimportantquestionsconcerningtheir 6 5 “controllability”,whichwouldbeofparticularrelevanceregardingpotentialapplications. Ideally,suchacontrollability 7 would allow “manipulation” of the condensates at will, e.g., sustaining condensates which may expand towards 1 complete dispersion, contract towards a complete collapse, or perform stable breathing oscillations. Moreover, such 0 processes could involve a targeted growth of a condensate up to a certain width or shrinkage down a desired size. 7 Our purpose in the present work is to illustrate how one can achieve this goal, by taking advantage of one of the few 0 analytically tractable tools in higher-dimensional condensates, namely of the analysis of their moments [13, 14, 15] / t (see also references therein). In particular, we focus on the case of quasi two-dimensional (2D), so-called, pancake- a m shaped condensates [16], under the effect of time-varying harmonic trapping [14, 15] and also time-dependent s-wave scattering length (which controls the inter-particle interactions). Note that the controllability over the harmonic - d trapping is straightforwardly realizable under AC-variation of the atom trapping frequencies, while the controllabity n over the interatomic interactions can be realized by using the so-called Feshbach resonance [17], connecting the s- o wave scattering length to external magnetic fields. This latter technique (which is usually called Feshbach resonance c management (FRM)) has been proposed theoretically as a means of avoiding collapse [18, 19, 20], but also as a way : v of producing robust coherent nonlinear matter-waves [21, 22, 23, 24]. i Forthequasi-2Dcondensatesdiscussedabove,uponpresentingthemomentanalysis,weobtainareduceddynamical X description effectively involving only a variable associated with the width of the condensate wavefunction. This r a ordinary differential equation (ODE) can be reduced to one of the Ermakov-Pinney(EP) type [25, 26, 27, 28], whose solutions can be obtained analytically, provided that the solutions of the underlying linear Schro¨dinger equation can be obtained. We use this feature and the freedom in selecting the time-dependent trapping and interactions of the condensates to illustrate that one can construct analytical solutions to this ODE that completely characterize the temporal evolution of the width of the wavefunction. In so doing, we fully prescribe the dynamical evolution of the condensate. We show three prototypical examples in applying this idea: one in which the width grows in time (leading to condensate expansion), one in which it decreases (leading to focusing), and one in which it periodically “breathes”betweenaminimumandamaximumvalue. Inallthreecases,wetesttheanalyticalpredictionagainstthe full numerical simulationof the mean-field partial differential equation (PDE) model fully describing the condensate. We find excellent agreement between the two, showcasing the accuracy of our theoretical approach. The analysis of the ensuing EP equation in this setting is of interest in its own right. This is because, as we will 2 see, the resulting EP equation has a time-dependent nonlinearity in the right hand side (contrary to what is the case forthe “standard”EPframework;seee.g. [25,26,27,28]andreferencestherein). Inthiscontext,andasanaside,we present a secondphysically relevantexample where such EP equations with time-dependent nonlinearities may arise, by studying anisotropic scalar field cosmologies of a particular anisotropic geometry. This generalizes the example of [25] parallelizing BECs without FRM and scalar field cosmologies in the isotropic case. Ourpresentationwillbe structuredasfollows: InsectionIII,weprovideabriefsynopsisofthe mainfeaturesofthe EPequation. InsectionIV,weproceedtoanoverviewofthemomentanalysis,followingtheearlierworks[13,14,15]. InsectionIV, wepresentthe analyticalsolutionsofthe EPequationthatwe developforeachofthe above-mentioned three cases. Section V tests these analytical results against full simulations of the PDE describing the condensate. SectionVIillustratesthesecondexampleoftime-dependentEPequationsinthedescriptionofanisotropicscalarfield cosmologies. Finally, in section VII, we summarize our findings and present our conclusions. II. THE ERMAKOV-PINNEY EQUATION The Ermakov-Pinney(EP) equation is a remarkable nonlinear ODE of the form: κ ′′ Y +Q(τ)Y = . (1) Y3 TheparticularlyattractivefeatureofthisnonlinearODEisthatitsgeneralsolutioncanbeobtained,providedthatone ′′ is able to solvethe time-independent linear Schr¨odinger(LS) equationY +Q(τ)Y =0. For details onthe properties of the EP equation, the interested reader is referred to [26, 27] and references therein. Here we just mention its basic superpositionprincipleproperty. Namely,ifthelinearlyindependentsolutionsoftheLSequationareY (τ)andY (τ), 1 2 then the most general possible solution of the EP equation is given by Y(τ)= AY2+BY2+2CY Y 1/2 (2) 1 2 1 2 (cid:0) (cid:1) where A,B and C are constants connected through κ AB−C2 = (3) W2 ′ ′ where W =Y Y −Y Y is the Wronskian of Y (τ) and Y (τ). 1 2 2 1 1 2 III. MOMENT ANALYSIS FOR BECS One of the interesting variants of the “regular” EP equation of the form (1) arises in the study of BECs, albeit in a somewhat modified form (see below). The relevant mean-field model for studying atomic Bose-Einstein condensates at zero temperature consists of the so-called Gross-Pitaevskiiequation [1] of the following dimensionless form: 1 i∂ u=− ∆u+ λ(t)r2+ν(t)|u|2 u, (4) t 2 (cid:0) (cid:1) whereurepresentsthewavefunctionofthecondensate,V(r)=λ(t)r2 denotestheharmonictrapconfiningthebosons [with a time-dependent frequency determined by λ(t)] and ν(t) is the coefficient of the nonlinear term proportional to the s-wave scattering length, characterizing the interaction between the particles. We will take advantage of the Feshbach resonance [17] to consider that the latter is time-dependent, as well. One of the popular approaches to studying Eq. (4) is through the use of moment methods [13, 14, 15]. The latter allow us to write ODEs for the moments of the mean-field wavefunction u as follows. We define ∞ I(d) = ra|u|2rd−1dr, (5) 2,a Z 0 ∞ I(d) = i ra(uu⋆−u⋆u )rd−1dr, (6) 3,a Z r r 0 ∞ 2 ∂u I(d) = ra rd−1dr, (7) 4,a Z (cid:12)∂r(cid:12) 0 (cid:12) (cid:12) ∞ (cid:12) (cid:12) I(d) = ra|(cid:12)u|4r(cid:12)d−1dr, (8) 5,a Z 0 3 wheresubscriptsdenotepartialdifferentiation,“⋆denotescomplexconjugateand(d)indexesthedimension(similarly to [13, 15]). The Hamiltonian of Eq. (4) ∞ 1 H = |∇u|2+ν(t)|u|4+2λ(t)r2|u|2 rd−1dr, (9) 2Z 0 (cid:2) (cid:3) can then be written as: 1 1 H = I(d)+ ν(t)I(d)+λ(t)I(d). (10) 2 4,0 2 5,0 2,2 One can then infer from the dynamics of Eq. (4) that: ′ ν (t) H˙ = I(d)+λ′(t)I(d) (11) 2 5,0 2,2 However, as derived in [13, 15], I¨(d) =4H −8λ(t)I(d)+2ν(t)(d−2)I(d). (12) 2,2 2,2 5,0 We focus on I(d), since this moment is associatedwith the width of the spatial profile of the wavefunction. From the 2,2 latter, we can infer a number of useful pieces of information concerning certain asymptotic values; in particular, if I(d) → 0, the condensate collapses, while if I(d) → ∞, the BEC disperses. Furthermore, since the L2 norm of the 2,2 2,2 wavefunction is conserved, the estimate of I(d) on the (square) width of the wavefunction can be used together with 2,2 thisconservationlawtoprovideinformationontheamplitudeAofthewavefunction(i.e.,approximatelyA∼1/I(d)). 2,2 Combining the two equations (11) and (12), we obtain a single equation for the time-dependence of I(d) as 2,2 ·I··(d) =2ν′(t)I(d)−4λ′(t)I(d)−8λ(t)I˙(d)+(d−2) 2ν′(t)I(d)+2ν(t)I˙(d) (13) 2,2 5,0 2,2 2,2 5,0 5,0 h i Following [20], we consider a quadratic phase for the solution (an assumption most relevant to dimension d=2 as discussed in [20], but which we will also consider more generally), and obtain (cf. Eq. (6d) of [20]) that K I(d) = , (14) 5,0 (d) I 2,2 where the constant K is determined by initial conditions. Notice that this is the only assumption in our calculations herein,whosevaliditywillbeexaminedaposterioribycomparingouranalyticalresultswithnumericalcomputations. Denoting for simplicity I(d) =y, the resulting ordinary differential equation for y can be written as: 2,2 d 1 d ν(t) yy¨− y˙2−2Kν(t)+4λ(t)y2 =2K(d−2)y (15) dt(cid:18) 2 (cid:19) dt(cid:18) y (cid:19) Clearly, from the above exposition, the most straightforward case is the one with d = 2, on which we will focus next. When d=2, the equation can be directly integrated, yielding 1 2Kν(t)+C y¨− y˙2+4λ(t)y = , (16) 2y y where C is an integration constant that can be computed from the initial conditions as: 1 C =4I(2)(0) H(0)−λ(0)I(2)(0) − (I(2)(0))2−2Kν(0) (17) 2,2 2,2 2 3,1 (cid:16) (cid:17) If we now use the transformationY(t)= y(t) [14, 25], then we obtainan Ermakov-Pinney(EP) type equationof the form: p Kν(t)+ C Y¨ +2λ(t)Y = 2 (18) Y3 Hence, Eq. (18) is the equation that describes the dynamics of two-dimensional Bose-Einstein condensates in the presence of a time-dependent trap [13, 14, 25], as well as in the scenario of Feshbach resonance management [18, 19, 20, 21, 22, 23]. One can now examine particular cases of time-dependence of λ(t) and ν(t), using the analytical tractability of the ensuing EP equation, in order to obtain completely analytical solutions for Y(t), and hence for I(2). 2,2 4 IV. ANALYTICAL RESULTS There are numerous possiblities in the case of a time-dependent λ and ν, thus, as explained above, we limited our investigation to the three cases of an expanding waveform, a collapsing waveform and an oscillatory waveform. In order to derive examples for each of these cases, ν(t) was chosen to be independent of the choice of λ(t) and y(t). Solutions to Eqn. (18) were determined by “reverseengineering”[28] (based on the desired behaviorof the dynamics of the wavefunction), and the details of the EP functional form. Below are the three cases examined. Notice that in all three cases, the function f(t) used below is given by f(t)=Kν(t)+C/2 i.e., by the numerator of the right hand side of Eq. (18). 1. The first case explores the possibility of an expanding width wavefunction. In orderfor the wavefunctionwidth to increase, λ(t) must be decreasing, thus test functions for λ(t) and y(t) were chosen to reflect this. A similar methodology was also used in the two other cases to determine test functions. y(t) = B2 A2+t2 (19) f((cid:0)t)−A2B(cid:1)4 λ(t) = (20) 2B4(A2+t2)2 z(t) = B4A2 =const. (21) 2. The second case, a focusing wavefunction, required an additional condition (y(0) 6= ∞) in order to perform the numerical simulations described below. Additionally, in this case and the next, an initially undetermined functionz(t)wasusedinλ(t) tosimplify the derivationofλ(t)andy(t). After thesewerefound, the exactform of z(t) was calculated by plugging λ(t) and y(t) into Eqn. (18) and solving for z(t). This results in B2 y(t) = (22) A2+tp Bq λ(t) = [f(t)−z(t)] (23) [2(A2+tp)]q/2 B5 B4ptp−2(ptp−2pA2+2A2+2tp f(t)(A2+tp)5/2 z(t) = + −f(t) (24) (A2+tp)(cid:20) 4(A2+tp)4 B5 (cid:21) 3. Finally, in the last case, equations for λ(t) and y(t) were calculated in order to generate a wavefunction with oscillatory width. The corresponding functions in this case read y(t) = b+sin2(ct) (25) f(t)−z(t) λ(t) = (26) 2 b+sin2(ct) 2 (cid:2)c2 (cid:3) z(t) = − sin2(2ct)+c2cos(2ct) b+sin2(ct) (27) 4 (cid:2) (cid:3) The resulting form of the normalized confining frequency λ(t) for the three different cases is shown in Fig. 1. We note in passing that such schemes are definitely realizable within the unprecedented control that exists over the magnetic confinement of the condensates. A very recent example illustrating this point can be found in the recent work of [29]; this experiment clearly realizes a form of the third type of confinement (an oscillatory one), by using a periodic modulation of the transverse confinement of the condensate, in order to induce longitudinal oscillations and the emergence of Faraday patterns. This directly shows that temporal modulation of the trapping frequencies is feasible. Whatweproposehereconsistsofthethreeprincipaltypesoftrappingfrequencydependence,i.e.,monotonic decrease, monotonic increase and oscillatory dependence, that allow an explicit analytical handle on the dynamics of the BEC, in excellent agreementwith our numerical findings (see below). Hence, we expect that such schemes would be directly applicable in experimental settings similar to those of the above experiment. V. NUMERICAL SETUP AND RESULTS In order to validate our analytical results based on the predictions made above, the evolution of the condensates under the time-dependent trappings imposed by Eqs. (20), (23) and (26) were tested by numerically solving the 5 2.5 500 0.5 2 400 1.5 300 λ λ λ 0 1 200 0.5 100 0 0 −0.5 0 1 2 3 4 5 0 1 2 3 4 5 0 2 4 6 8 10 Time(t) Time(t) Time(t) FIG.1: Plots of λ(t)forCases 1-3, i.e., respectively expanding(left), contracting(middle) andoscillatory (right)condensates. original PDE model of Eq. (4). The numerical simulations were run using spectral methods [30], with a spatial grid composed of Chebychev nodes (x = cos(πn/N) where N is the total number of nodes). Since the Chebychev n nodes are contained within the interval [−1,1], the spatial variable was normalized in order for the solution u(t) to be supported on the interval [−1,1] from t = 0 to some evolution horizon t = T > 0. The temporal integration was implemented using a 4th-order Runga-Kutta scheme. It is worth noting here that in all three cases examined, the choice of ν(t) was independent of the equations for λ(t) and y(t), thus ν(t) = sin2(t) was chosen to be used for all cases, in consonance with the Feshbach resonance management scheme discussed in the Introduction; in that sense, however, notice that the key temporal variation in the schemes presentedherein is that of the magnetic trapping frequency. As an initial conditionto Eqn. (4), we used a generic Gaussian profile of the form u(t) = Aexp(−Br2/2) with the values of A and B subject to the condition y(0)=A2/(2B2). For the cases of the expanding and collapsing width wavefunctions, the numerical results are in perfect agreement with the analyticalcalculations as shownin Figs. (2) and (3). Notice in the lower left-hand plot for both figures how the lines for the analytical and numerical plots are completely indistinguishable from each other. For the last case of the oscillatory width wavefunction, there was a very small discrepancy between the analytical calculations and numerical simulation, as seen in the lower left-hand plot of Fig. (4). This discrepancy is only seen at the maximum values of the wavefunction width, starting with the second and subsequent maximums. While this small discrepancy may be triggered by the sole approximation of our approach, namely Eq. (14), further numerical experiments (not shown here) seem to indicate that it is more likely to be the result of the numerical approximation to the corresponding moment. In any case, the overall excellent numerical agreement between the analytically obtained moment I(2) and its 2,2 numericallyfoundcounterpartillustratetherelevanceofthisapproachanditsusefulnessinsystematicallyprescribing the wavefunction behavior, in an analytically tractable way. VI. ANISOTROPIC SCALAR FIELD COSMOLOGIES An interesting analogy has been recently made between Bose-Einstein condensates and isotropic scalar field cos- mologies in the absence of time dependence in the coefficient of the nonlinearity [25]. The so-called Friedmann- Robertson-Walker (FRW) metric, coupled to a scalar field had been shown earlier [26, 28] to be described by an EP equation and that viewpoint was exploited in [25] to illustrate an analogy between cosmic dynamics and BECs. Here, we extend this analogy to EP equations with time-dependent nonlinearities. In particular, we show that anisotropic cosmic dynamics can be parallelized with Feshbach-managedBose condensates. If we consider the spatially homogeneous, yet anisotropic geometry of Bianchi type I, we have a line element of: ds2 =−N(t)2dt2+A(t)2dx2+B(t)2dy2+Γ(t)2dz2. (28) We can define the tensor: F = G − 8πT where G = R − 1g R is the Einstein tensor and T = µν µν µν µν µν 2 µν µν φ φ −1g (φ;αφ +m2φ2)istheenergymomentumtensor. Then,thequadraticconstraintistheequationF0 =0, ;µ ;ν 2 µν ;α 0 the kinematic equation is given by F1 = 0, while the Klein-Gordon equation for the scalar field (that is coupled to 1 gravity) is given by φ;µ −m2φ = 0 ∝ Tµν = 0. Notice additionally, that the two integrals of the motion, namely ;µ ;ν I =F1−F2 =0 and I =F1−F3 =0 yield B(t)=A(t)eκt/2 and Γ(t)=A(t)eλt/2. 1 1 2 2 1 3 6 t = 0 t = 1 −5 −5 Y 0 Y 0 5 5 −5 0 5 −5 0 5 X X 1.5 t = 2 t = 3 Analytical −5 −5 Numerical 1 h dt Wi Y 0 Y 0 0.5 0 5 5 0 1 2 3 4 5 −5 0 5 −5 0 5 Time(t) X X FIG. 2: Evolution plots of the expanding width wavefunction of Eqs. (19)-(20). The upper left figure shows the spatio- temporal evolution of a fixed contour of the solution over time. The lower left figure compares the analytically determined moment associated with the wavefunction width with its numerically computed counterpart. The remaining figures show the contourofthesolutionatvarioustimesduringtheevolution,clearlyindicatingtheexpandingnatureoftheBECwavefunction. SolvingtheKlein-Gordonequationforφ′′(t)andsubstitutingtheresultinto∂ F0 =0(aswellassolvingF0 =0for t 0 0 φ and using the resulting expression in ∂ F0 = 0), one is led to a dynamical equation for the remaining scale factor t 0 A(t) in the form: κλ κA′(t) λA′(t) 4A′(t)2 φ′(t)2 A′′(t) + + + − − =0 (29) 4 A(t) A(t) A(t)2 2 A(t) Using now: A(t)=Y(t)2/n and a change of variable τ = tΩ(t′)dt′, we obtain: R Ω′(t) Y˙(τ)2(6+n) Y˙(τ) nκλ Y(τ) nY(τ)φ˙(τ)2 Y¨(τ)+Y˙(τ) − −(κ+λ) − + =0 (30) Ω(t)2 Y(τ) n Ω(t) 8 Ω(t)2 4 Hence, a choice of time reparametrization according to: ′ ′ Ω(t) (6+n)Y (t) =κ+λ+ (31) Ω(t) n Y(t) (which leads to Ω(t)=θe(κ+λ)tY(t)(6+n)/n, where θ >0 is a constant of integration), results in the form: Γ(t(τ)) φ˙(τ)2 nκλe−2(κ+λ)t Y¨(τ)+Q(τ)Y(τ)= , Q=n , Γ= . (32) Y(τ)1+12/n 4 8θ2 Eq. (32)becomesanEPequationforthechoiceofn=6. NoticethatinEq. (32),anontrivialcomplicationisthatthe time dependentcoefficientΓdepends onthe reparametrizationoftime throughΩ. Nevertheless,theoriginalequation (29) has been solved in [31]. This may, in turn, provide valuable insights in the solution of such time-dependent EP equations; this would be an interesting direction for future studies in such time-dependent EP settings. We close this section by noting that one can therefore generalize the analogy of Bose condensates and scalar field cosmologiesintheFeshbach-managedviz. anisotropiccaseoftheEPequationwithtime-dependentnonlinearity. The 7 t = 0 t = 1 −5 −5 Y 0 Y 0 5 5 −5 0 5 −5 0 5 X X 0.5 t = 2 t = 3 Analytical −5 −5 Numerical 0.4 0.3 h dt Wi Y 0 Y 0 0.2 0.1 0 5 5 0 1 2 3 4 5 −5 0 5 −5 0 5 Time(t) X X FIG. 3: Evolution plots of the focusing case of Eqs. (22)-(24), with p=2 and q=−9. The upperleft figureshows thespatio- temporal evolution of a fixed contour of the solution over time. The lower left figure compares the analytically determined wavefunction width with its numerically computed analog. The remaining figures show the contour of the solution at various times, illustrating thefocusing effect. quantities that now bear the analogy is the scale factor of Bianchi type I with the second moment of the condensate wavefunction, the time-dependent magnetic trap strength with the time-dependent scalar field and finally the time dependent nonlinearity coefficient (i.e., the scattering length in the condensate dynamics) is (indirectly) connected with the reparametrizationof time in the cosmologicalproblem. VII. CONCLUSIONS AND FUTURE CHALLENGES In this short communication, we revisited the theme of higher dimensional Bose-Einstein condensates under the presence of magnetic trapping and Feshbach Resonance management. We used the moment method to develop an Ermakov-PinneyODE for the second moment of the distribution function, which is associated with the width of the condensate. ThisEPequationisanalyticallytractable,inanumberofcases. Inparticular,onecan“reverseengineer” magnetic trappings that will induce the expansion, contraction or oscillatory behavior of the condensate at will. For such scenaria, the relevant moment of the wavefunction can be obtained analytically and is found to be in excellent agreementwithournumericalsimulationsofthefulloriginalPDEmodel. The approachpermitsadetailedanalytical handle on the behavior of higher dimensional condensates which is usually quite difficult to acquire with different methods based on nonlinear PDEs. It may also, in turn, permit to appropriately craft experiments, based on the reshaping “operation” that is desirable to perform on the condensate. As an aside example, we have illustrated the relevance of such time-dependent Ermakov-Pinney equations in a completely different physical system, namely in anisotropic scalar field cosmologies of Bianchi type I. This approach also suggests a number of interesting questions. It would be, in particular, relevant to examine whether a general solution of the time-dependent EP equation developed herein can be obtained on the basis of its linear Schr¨odinger counterpart. This would bear direct consequences both in the atomic physics and in the cosmologicalproblem, allowing for a general analytical description of their respective time-dependent properties (for the condensate wavefunction width or the cosmological scale factor, respectively). Furthermore, the approach was presentedherefor two-dimensionalsettings forreasonsthathaveto do withthe simplification/closureof the moment 8 t = 0 t = 1.6 −5 −5 Y 0 Y 0 5 5 −5 0 5 −5 0 5 X X t = 3.1 t = 4.7 1.8 −5 −5 Analytical 1.6 Numerical 1.4 h1.2 dt Y 0 Y 0 Wi 1 0.8 0.6 0.4 5 5 0 2 4 6 8 10 −5 0 5 −5 0 5 Time(t) X X FIG.4: EvolutionplotsfortheoscillatorywidthwavefunctioncaseofEqs. (25)-(27). Thepanelsaresimilartotheprevioustwo figures. Themiddleandrightpanelsshowthecontourofthesolution atthefirsttwomaximaandminimaofthewavefunction width. approachinthatcase. However,itwouldbe ofparticularinterestto developsimilarapproachesandpotentialclosure schemesforone-dimensionalorthree-dimensionalsettings. Suchstudiesarecurrentlyinprogressandwillbereported in future publications. [1] F. Dalfovo et al., Rev.Mod. Phys. 71, 463 (1999). [2] M.H. Anderson et al, Science 269, 198 (1995); K.B. Davis et al, Phys. Rev. Lett. 75, 3969 (1995); C.C. Bradley et al, Phys. Rev. Lett. 75, 1687 (1995); see also the reviews E. A. Cornell and C. E. Weiman, Rev. Mod. Phys. 74, 875 (2002) and W. Ketterle, Rev.Mod. Phys. 74, 1131 (2002). [3] S.Burger et al.,Phys.Rev.Lett.83, 5198 (1999); J.Denschlaget al.,Science287, 97(2000); B.P.Andersonet al.,Phys. Rev.Lett. 86, 2926 (2001). [4] K.E. Strecker et al.,Nature 417, 150 (2002); L. Khaykovichet al., Science296, 1290 (2002). [5] B. Eiermann et al., Phys.Rev.Lett. 92, 230401 (2004). [6] M.R. Matthews et al., Phys.Rev.Lett. 83, 2498 (1999); K.W. Madison, et al. Phys. Rev.Lett.84, 806 (2000); S.Inouye et al.,Phys. Rev.Lett. 87, 080402 (2001). [7] J.R.Abo-Shaeeretal.,Science292,476(2001);J.R.Abo-Shaeer,C.RamanandW.Ketterle,Phys.Rev.Lett.88,070409 (2002); P. Engels et al.,Phys.Rev. Lett.89, 100403 (2002). [8] K.Staliunas, S.Longhi and G.J. deValc´arcel, Phys.Rev.Lett. 89, 210406 (2002). [9] G. Theocharis et al.,Phys. Rev.Lett. 90, 120403 (2003). [10] B.B. Baizakov, B.A. Malomed, and M. Salerno, Europhys. Lett. 63, 642 (2003); D. Mihalache et al., Phys. Rev. A 72, 021601 (2005). [11] E.A. Ostrovskayaand Yu.S.Kivshar, Phys.Rev.Lett. 93, 160405 (2004). [12] P.Pedri and L. Santos, Phys.Rev.Lett. 95, 200404 (2005). [13] J.J. Garcia-Ripoll and V.M. P´erez-Garc´ıa, patt-sol/9904006; see also V.M. P´erez-Garc´ıa et al., J. Opt.B 2, 353 (2000). [14] J.J. Garcia-Ripoll, V.M. P´erez-Garc´ıa, and P. Torres, Phys.Rev.Lett. 83, 1715 (1999). [15] P.G. Kevrekidis, A.R.Bishop and K.Ø.Rasmussen, J. Low Temp. Phys. 120, 205-212 (2000). [16] A. G¨orlitz et al., Phys. Rev. Lett. 87, 130402 (2001); M. Greiner et al., Phys. Rev. Lett. 87, 160405 (2001); S. Burger et al., Europhys. Lett.57, 1 (2002). 9 [17] S. Inouye et al., Nature 392, 151 (1998); J. L. Roberts et al., Phys. Rev. Lett. 81, 5109 (1998); J. Stenger et al., Phys. Rev.Lett. 82, 2422 (1999); S. L. Cornish et al.,Phys.Rev. Lett.85, 1795 (2000). [18] F.Kh.Abdullaev et al.,Phys. Rev.A 67, 013605 (2003); H. Saito and M. Ueda, Phys.Rev.Lett. 90, 040403 (2003). [19] F. Kh.Abdullaev et al., Phys.Rev.A 68, 053606 (2003). [20] G.D. Montesinos, V.M. P´erez-Garc´ıa, and P.J. Torres, Physica D 191, 93 (2004). [21] P.G. Kevrekidis et al., Phys. Rev.Lett. 90, 230401 (2003). [22] D.E. Pelinovsky,P.G. Kevrekidisand D.J. Frantzeskakis,Phys. Rev.Lett. 91, 240201 (2003). [23] Z. Raptiet al.,Phys. Scripta T107, 27-31 (2004). [24] F.Kh.Abdullaev et al.,Phys. Rev.Lett.90, 230402 (2003). [25] J.E. Lidsey,Class. QuantumGrav. 21, 777 (2004). [26] R.M. Hawkins and J.E. Lidsey,Phys. Rev.D 66, 023523 (2002). [27] P.B. Espinoza Padilla, math-ph/0002005. [28] F.L. Williams and P.G. Kevrekidis, Class. QuantumGrav. 20, L177 (2003). [29] P.Engels, C. Atherton and M.A. Hoefer, cond-mat/0701028. [30] L.N. Trefethen, Spectral Methods in MATLAB, SIAM(Philadelphia, 2000). [31] T. Christodoulakis, Th. Grammenos, Ch. Helias, P.G. Kevrekidis and A. Spanou,J. Math. Phys. 47, 042505 (2006).