Dynamical stabilization of matter-wave solitons revisited

1.5 MB·English

by  Alexander Itin

Dynamical stabilization of matter-wave solitons revisited
Alexander Itin1,2, and Shinichi Watanabe1, Toru Morishita1
1The University of Electro-Communications,
6
1-5-1 Chofu-ga-oka, Chofu-shi, Tokyo 182-8585,
0
0 Japan and 2Space Research Institute, RAS,
2
Profsoyuznaya Str. 84/32, 117997 Moscow, Russia
n
a
J
We consider dynamical stabilization of Bose-Einstein condensates (BEC) by time-

7
dependent modulation of the scattering length. The problem has been studied before by

]
r several methods: Gaussian variational approximation, the method of moments, method of
e
h
modulated Townes soliton, and the direct averaging of the Gross-Pitaevskii(GP) equation.
t
o
We summarize these methods and find that the numerically obtained stabilized solution
.
t
a
has different configuration than that assumed by the theoretical methods (in particular a
m
phase of the wavefunction is not quadratic with r). We show that there is presently no
-
d
n clear evidence for stabilization in a strict sense, because in the numerical experiments only
o
metastable (slowly decaying) solutions have been obtained. In other words, neither numer-
c
[
ical nor mathematical evidence for a new kind of soliton solutions have been revealed so

2
far. The existenceofthe metastablesolutionsis neverthelessaninterestingandcomplicated
v
2
phenomenon on its own. We try some non-Gaussian variational trial functions to obtain
7
4 better predictions for the critical nonlinearityg for metastabilization but other dynamical
cr
6
0 properties of the solutions remain difficult to predict.
5
0
/
t
a I. INTRODUCTION the NLSE, since it is known that NLSE (of-
m
ten called the Gross-Pitaevskii (GP) equation in
-
d
n The nonlinear Schrodinger equation (NLSE) that context) describes the dynamics of BEC at
o
c appearsinmanymodelsofmathematicalphysics zero temperature very well [1].
:
v
i and has numerous applications. The one- While early analytical studies of BECs were
X
r dimensional NLSE is famous due to its in- concentrated on (quasi-)one-dimensional sys-
a
tegrability and soliton solutions. The two-
tems, (quasi-)2D and 3D systems are more im-
dimensional and three-dimensional versions do
portant for real experiments. In 2D and 3D sys-
not have such properties and are much less ex-
tems analytical treatment of NLSE is very diffi-
plored.
cult and one has to use approximate methods.
In the last decade dynamics of BECs has at- One of the very interesting and complicated
tracted enormous amount of interest which in phenomena being studied recently is stabiliza-
turn is causing a renewed growth of interest in tion of BEC by the oscillating scattering length
2
in two and three dimensions. cluded from the dynamics, say, due to a tight
In 1D geometry, bright solitons can be stable confinement). The system is described by the
without trapping potential if nonlinearity is at- GP equation:
tractive and sufficiently strong. In NLSE with ∂ψ 1 ω2(t)
i = 2ψ+ r r2ψ+g(t)ψ 2ψ, (1)
attractive interaction (corresponding to BEC ∂t −2∇ 2 | |
withnegativescatteringlength)in2Dfreespace, where r2 = x2 + y2 and g(t) =
kinetic energy can balance interaction energy at (8πmω /¯h)1/2Na(t) describes the strength
z
certain critical value of nonlinearity gcr, but the of the two-body interaction. The interaction
resulting solution (Townes soliton) is unstable. g(t) is rapidly oscillating: g(t) = g +g sin(Ωt),
0 1
That is, if nonlinearity is either increased or de- while the confinement trap described by ω (t)
r
creased (and kept fixed afterwards), the solu- is slowly turned off. Refs. [4, 5, 6, 7] suggest
tion either expands or collapses correspondingly. it is possible to obtain a dynamically stabilized
It was shown by several authors that stabilized bright soliton in free space in such a way.
solutions are possible with the oscillating scat- Interactions between such objects were very
tering length. The oscillations of the scattering recently studied in Ref. [9]. This is a very
length lead to creation of pulsating condensate, interesting phenomenon not only in the context
i.e. some kind of breather solution. One can of BECs but also from a broader scope of
drawananalogy withKapitzapendulum(apen- nonlinear physics.
dulum with a rapidly oscillating pivot), where Suchkindof stabilization in3Dhas alsobeen
unstableequilibria ofunperturbedsystemis sta- reported [10]. The latter finding is, however, in
bilized by means of fast modulation. This idea some disagreement with other investigations on
was already applied to stabilization of beams in this topic (for example, Ref. [6]). In Ref. [11] it
nonlinear media [2]. Among many other appli- wasshownthatthescatteringlengthmodulation
cations in related fields, the atom wire trap sug- may indeed provide for the stabilization in 3D,
gested in Ref. [3] should be mentioned. In Refs. but only in combination with a quasi-1D peri-
[4, 5] the novel application of this stabilization odicpotential. So3Dgeometrymightneedaddi-
mechanism to BEC physics was presented which tional careful examination. In the present paper
in turn encouraged several other works on that weconcentrate onquasi-2Dcaseonly,wherealso
subject [6, 7, 8, 9, 10]. not everything is clear yet. Unlike conventional
We consider herethe problemof stabilization 1D solitons, higher-dimensional solitonic objects
of BEC in 2D free space by means of rapid os- may decay. Therefore, it is interesting to inves-
cillations of the scattering length in a greater tigate the following question: is there indeed a
detail (the third dimension is assumed to be ex- novel genuine breather solution behind the phe-
3
nomenon of stabilization? As we show in this lized wavefunction does not have such parabolic
paper, it turns out that the phenomenon does phase factor and does not fit into self-similar
not fit into simple models being suggested ear- patterns implied by the above-mentioned meth-
lier. For theoretical description of the process, ods. This qualitative difference between the ex-
several methods were used by different groups of act numerical solution and all theoretical mod-
authors: variational approximation basedonthe els considered so far was not mentioned earlier.
Gaussian anzatz [4, 6], direct averaging of the Besides, we noticed presence of steady outgoing
GP equation [6], a method based on modulated flux of atoms in numerical stabilized solutions.
Townes soliton [6], and the method of moments So, even numerically there is no 2D soliton so
[8]. Surprisingly,we findall themethods are not far, but some slowly decaying object instead.
verysatisfactoryevenforqualitativepredictions. Section 2 reviews the abovementioned theoret-
In brief, the direct averaging of the GP equation ical methods. In Section 3 we give some re-
has the disadvantage of ommiting terms which sults obtained using the variational approxima-
are of the same order as those responsible for tion with non-gaussian trial functions, including
creation of the effective potential, while three ”supergaussian anzatz”. It is shown that a bet-
other methods, although very different, all rely ter accuracy can be obtained for predicting crit-
on the unwarranted assumption of parabolic de- ical nonlinearity g , butwe were not able to de-
cr
pendenceofthephaseofthestabilizedwavefunc- termine accurately such dynamical properties as
tion on r: arg ψ = α(t)+β(t)r2. We find that the frequency of slow oscillations. Additionally,
the behavior of the exact numerical wavefunc- we checked the supergaussiananzatz for another
tion is, however, completely different (see Fig. problem: determination of critical number of at-
2). The above-mentioned parabolic approxima- tractive BEC ina parabolicwell, andfoundit to
tion (PA) of the phase factor is very popular be much more accurate than the usual gaussian
because it is appealingly simple and indeed of- anzatz. This example also demonstrates that
ten appears in solutions of the time-dependent the stabilization mechanism is essentially more
GP equation [13]. Usually it comes from self- complicated than that assumed by the present
similar time evolution of the condensate den- (PA-based) methods, because predictions of the
sity, for example in 3D the following dynamics supergaussian anzatz for dynamical properties
of the condensate density is possible ρ(x,y,z) = of the stabilized solution are much less accurate
[λ (t))λ (t)λ (t)] 1ρ(x/λ (t),y/λ (t),z/λ (t)) , than in static problems.
1 2 3 − 1 2 3
where coefficients λi are coupled by nonlinear In Section 4 numerical results are presented
differentialequations. Itistheimportantfinding andcomparedwithpredictionsofthetheoretical
of thepresent paperthatin our problemastabi- methods discussed in Sections 2 and 3. Configu-
4
ration of stabilized solution is discussed and dy- error, and similar values of discrepancy for other
namicsofsomeintegralquantitiesofthesolution dynamical quantities (as a useful test, in the
is investigated. Appendix we provide corresponding results ob-
In Section 5 concluding remarks are given. tained with a supergaussian variational ansatz).
Wementiontherelation betweentheBECstabi- However, it seems that in this example GA en-
lization problemand stabilization of optical soli- ables toreproduceimportantfeatures of thesys-
tons in a layered medium with sign-alternating tematleastqualitatively. TheGAwasalsoused
Kerr nonlinearity. in many other treatments of the GP equation
using a variational technique. In particular, it
II. SEVERAL APPROXIMATE was applied to the problem of BEC stabiliza-
METHODS TO STUDY THE PROBLEM: tionbytheoscillatingscatteringlength. TheLa-
PA-BASED METHODS (GAUSSIAN grangian density corresponding to the GP equa-
VARIATIONAL APPROXIMATION, THE
tion (1) is
MODULATED TOWNES SOLITON, THE
METHOD OF MOMENTS), AND THE
DIRECT AVERAGING OF THE GP i ∂ψ ∂ψ 1 ∂ψ 2 1
L[ψ] = ψ ∗ψ g(t)ψ 4.
∗
EQUATION. 2 ∂t − ∂t −2 ∂r −2 | |
(cid:18) (cid:19) (cid:12) (cid:12)
(cid:12) (cid:12) (2)
(cid:12) (cid:12)
(cid:12) (cid:12)
A. PA-based methods
The normalization condition for the wave-
1. Gaussian variational approximation function is 2π 0∞|ψ|2rdr = 1.
R
In Ref. [4], a variational method with the
Thevariational approach based on the Gaus- following Gaussian anzatz was used,
sian approximation (GA) is one of the most of-
1 r2 R˙(t)
ψ(r,t) = exp +i r2 ,
ten used in studying dynamics of the GP equa-
√πR(t) "−2R2(t) 2R(t) #
tion. In actual calculations this approximation (3)
however often gives a large error as compared where R(t) is the variational parameter that
to exact numerical results [7, 12]. For example, characterizes the size of the condensate, and the
in Ref. [12] the Gaussian approximation in dy- phase factor of the wavefunction describes the
namics of attractive BEC was compared to ex- mass current [4, 5, 16].
act numerical solution of the GP equation. It After substitution of expression (3) into the
was found that in estimating the critical num- Lagrangian density (2) one obtains the effective
ber Nc of the condensate (the maximal number Lagrangian L = 2π 0∞rL[ψ]dr and the corre-
of condensed particles in a trap before collapse sponding Euler-LagrRange equations of motion.
occurs)theGaussian approximation gives a 17% One can obtain then the equation of motion for
5
R(t) as dictions (5) and (6) based on the Gaussian ap-
proximation can catch (g /Ω)1/2 dependence of
1 g +g sinΩt 1
R¨(t) = + 0 1 . (4)
R3(t) 2πR3(t) the monopole moment < r > and (Ω/g ) de-
1
So the gist of the model is to represent the pendence of the breathing-mode frequency ω
br
2D BEC as a classical nonlinear pendulum with but cannot determine the corresponding coeffi-
modulated parameters. It is important that cients of proportionality, of which the one in (5)
other one-parameter PA-based anzatzes also becomes infinity while the one in (6) becomes
give the same nonlinear pendulum (R¨ = (a + zero for g = 2π, the value actually used in
0
−
bsinΩt)/R3, where a,b depend on the parame- the numerical calculations. On the other hand,
ters g ,g ,Ω), but with different functional de- from numerical calculations they were able to
1 0
pendence of a,b on the parameters. determine the coefficients as 1.06 and 0.32 cor-
The authors of Ref. [4] use then the Kapitza respondingly (see Fig. 2 of Ref. [4]). It was also
averaging method to study behavior of the determined in Ref. [4] that in order to stabilize
system with the rapidly oscillating scattering the bright soliton, g must exceed the critical
0
| |
length. They assume the dynamics of R can be value of collapse gcr . Their numerical estimate
| |
separatedintoaslowpartR0 andasmallrapidly for |gcr|is ≈ 5.8 whiletheoretical estimate based
oscillating component ρ: R = R (t) + ρ(Ωt). on Gaussian approximation is 2π 6.28. The
0
≈
From the equations of motion for R and ρ one 2π estimate in fact corresponds to fitting the so-
0
extracts the effective potential for the slow vari- called Townes soliton by a Gaussian trial func-
able U(R ) A2 + A6 and determines its mini- tion as will be discussed below.
0 ≈ R20 R60
mum Inspired by the idea of comparing a numer-
3 1/4 g 1/2 ical solution with simple model nonlinear pen-
1
R = − . (5)
min
(cid:18)4π(g0 +2π)(cid:19) (cid:18)Ω(cid:19) dulum, one may ask if it is possible to obtain
From the expression for the effective poten- a better accord with the numerical experiments
tial for R they obtained dependence of the using different ansatzes. We study this question
0
monopole moment < r > and the breathing- in Section 3, and it seems that only the station-
mode frequency ω on parameters g ,Ω. The ary Townes soliton can befit accurately, butnot
br 1
frequency of small oscillations (breathing mode) the stabilized breather solutions.
around the minimum is given by [4]
8Ω2
ω2 = (g +2π)2. (6) 2. Modulated Townes soliton
br 3g2 0
1
Their numerical calculations were done for A method based on modulated Townes soli-
g = 2π. One can see that theoretical pre- ton used in Ref. [6, 8] should be mentioned.
0
−
6
The Townes soliton is a stationary solution to where n = 2,3 is the dimension of the problem.
the 2D NLS equation with constant nonlinearity In 2D, dr = 2πrdr, and in 3D dr = 4πr2dr.
g . Inournotations g 1.862π 5.85. This For all t, we have I = 1. For the remaining
cr cr 1
| |≈ ≈
solution is unstable: if g is slightly increased I one can write down the dynamical equations
i
| |
or decreased, the solution will start to collapse of motion as [8]:
or expand correspondingly. If the value of g is
close to gcr, one may search for a solution of the I˙ = I , I˙ = I , I˙ = gn−2I +g˙I , (9)
2 3 3 4 4 5 5
n
problem with fast oscillating g in the form of a
nπ2n ∂ ψ 4 ∂argψ
modulated Townes soliton, as described in Refs. I˙5 = ∞ | | rn−1dr(10)
8 ∂r ∂r
Z0
[6, 8]. A solution is sought in the form of
The system of equations for the momenta is not
closed because of I , and one should make some
Ψ(r,t) [a(t)] 1R [r/a(t)]eiS, 5
− T
≈
r2a˙ approximation in order to close it. In Ref. [8] it
S = σ(t)+ , σ˙ = a 2, (7)
−
4a was assumed that
where R represents amplitude of the Townes I r2
T 3
argψ = , (11)
4I
soliton. Then, starting from the approxima- 2
i.e. thephasefactorisproportionaltor2(sothat
tion (7), one can derive the evolution equation
again it is a PA-based method) and the coeffi-
for a(t) and so determine the dynamics of the
cient of proportionality is given by the ratio of
system. Note that the approach is also PA-
I and I . Then the system (9) poses dynamical
based. It is inevitable if we are to use one- 3 2
invariants [8]:
parameter self-similar trial function in the form
of |ψ(r,t)| = Af(r/a,t). Q1 = 2(I4−gI5)I2− 14I32, (12)
n/2
Q = 2I I . (13)
2 2 5
3. The method of moments
With the help of these invariants, the system
Another PA-based method we would like to becomes
mention hereis themethodof moments [8]. One I¨ 1 I˙ 2 = 2 Q1 +g Q2 . (14)
introduces integral quantities I ,I ,I ,.. as 2− 2I2 2 I2 In/2!
1 2 3 (cid:16) (cid:17) 2
Introducing X(t) = I (t) one obtains [8]
2
I = ∞ ψ 2dr, I = ∞r2 ψ 2dr,
1 2
Z0 | | ∂ψ ∂ψZ0 | | X¨ = Qp1 +g(t) Q2 . (15)
I3 = i ∞ ψ ∗ ψ∗ rdr, (8) X3 Xn+1
∂r − ∂r
Z0 (cid:18) (cid:19)
1 n The equation is analogous to that obtained
I = ∞ ψ 2 + g(t)ψ 4 dr,
4
2 |∇ | 2 | |
Z0 (cid:18) (cid:19) by other PA-based methods. One can investi-
n
I = ∞ ψ 4dr,
5 4 | | gate the obtained equation (15) using various
Z0
7
methods of nonlinear dynamics. The simplest of optical solitons) [14]. In Ref. [6] the solution
Kapitza averaging method can be used again, is sought as an expansion in powers of 1/Ω (in
but of course it is better to use rigorous averag- our notation):
ing technique since modern averaging methods
ψ(r,t) = A(r,T )+Ω 1u (A,ζ)+Ω 2u (A,ζ)+....,
k − 1 − 2
are available [20] which have been extensively
(16)
used already in plasma physics, hydrodynamics,
with < u >= 0, where < ... > stands for the
k
classical mechanics [21]. The authors of Ref. [8]
average over the periodof therapid modulation,
fulfilled rigorous analysis of model Eq. 15 us-
T Ω kt are the slow temporal variables (k =
k −
ing results of Ref. [15]. It is important to have ≡
0,1,2,...), while the fast time is ζ = Ωt. Then,
in mind that the relation between the exact dy-
for the first and second corrections the following
namics of the full system and that of the model
formulas were obtained:
(15) of the method of moments remains unclear,
therefore one cannot determine sufficient condi- u = i[µ < µ >]A2A,
1 1 1
− − | |
ζ
tions for stabilization, etc. In Ref. [8] it was no-
µ [g(τ) < g >]dτ, (17)
1 1
≡ −
ticedthatthecorrespondencebetweennumerical Z0
u = [µ < µ >][2iA2A +iA2A +∆(A2A)]
2 2− 2 | | t ∗t | |
simulation of full 2D GP equation and dynam- 1
A4A( [(µ < µ >)2 2M]+
1 1
ics of the model system (15) is not good. As it − | | 2 − −
+ < g >(µ < µ >))],
2 2
is seen from Fig. 3 of Ref. [8], neither the fre- −
ζ 1
µ = (µ < µ >)ds, M = (< µ2 > < µ >2).
quencyofslowoscillationsnorthepositionofthe 2 1− 1 2 1 − 1
Z0
minimum of the effective potential is predicted
Using these results, the following equation was
correctly. Nevertheless, we foundthat in numer-
obtainedfortheslowlyvaryingfieldA(r,T ),de-
0
ical stabilized solutions magnitudes of Q and
1 rived up to the order of Ω 2:
−
Q are often well-conserved, i.e. they oscillate
2
∂A g 2
i = ∆A+ A2A+2M 1 [A6A
about some mean value (see Section 4).
− ∂t | | Ω | |
(cid:18) (cid:19)
3A4∆A + 2A2∆(A2A)+A2∆(A2A )].
∗
− | | | | | | | |
(18)
B. Direct averaging of the GP equation
Theaboveequationwasrepresentedinthequasi-
Ref. [6]alsoexplorestheGaussianvariational
Hamiltonian form
approximation. Beside that, a very promising
g 2 ∂A δH
method of directly averaging the GP equation 1 + 6M 1 A4 = q,
Ω | | ∂t −δA
h (cid:18) (cid:19)i ∗
was investigated. It is based on an analogous g 2
H = dV A2 2M 1 A8
q
|∇ | − Ω | |
methodusedfortheone-dimensionalNLSEwith Z h (cid:18) (cid:19)
1 g
A4+4M 1 (A2A)2)2 . (19)
periodically managed dispersion (in the context
− 2| | Ω|∇ | | |
i
8
However, somecontributionwasmissedwhile III. VARIATIONAL APPROXIMATION
deriving Eq. (18). Let us take into account the WITH NON-GAUSSIAN ANSATZES
third correction u (A,ζ):
3
Here we try to investigate the system more
accurately using some non-Gaussian ansatzes
and see if it is possible to get more accurate
ψ(r,t) = A(r,T )+Ω 1u (A,ζ)+Ω 2u (A,ζ)
k − 1 − 2 theoretical estimates. One may be interested in
+ Ω−3u3(A,ζ)+.... (20) three dynamical quantities of the system: the
value of critical nonlinearity g , slow frequency
cr
of breathing oscillations of the stabilized soli-
ton ω , and minimum of the effective potential
Then, up to terms of order Ω 2 it changes noth- br
−
R about which the expectation value of the
min
ing in r.h.s of Eq.(18) (spatial part), but it adds
monopole moment < r > oscillates slowly.
to l.h.s. of Eq. (18) an undetermined term
Ω 2∂u /∂ζ . This term has the same order Ω 2 Table1summarizesresultsofvariational pre-
− 3 −
dictions for the critical nonlinearity g and fre-
as the terms from the second correction. So we cr
quency of small breathing oscillations using sev-
donotgethereaconsistentequationfortheslow
eral different ansatzes. Note that the phase de-
field A because we do not have a closed set of
pendence of a one-parameter trial function is
equations for thesecond-order corrections (third
not important for calculating g . It is under-
ordercorrectionbecomessecondordercorrection cr
stood that if we choose a trial wavefunction
after differentiating in time), and so the quasi-
with its amplitude in the form of ψ(r,t) =
Hamiltonian (19) contains an undetermined er-
| |
ror of the second order in Ω 1. The influence Af[r/a(t)], then we need to use a phase factor
−
with quadratic r dependence in order for the
of the contribution is not very clear but require
−
ansatz to be self-consistent (i.e., the mass cur-
additional investigation. Nevertheless, formally
rentgenerated bythechangingparameter would
the omitted terms have the same order as those
beincorporatedinthephasefactorofanansatz).
responsible for the creation of the effective po-
On the other hand, since amplitude part of the
tential. Having in mind how many difficulties
trial function is just an approximation, one may
ariseinaveragingofsystemsofordinarydifferen-
try to use other forms of phase factor with the
tial equations [20], the rigorous direct averaging
same functional form of the amplitude.
oftheGPequation constitutes averyinteresting
and challenging open problem, since in principle When predicting the frequency of breathing
it could reveal a true periodic solutions in such oscillations from the corresponding effective po-
oscillating objects. tential, it is easy to obtain the result for small
9
amplitude linear breathing oscillations (given in g just because it is decaying too fast at large
cr
Table1),butinactualstabilizedsolutionsampli- r. The supergaussian trial function provide a
1
tudes of breathing oscillations are not so small. better approximation, namely gcr = π2ln2ln2
It is possible to take into account anhar- which corresponds to the supergaussian wave-
monicity of breathing oscillations. As was men- function with η = ηT = 2ln2 < 2. Previously
tioned earlier, all PA-based anzatzes produce the supergaussian ansatz was used to fit station-
the nonlinear pendulum R¨ + (a + bsinΩt)/R3, ary solutions of some nonlinear problems includ-
with a corresponding effective potential having ing NLS equation in the context of BECs [19].
R = 3b2 1/4 , ω = 8Ω a/b, where The superposition of two Gaussians in the form
min −2Ω2a br 3 | |
ωbr is th(cid:16)e frequ(cid:17)ency of the qsmall amplitude Aexp(−2rR22)Cosh(γ2rR22) also enables one to ob-
breathing oscillations (near the bottom of tain some improvement: gcr 5.883. The Se-
≈
the effective potential). For larger breathing canth ansatz
oscillations the (anharmonic) breathing fre- A
ψ = exp[iS(R˙,R)r2]
cosh(r/R)
quency will be amplitude-dependent: ωanh =
br
1 works better, with only one parameter it over-
2π 2 x3 K(k)+ x2√x x E(k) − ,
−h √x2 x3 x1 2− 3
wit(cid:16)hq k =h −x2 x1, where x = R2, x i=(cid:17)R2 comes the above-mentioned two-parameter trial
x2−x3 1 1 2 2
q − functions. A very good approximation is pro-
(R ,R being the turning points), x is the
1 2 3
vided by the simplest ansatz among all consid-
third root of the equation h = a + b . The
2x 4Ω2x3
ered:
magnitudes of x ,x ,x ,h can be determined
1 2 3
1 r r
fromnumerically obtained breathingoscillations ψ = 1+ exp +iS(R˙,R)r2 .
3R√π 2R −2R
(cid:18) (cid:19) (cid:26) (cid:27)
(but results depend on the choice of a particular (21)
anzatz). Even this improvement is not helpful, It fits the Townes soliton adequately both at the
simply because the parabolic approximation is origin and asymptotically at infinite r ( a pre-
not valid. exponential multiplier is notso important as the
Finding g only might be considered as an exponential factor ). The pre-exponential factor
cr
approximation to the stationary Townes soli- is needed in order to fulfill the boundary condi-
ton by a trial function so that the mass cur- tion in the origin lim 1ψ < . Note that
r→0 r r ∞
rent term equals zero and that a phase factor inthesupergaussianansatztheformercondition
may be skipped from the calculations. It is is not fulfilled, otherwise (if one included it in
known that the Townes soliton ψ = eitR (r,t) a similar way) the result would be better at the
t T
at large r has asymptotic behavior for its am- cost ofmorebulkycalculations. Theaccuracy of
plitude in the form R e r/√r . So that thepredictionimpliesthatansatz(21)providesa
T −
∼
Gaussian ansatz is not very good for finding very good approximation to the Townes soliton
10
at fixed R, and could approximately represent sion for the phase factor. We also try the su-
the modulated Townes soliton when R is time- pergaussian ansatz with fixed η and with non-
dependent and the phase factor with parabolic quadratic phase dependence (which is unfortu-
r dependence is used in accordance with the natelynotself-consistenttrialfunction)ψ(r,t) =
−
continuity condition. Aexp (a+ib)rηT , where A,a,b, and η are all
− 2
h i
After obtaining estimates for g , one can use functions of time, parameter η is fixed at the
cr
theabove-mentioned ansatzesinordertofindan valueofitsTownessoliton-likesolutionη = ηT =
effective potential, its minimum and frequency 2ln2. We find that such modification drasti-
of the breathing oscillations of the monopole cally changes dynamical parameters of the sys-
moment about this minimum in the same way tem. Still, the resulting model is the same clas-
as it was done for the Gaussian ansatz. We sical nonlinear pendulum as in the Gaussian ap-
checked the Sech ansatz and the supergaussian proximation, butwithdifferentparameters. The
with quadratic phase dependence. In the su- rigorous way to employ the two-parameter su-
pergaussian ansatz the parameter η was fixed pergaussian ansatz is to let η be a dynamical
at the value of its ”Townes soliton-like” solu- variable and construct a phase factor fulfilling
tion η = η = 2ln2. In such a way the varia- continuity condition for the trial function. One
T
tional approximation with supergaussian ansatz could then obtain the two-dimensional effective
resembles method of modulated Townes soli- potential within the same Kapitza approach.
ton. However, we find that such trial function As a useful test of applicability of the su-
seriously underestimate minimum of the effec- pergaussian anzatz, we determine the critical
tive potential (i.e. the mean value about which number of attractive BEC in the 3D parabolic
the monopole moment oscillates). Nevertheless, trap studied in Ref. [12]. Their numerical re-
the result of the Gaussian ansatz is even worse sult was Ncr = 1258.5, while the gaussian ap-
since for g = 2π it gives the diverging expres- proximation yields NG = 1467.7. We found
0 cr
sion for R and zero for frequency of slow thesupergaussianprediction tobevery accurate
min
breathing oscillations ω , as mentioned in Sec- NSG = 1236.1.
br cr
tion 1 and [4]. A natural idea for remedy is
to use two-parameter trial functions to repro-
IV. NUMERICAL RESULTS
duce the non-parabolic phase factor dependence
on r. In the supergaussian ansatz it can be Numerical calculations reveal the fact that
done by considering η as a dynamical (time- stabilized solutions do not have parabolic phase
dependent) parameter. The problem is that it factors in contradiction to all the methods con-
is difficult to obtain the self-consistent expres- sideredinSection 2(except themethodof direct