Dark spherical shell solitons in three-dimensional Bose-Einstein condensates: Existence, stability and dynamics

0.99 MB·

by  Wenlong Wang

Dark spherical shell solitons in three-dimensional Bose-Einstein condensates:
Existence, stability and dynamics
Wenlong Wang,1,2,∗ P.G. Kevrekidis,3,4,† R. Carretero-Gonza´lez,5 and D. J. Frantzeskakis6
1Department of Physics and Astronomy, Texas A&M University, College Station, Texas 77843-4242, USA
2Department of Physics, University of Massachusetts, Amherst, Massachusetts 01003 USA
3Department of Mathematics and Statistics, University of Massachusetts, Amherst, Massachusetts 01003-4515 USA
4Center for Nonlinear Studies and Theoretical Division,
6
Los Alamos National Laboratory, Los Alamos, NM 87544
1 5Nonlinear Dynamical Systems Group,‡ Computational Sciences Research Center,
0
and Department of Mathematics and Statistics, San Diego State University, San Diego, California 92182-7720, USA
2
6Department of Physics, University of Athens, Panepistimiopolis, Zografos, Athens 15784, Greece

n
a Inthisworkwestudysphericalshelldarksolitonstatesinthree-dimensionalatomicBose-Einstein
J condensates. Their symmetry is exploited in order to analyze their existence, as well as that of
topologically charged variants of the structures, and, importantly, to identify their linear stability
0
1 Bogolyubov-de Gennes spectrum. We compare our effective 1D spherical and 2D cylindrical com-
putations with the full 3D numerics. An important conclusion is that such spherical shell solitons

] can be stable sufficiently close to the linear limit of the isotropic condensates considered herein.
s We have also identified their instabilities leading to the emergence of vortex line and vortex ring
a
cages. In addition, we generalize effective particle pictures of lower dimensional dark solitons and
g
ring dark solitons to the spherical shell solitons concerning their equilibrium radius and effective
-
t dynamics around it. In this case too, we favorably compare the resulting predictions such as the
n
shell equilibrium radius, qualitatively and quantitatively,with full numerical solutions in 3D.
a
u
q
I. INTRODUCTION such structures are radially symmetric solitons, such as
.
t
a theringdarksolitons(RDS),initiallyproposedinatomic
m BECs in Ref. [26] (see also Refs. [27–29] for subsequent
The pristine setting of atomic Bose-Einstein conden-
studies). Nevertheless, such structures were studied to a
- sates(BECs)has offereda significantplaygroundfor the
d far lesserextend, arguablydue to their generic identified
examination of numerous physical concepts [1, 2]. One
n instability, which was found to give rise to states such
o of the flourishing directions has been at the interface of
as alternating charge vortex polygons. Notice that, as
c the theory of nonlinear waves and such atomic (as well
was recently shown in Ref. [30], stabilization of RDS is
[ as optical) systems, concerning, in particular, the study
possible upon employing a proper potential barrier.
of matter-wave solitons [3, 4]. Such coherent nonlinear
1
v structureshavebeenobservedinexperimentsandstudied The three-dimensional non-planar generalization of
6 extensivelyintheory. Prototypicalexamplesofpertinent darksolitons,namelysphericalshellsolitons,hasreceived
7 studies include —but are not limited to— bright [5–7], very limited attention. The analysis of Ref. [27] in the
1 dark [8] and gap [9] matter-wavesolitons, as well as vor- case of box potentials (rather than the typical parabolic
2 tices [10, 11], solitonic vortices and vortex rings [12]. trap setting of BECs) revealed their potential existence,
0 butdidnotpursuetheirstabilityanalysis. Moreover,the
. Among these diverse excitations, dark solitons in re-
1 adiabatic dynamics of these structures was studied ana-
pulsive BECs have enjoyed a considerable amount of at-
0 lytically via a variational approach in Ref. [31] without
tentionin effectively one-dimensionalsettings due to nu-
6 providing, however, relevant numerical results. Never-
merousexperimentsleading(especially,morerecently)to
1
theless, it is important to note that in an experimental
: their well-controlled creation [13–18]. Moreover, numer-
v realizationofvortexringstates,thetransientobservation
ousworksconsidereddarksolitonsinhigher-dimensional
i of spherical shell solitons was reported [32]. In addition,
X settingsinordertoexperimentallyexploretheirinstabil-
they were also suggested as a potential outcome of colli-
ity leading to the formation of vortex rings and vortex
r
sions of vortex rings in the work of Ref. [33].
a lines, as illustrated, e.g., in Refs. [19–23]. More recently,
suchexcitationshavealsobeenobservedinfermionicsu- It is the purpose of the present manuscript to re-
perfluids [24, 25]. visit these states and explore not only their existence
in isotropic traps, but also perform their systematic
On the other hand, multi-dimensional (i.e., non-
Bogolyubov-de Gennes (BdG) stability analysis. This is
planar) variants of dark solitons, were studied as well
done by projecting perturbations onto eigenfunctions of
—cf. e.g., Chapters 7 and 8 of Ref. [3]. Examples of
the linear problem in the angular directions. Using this
basis, the relevant computations are greatly simplified,
and the full 3D spectrum computation is reduced to an
effectivelysmallsetofone-dimensionalproblems(forthe
‡URL:http://nlds.sdsu.edu
∗Electronicaddress: [email protected] different modes). It should be mentioned here that the
†Electronicaddress: [email protected] use of sphericalharmonicsas a basis in the angularvari-
2
ables for states bearing radial symmetry is rather com- z, respectively. Note that the potential has rotational
mon in terms of numerical schemes for linear and non- symmetry with respect to the z-axis. In our numerical
linear Schr¨odinger (NLS) equations —see, e.g., Ref. [34] simulations, we focus on the fully symmetric (isotropic,
foraarecentexample. Earliereffortsalongtheselinesin spherically symmetric) case with ω = ω = ω = 1, in
ρ z
thecontextoflinearSchr¨odingerequationscanbefound, which we can benchmark our numerical methods in ef-
e.g., in the works of Refs. [35, 36]. In 2D BECs, analo- fective 1D and 2D against those of the fully 3D setting.
gousdecompositionsofazimuthalmodesforradialstates Weexaminethedarksphericalshellsolitons,hereafter
including RDSs and vortices [27, 28, 37] have appeared, referred to as DSS. This single radial node state exists
andin3Dsimilarpossibilitiesarestartingtoemerge[38]. in this isotropic case from the linear limit onwards as a
In this work, we find that indeed spherical shell soli- stationary state of the form ψ(x,y,z,t) = e−iµtψ (r)
DSS
tonscanbeidentifiedasstable3DsolutionsoftheGross- wherer = x2+y2+z2 isthe sphericalradialvariable.
Pitaevskiiequation,sufficientlyclosetothelinear(small- At the linear limit, the waveformis an eigenmode of the
p
amplitude) limit of small chemical potential from which quantumharmonicoscillatorwithchemicalpotentialµ=
they bifurcate, as is shown below. Moreover,we provide 7ω/2, and spatial profile:
atheoreticalanalysisoftheir dynamicsas“effective par-
1
ticles”. Thisapproachallowstocharacterizetheopposite ψ = (200 + 020 + 002 ) (3)
DSS linear
limit of large chemical potential. Our numerical compu- | i √3 | i | i | i
tations of existence and stability, complemented by di- ωr2 3 e−ωr2/2, (4)
rect numerical simulations of the spherical shell soliton ∝ − 2
(cid:18) (cid:19)
dynamics,notonlycorroboratetheabovetwolimits,but
wherethe basis n n n denotesthe 3Dharmonicos-
alsoprovideasystematicwaytointerpolatebetweenthe x y z
{| i}
cillator quantum states in Cartesian coordinates.
two analytically tractable regimes.
The DSS is generally expected to be unstable for
Our presentation is organized as follows. First, in
high chemical potentials due to transverse modulational
Sec. II, we introduce the model and describe the imple-
(snaking-type) instability, in a way similar to its pla-
mentationof the linear stability analysis. Ouranalytical
nar and RDS counterparts, as summarized, e.g., in
approaches for the spherical shell dark solitons in the
Refs. [4, 8]. However, for smaller chemical potentials,
above mentioned analytically tractable limits are then
and especially near the linear limit, it is relevant to an-
presentedinSec.III.Next,inSec.IV,wepresentoursys-
alyze stability of the DSS and identify corresponding in-
tematic numerical results. Finally, our conclusions and
stabilities and where they may arise. In what follows
a number of open problems for future consideration are
(cf. Sec. IV below), we will turn to numerical computa-
giveninSec.V.IntheAppendix, asimilarnumericalde-
tions, based on a fixed point iteration and our spherical
composition method but using cylindrical (rather than
harmonic decomposition of the spectral stability BdG
spherical) symmetry is also discussed whose results are
problem, in order to efficiently identify the intervals of
also given for comparison in Sec. IV.
existence and stability of the relevant modes.
II. MODEL AND COMPUTATIONAL SETUP
B. The Linear Stability Problem: Spherical
Harmonic Decomposition
A. The Gross-Pitaevskii equation
In this section, we discuss our analysis of the linear
In the framework of lowest-order mean-field theory,
stability of the DSS states utilizing their effective 1D
and for sufficiently low-temperatures, the dynamics of a
nature, namely their spherical symmetry. In spherical
3D repulsive BEC, confined in a time-independent trap
coordinates (r,θ,φ), the Laplacian can be decomposed
V, is described by the following dimensionless Gross-
into radial and angular parts, denoted by ∆ and ∆
R S
Pitaevskii equation (GPE) [1–4]:
respectively, namely
iψt =−12∇2ψ+Vψ+|ψ|2ψ−µψ, (1) ∇2f =∆Rf + ∆rS2f, (5)
where
whereψ(x,y,z,t)isthe macroscopicwavefunctionofthe
BEC and µ is the chemical potential (subscripts denote 1 ∂ ∂f
∆ f = r2 , (6)
partial derivatives). Here, we consider a harmonic trap R r2∂r ∂r
(cid:18) (cid:19)
of the form: 1 ∂ ∂f 1 ∂2f
∆ f = sinθ + , (7)
1 1 S sinθ∂θ ∂θ sin2θ∂φ2
V = ω2ρ2+ ω2z2, (2) (cid:18) (cid:19)
2 ρ 2 z The∆ operatorhaseigenstatesgivenbythespherical
S
harmonics Y with eigenvalues ℓ(ℓ+1), i.e.,
where ρ = x2+y2, ω and ω are the trapping fre- { ℓm} −
ρ z
quencies along the (x,y) plane and the verticaldirection ∆ Y = ℓ(ℓ+1)Y . (8)
S ℓm ℓm
p −
3
For states with sphericalsymmetry (ℓ=0), the ∆ part
S
TABLEI: Relevantstates foreigenvalues nearthespectrum
of the Laplacian is not relevant for identifying the sta-
of Im(λ)=0,1 and 2 for the degenerate perturbation theory.
tionary state. Thus, the stationary DSS state ψ (r) of
0 Here mnp stands for eigenstates with all possible distinct
Eq. (1), defined in the domain r [0, ), satisfies the | i
permutations with thequantumnumbersm,n and p.
∈ ∞
following radial equation:
Im(λ) “up”states “down” states
1
∆ ψ +V(r)ψ + ψ 2ψ µψ =0. (9) 0 002 011 002 011
R 0 0 0 0 0 | i | i | i | i
− 2 | | − 1 003 012 111 001
| i | i | i | i
Now, let us consider the linear stability (BdG) spec- 2 004 013 022 112 000
| i | i | i | i | i
trum of such a stationary state. We consider small-
amplitude perturbations expanded using the complete
basis of Y (the spherical harmonic eigen-basis) as mentioning that a large class of BEC states belongs to
ℓm
{ }
follows: this class.
Finally, we discuss the application of the degenerate
ψ(r,t)=ψ0+ [aℓm(r,t)Yℓm+b∗ℓm(r,t)Yℓ∗m], (10) perturbationmethod (DPM) [39] to the DSS state. This
ℓm method is employed for the determination of the DSS’
X
spectrum near the linear limit, and also for the identi-
where the space- and time-dependent coefficients
fication of the nature of instabilities. The method can
a (r,t) and b (r,t) determine the radial and tempo-
ℓm ℓm be most readily understood from the 3D version of the
ral evolution of the perturbation. Then, substituting
matrix M, where one can separate the free (linear) part
Eq. (10) into Eq. (1) and retaining up to linear terms
with known eigenstates and eigenvalues and the rest of
in the expansion,one notes that all the ℓ modes are mu-
thetermscanbetreatedasperturbationsnearthelinear
tually independent due to the sphericalsymmetry of the
limit. Note that there are 2 2 block matrices within M
steady state. For mode ℓ, the coefficients a and b of the ×
and the complete basis consists of both “positive” and
expansion obey the following evolution equations:
“negative” eigenstates. This terminology is associated
with quantum harmonic oscillator eigenstates leading to
ia = 1∆ a+ ℓ(ℓ+1)a+Va+2ψ 2a+ψ2b µa,
t −2 R 2r2 | 0| 0 − positive or negative eigenvalues for the linear analogs
−ibt = −21∆Rb+ ℓ(2ℓr+21)b+Vb+2|ψ0|2b+ψ0∗2a−µb, 0o.f tIhneDoPpeMra,ttohres fMuniictiinonEaql.d(e1p2e)ndi.een.,cewoitfhth|ψe0e|2igesentvatlo-
where we have dropped for simplicity the subscripts ues on the chemical potential µ is dominated by the in-
(ℓ,m). Asweareinterestedinthestabilityofthesteady- terplay of equi-energetic, degenerate states. The result
state, we now expand the coefficients using of the nonlinear perturbations ( ψ 2) in M within
0 ij
∝ | |
Eq.(12) is to cause the degenerateeigenvalues to depart
a = a0 eλt, (11) from their respective linear limit, as the norm of the so-
b b0 lution (mass) increases. The free part provides the basis
(cid:18) (cid:19) (cid:18) (cid:19)
that is formed from “up” states and “down” states, of
where λ is an eigenvalue of the following matrix
n n n 0
the form | x y zi and , with eigenval-
0 n n n
M = M11 M12 , (12) ues E (cid:18) E (cid:19)and E(cid:18)| x Ey zi(cid:19) respectively,
(cid:18)M21 M22(cid:19) where|nExΛnyinsztih−eeigDeSnSenergyDofSSth−e3D|nxqnuyannztiumharmonic
oscillator for eigenstate Λ. The relevant states for eigen-
with
values near the spectrum of Im(λ)=0, 1 and 2 are listed
M = i 1∆ + ℓ(ℓ+1) +V +2ψ 2 µ , in Table. I. Then, the effect of the nonlinearity on these
11 − −2 R 2r2 | 0| − eigenmodesandthecorrespondingeigenvaluecorrections
M = iψ(cid:16)2, (cid:17) are calculated, as discussed in Ref. [39] (and originally
12 − 0
spearheaded in this context in Ref. [40]). Below, this is
M = iψ∗2,
21 0 compared to the numerical spectrum in Sec. IV.
M = i 1∆ + ℓ(ℓ+1) +V +2ψ 2 µ .
22 −2 R 2r2 | 0| −
(cid:16) (cid:17)
One can, therefore, compute the full spectrum by com- III. THE PARTICLE PICTURE FOR THE DSS
puting that of each mode ℓ independently, and then
putting them all together to “reconstruct” the full spec- We now turn to the limit of large chemical poten-
trum. In what follows we use ℓ=0,1,2,...,10 for values tial µ. In this so-called Thomas-Fermi (TF) limit of
of the chemical potential up to µ 12. The reduction largeµ, the groundstate ofthe GPE is approximatedas
≈
of dimensionality with the decomposition on a complete ψ = max(µ V,0)[1,2]. Anaturalwayto obtaina
TF
−
basis set is also availablein 2D for states with rotational reduceddynamicaldescriptionofthe DSS —analogously
p
symmetry up to a topological charge along the z axis. to what is done in lower-dimensionalsettings— is to de-
For clarity, this is shown in the Appendix. It is worth velop an effective particle picture for the DSS in the TF
4
10
Ψ 0
−10
0 5 10 15
r
5
x 10
2
N
1
0
0 20 40 60 80
µ
FIG.2: (Color online) Top panel: aradial profileofthedark
soliton shell state in the TF limit for µ=80. Bottom panel:
number of particles as a function of chemical potential µ,
showing the continuation of states from the linear limit to
thenonlinear regime.
andprimesdenotederivativeswithrespecttor. Inspired
by the profile of the 1D and 2D equivalents of the DSS,
namely the darksolitonand the RDS, we seek a station-
ary DSS solution in the form of q(r)=tanh[õ(r rc)].
−
Multiplying both sides of Eq. (13) by q′ and integrat-
ing over r from to (bearing in mind that the
−∞ ∞
contribution of the integral from r = to r = 0 is
−∞
exponentiallysmall), we find thatthe equilibriumradius
is given by
FIG. 1: (color online) The two pairs of the top panels show
the modulus Ψ (left panels) and argument κ (right panels) √αµ
of adarksolit|on|shellstate at µ=12; thetoprow illustrates rc = , (14)
ω
the z =0 plane, while the middle row the x=0 plane. The
bottom panel shows an isocontour plot of thesame state. where α=5 √17 0.8769. This result, as in the case
− ≈
of the dark soliton ring in 2D [30], slightly overestimates
the actualvalue for α that, as we will see below,is accu-
ratelypredictedbythe secondparticlepicture approach,
limit. Here, following Ref. [30], we focus on the equi-
which retrieves the precise value for α to be α = 4/5;
librium radius r of the DSS, as a prototypical diagnos-
c this is confirmed by the numerical results shown below,
tic that we can compare to our numerical results. This
in Sec. IV.
equilibrium radius is determined as the critical value of
Thesecondapproachreliesonenergyconservationand
the DSS radius for which the restoring force due to the
it is based on the analysis of Ref. [42]. In this approach,
harmonic trap is counterbalanced by the effective force
it is argued that the equation of motion can be derived
exerted due to the curvature of the DSS.
by a local conservation law (i.e., an adiabatic invariant)
The first approach (motivated also by the work of
in the form of the energy of a dark soliton under the
Ref.[41])focusesonthestationarystateusingtheansatz
effect of curvature and of the density variation associ-
ψ(r)=ψ (r)q(r),whereoneobtainsanequationforthe
TF ated with it. More specifically, knowing that the en-
DSS profile q(r) (on top of the TF ground-state ψ )
TF ergy of the 1D dark soliton centered at x is given by
c
given by: = (4/3)(µ x˙ )3/2 [8], the generalization of the rele-
c
E −
vantadiabaticinvariantquantityina3Ddomainbearing
1
q′′+µq(1 q2)=P(r), (13) density modulations reads:
2 −
= 4πr2 4(µ V(r) r˙2)3/2
where E 3 − − (15)
2 ψ′′ ψ′ 1ψ′ = 4πr02(cid:2)43(µ−V(r0)−r˙02)3/2(cid:3) ,
P(r)=Vq(1 q2) q′ TF q TFq′ TFq,
− − r − 2ψ − ψ − rψ wherer(t) andr˙(t) ar(cid:2)ethe (radial)location(cid:3)andvelocity
TF TF TF
5
of the DSS with initial conditions r0 and r˙0. Taking a 3 (a) 1D/2D
time derivative on both sides of Eq. (15) and assuming
that the DSS has no initial speed (r˙ = 0), we obtain
0
a Newtonian particle equation of motion for the DSS of 2
the following form: λ
1∂V 2 r 4/3 1
0
r¨= + [µ V(r )]. (16)
0
−2 ∂r 3r r −
(cid:16) (cid:17)
Fromthisequation,the equilibriumpositionforthe DSS 0
4 6 8 10 12 14
4µ µ
yields rc = /ω,a resultconsistent with the numeri- 2
r 5 (b) 1D/2D
cal findings of the next Section. We now proceed to test
these predictions, the stability analysis spectrum, and 1.5
the equilibrium positions for the DSS.
λ 1
IV. NUMERICAL VS. ANALYTICAL RESULTS
0.5
Let us start by providing some of the basic features of 0
theDSS.AtypicalDSSstateatµ=12isshowninFig.1, 2
(c) 3D
where its spherical symmetry is apparent. A radial plot
ofthewavefunctionintheTFlimitisdepictedinthetop 1.5
panelofFig.2. Asthechemicalpotentialisincreasedthe
number of atoms (“mass”)
λ 1
∞
N = ψ 2dxdydz =4π r2 ψ 2dr, (17) 0.5
| | | |
Z Z0
increases as depicted in the bottom panel of Fig. 2. 0
4 5 6 7
We now examine the stability of the DSS. The corre- µ
sponding spectra illustrating the imaginary part of the
eigenvalues [dark (blue) lines) and the real part [light FIG. 3: (Color online) Stability spectra for the dark soliton
shellasafunctionofthechemicalpotentialµ. (a)-(b)Spectra
(orange) lines] are shown in Fig. 3. Panel (a) depicts
computedusingoureffectively1Dmethod(intheradialvari-
the spectrum obtained via our numerical method over a
able, using a spherical harmonic decomposition in the trans-
wide range of chemical potentials. Identical results were
verse directions). Identical results are obtained when using
obtainedwhenusingthemethodin2D,basedonacylin-
the2D variant of the method in cylindrical coordinates (and
drical coordinate decomposition and a representation of
decomposing onlytheazimuthaldirection insuitable Fourier
the azimuthal variable dependence in the corresponding modes). Both1Dand2Dnumericsuseaspacingofh=0.01.
Fourier modes. For reasons of completeness, a summary (c)Spectracomputedusingthefull3Dnumericswithspacing
ofthisvariantofthemethodispresentedintheappendix. h = 0.2. Orange (grey) lines depict the real part of λ (i.e.,
Panels (b) and (c) depict a zoomed in regionto contrast unstablepartsoftheeigenvalues),whileblue(dark)linesthe
theresultsbetweentheselowerdimensionalmethodsand imaginary part of λ.
the full 3D numerics. Panel(b) corresponds to the spec-
trum as computed by means of the 1D spherical coor-
dinate decomposition method using spherical harmonics
with spatial spacing of h = 0.01, while panel (c) depicts Toconfirmthatthediscrepanciesareindeedcausedby
the results from direct full 3D numerics using an spatial themeshresolution,whenwefurtherdecreasethelattice
spacingofh=0.2. Aclosecomparisonbetweenthe pan- spacing to h = 0.15 (results not shown here) and recal-
els suggests small discrepancies of the imaginary parts culate the spectrum in 3D, the branches at Im(λ) = 1
at Im(λ) = 1 and Im(λ) 0.2, and a spurious unstable indeedsplit into twonon-degeneratecurvesandthe spu-
mode with Re(λ) 0.08.≈These discrepancies (and their riousmodesatIm(λ) 0.2andRe(λ) 0.08begintoget
amendment as the≈mesh resolutionh decreases—see be- suppressed. This sugg≈ests that the di≈screpancies are in-
low) are likely due to the large spatial spacing h —that deedattributable to resolutioneffects and that the spec-
needstobe chosensuchthatthe 3Dcalculationsarestill trumcomputedwiththeeffectively1Dnumericalmethod
“manageable”— in which case the DSS radial profile is is indeed accurate.
under-resolved. It is crucial to note that the 3D calcula- The comparison between the stability spectrum pre-
tionsarefarmoreexpensivethantheirlowerdimensional dicted by the degenerate perturbation method (DPM)
(1D or 2D) counterparts considered above. andthe one obtained numerically nearthe linear limit is
6
2
1.5
1
λ
0.5
0
3.5 4 4.5 5 5.5
µ
FIG. 4: (Color online) Comparison near the linear limit of
thespectrumofthedegenerateperturbationmethodandthe
numerical result of the effectively 1D method involving the
sphericalharmonicdecomposition. Thickblue(dark)andor-
ange (light) lines depict the imaginary and real parts of the
numerical eigenvalues while the thinner cyan and green lines
depict the theoretical results of the degenerate perturbation
method. The eigenvalue depicted in green is the one respon-
sible for the destabilization of the DSS. The large red dots 2 t=t0
t=t0+T/6
ctioornrsesopfotnhdetpoertthuerbtwedoDfrSeqSudeenpciicetseedxitnraFcitge.d5f.roTmhethlaergoesrciflrlae-- r,t)|1 tt==tt00++TT//32
(
quencycorrespondstotheoscillatory modeof theDSSwhile ψ
|
the smaller one corresponds to the breathing mode of the
background density.
0 1 2 r 3 4
1.7
r01.6
shown in Fig. 4. As the figure shows, the DPM captures 1.5
thebehaviorofthespectrumclosetothelinearlimit. Let
0 20 40 60 80 100
us nowconfirmthe stability ofthe DSSinthis limit. We t
have performed a direct numericalintegrationof the full
FIG. 5: (Color online) Top set of panels: 2D density cuts
3DdynamicsoftheDSSforµ=4.5,whereitispredicted
(z = 0) from the full 3D numerical solution for µ = 4.5.
tobestable—cf.Figs.3and4. Thisevolutionisdepicted
TheinitialconditioncorrespondstothestationaryDSSatits
in the top set of panels of Fig. 5. Also, in this figure,
equilibrium radius. Second set of panels: Evolution of aDSS
we depict the stable oscillations of a perturbed DSS by shifted from its equilibrium position undergoing oscillatory
displacing its initial location away from the equilibrium motion over one oscillation period T. The white dashed cir-
radius. The evolutionfor such a DSS performing oscilla- cle denotesthestarting location of theDSSat thebeginning
tions is displayedin the secondsetof panels in Fig. 5. A of the period. Third panel: Radial profiles for the oscillat-
movie for this evolution is provided in the supplemental ing DSS over half a period. Bottom panel: Oscillations of
material. To more clearly visualize these oscillations, we the DSS radius vs. time. Blue dots represent the radius ex-
tracted from full 3D numerics while the solid red line is the
depict on the third panel of Fig. 5 the transverse den-
best (least-square) fit to a linear combination of oscillations
sity cuts at four different stages over half a period of
withfrequenciesw1 1.7581andw2 =2.0671(seereddotsin
the oscillation. These density cuts also reveal the exci- ≈
Fig. 4) corresponding, respectively, to the frequencies of the
tation of the background cloud which performs breath-
DSS oscillations and the breathing mode of the background
ing oscillations. From the density cuts we extracted the
cloud. ForamoviedepictingthestableoscillationsoftheDSS
(radial) location of the DSS along the evolution which please see thesupplemental material.
displays a beating-type quasiperiodic oscillation as de-
picted by the blue dots in the bottom panel of Fig. 5.
In fact, by (least-square) fitting a linear combination of
two harmonic oscillations (see red line) we extract the that they indeed belong to the stability spectrum of the
frequencies w 1.7581 and w = 2.0671. Upon closer stationary DSS state. It is important to stress that the
1 2
≈
inspection, these frequencies correspond, respectively, to observed DSS oscillations are supported by the stability
the frequency of the DSS oscillations and the frequency oftheDSSatitsequilibriumposition. Therefore,despite
of the breathing mode of the background cloud. These strong perturbations of the background that undergoes
frequencies are also depicted in Fig. 4 clearly showing breathingoscillations,theDSSrobustlypersistsandper-
7
forms stable oscillations about its equilibrium position.
Asthe spectrumindicates,seeFigs.3and4,forlarger
values of the chemicalpotential the DSS becomes unsta-
ble. To identify this instability, we have looked at the
eigenvectorsof the DPM. There arehigh degeneraciesin
the eigenvalues in this case, rendering harder the iden-
tification of the nature of the instability. Nevertheless,
the DPM is still helpful in that regard. For instance,
we know from the eigenvectors that the mode with the
largest slope decreasing at Im(λ) = 1, which causes in-
stabilities,isalinearcombinationof“up”states,andthe
“down”statesarenotinvolved. Therefore,theinstability
is due to a bifurcation rather than a collision with “neg- FIG. 6: (Color online) The VL6 cage (a) and the VR6 cage
ative” energy modes; the latter may lead to oscillatory (b)thatemergeasaresultoftheinstabilityoftheDSSstate.
instabilities, as discussed, e.g., in Ref. [39]. The former The states are captured from the numerical integration of
may be associatedwith symmetry breaking features and the GPE (1) dynamics at µ = 5.8 and 8.2, respectively. In
particular, the VR6 cage seems to be robust and appears in
exhibits the instability via real eigenvalue pairs. Direct
thedynamicsforawiderangeofchemicalpotentialsfromthe
numerical integration of Eq. (1) for chemical potentials
onset of the instability around µ=6.2.
µ=5.8,8.2 and 10.6, shows that the first instability is a
bifurcationofasixvortexlines(VL6)cage,whiletherest
of the instabilities are dominated by a cubic six vortex
rings(VR6)cage. Thefirstinstabilitysuggestsabifurca-
tion of VL6 from the interactions of the DSS and a dark 1
solitonstatewiththreeplanenodesperpendiculartothe
(x,y)-plane, separatedby anglesof π/3;the dark soliton 0.95
state (DS3) can be written as
0.9
ψ ρmcos(mφ)e−ωr2/2, (18)
DSm linear
| i ∝
µ0.85
with m = 3 and a linear eigenenergy (m + 3)ω. A √
2 /
two-mode analysis [39, 43] using this state and the DSS rωc 0.8
state yields a prediction of the bifurcation as occurring
at µc = 4.78, which is in good agreement with the nu- 0.75
mericalresultµc =4.85 0.01ofFig.3. Asimilartreat- PDE
±
ment of the VR6 cage turns out to be more challenging 0.7 α=0.8769
due to the difficulty in designing the correspondingdark α=4/5
soliton state. Instead, we have looked at the DS4 state, 0.65
0 20 40 60 80
whichcanformaneightvortexline(VL8)cage. Thetwo- µ
mode analysis then predicts a bifurcation at µ = 5.97,
c
which is in fair agreement with the second instability at FIG.7: (Coloronline)Theequilibriumradiusofthespherical
µ =6.25 0.01fromFig.3. Thenumericalresults,how- shell dark soliton, scaled by √µ/ω, as a function of µ (solid
c
ever,show±that the dynamically resulting state is a VR6 redline). Notethatthenumericalvaluesreachtheasymptotic
cage,notaVL8cage. PresumablytheVL8cagecanalso valueof p4/5 (dashedblack line), as perthesecond particle
picture, when µ is large. The first particle picture slightly
cause instabilities of the DSS but is less robust than the
VR6 cage,possibly due to the intersectionsofthe vortex overestimates rc (dash-dotted blueline).
lines, which are absent in VR6. Further analysis with
higher values of m (not shown here) also suggests the
progressively lesser relevance of the DSm state in caus-
inginstabilities. Ontheotherhand,theVR6stateseems
tobeveryrobustanddominatestheinstabilityinawide slightly larger equilibrium r , while the particle picture
c
range of chemical potentials from the onset of this in- basedonenergy conservationofthe DSS as anadiabatic
stability around µ = 6.2. The VL6 and VR6 cages are invariantagreesverywellwiththenumericalresults. Itis
shown in Fig. 6. interestingtonotethattheasymptoticbehavioronlysets
Finally, we present the numerical results for the equi- inataboutµ=30,deepinsidetheTFregime. Although
libriumlocationr ofthe DSS asafunctionofthe chem- thiswouldbe arathercomputationallydemandingpara-
c
ical potential µ and compare the results with those of metric range when using 3D (or even 2D) methods, our
the two particle pictures in the large density limit (see effective 1D approachbased on a spherical harmonic de-
Sec.III).Aplotofrcω/√µ asa functionofµ is depicted composition of the angular dependencies, allows us to
in Fig. 7. Note that the first particle picture predicts a reach such large values of µ.
8
V. CONCLUSIONS AND FUTURE A&M University for access to their Ada cluster.
CHALLENGES
APPENDIX: 2D DECOMPOSITION USING
In this work, we have revisited the theme of spher-
CYLINDRICAL COORDINATES
ical shell dark solitons. We have employed analytical
methods to study these structures both in the vicinity
of the linear limit, using a bifurcation analysis, and in Thegreatestadvantageofthe1Dmethodinvolvingthe
the Thomas-Fermi limit, where they can be treated as decomposition into spherical harmonics is that one can
effective particles. For the numerical simulations, we computethespectrumforlargechemicalpotentialswith-
have used continuation methods, as well as a quasi-one- out excessive computational costs. The reason that typ-
dimensional radial computation method, decomposing ical computations for large chemical potentials become
the angulardependence ofthe perturbationsinspherical computationallyprohibitive,isthatinthislimittheden-
harmonics, to determine their stability. We have found sityislargeand,thus,the widthoftherelevantlocalized
that spherical shell dark solitons can, in fact, be dynam- structures(proportionalto1/õ)becomesmuchsmaller
icallyrobust,spectrallystablesolutionsofthe 3DGross- thanthedomainsize(proportionaltoõ);hence,alarge
Pitaevskii equation within a certain interval of chemical number of mesh points is necessary to resolve properly
potentials (and atom numbers); this occurs sufficiently these configurations. Nonetheless, it is worth mention-
close to the linear limit of relatively small chemical po- ing that, for the projection method in 1D, the form of
tentials. Their instabilities, emerging when the chemical states that one can study is also restricted by the spher-
potential increases, have been elucidated and their dy- icalsymmetry; therefore, it is alsorelevant to employ an
namical outcome, namely the breakup into vortex ring analogof this method in 2D, so as to encompass a wider
and vortex line cages, has been showcased. We believe class of states. An example concerns states with a topo-
that these results, may not only be used to partially ex- logical charge S along the z axis, i.e., stationary states
plainthetransientobservationofthesestructuresinear- of the form ψ = ψ0eiSφ. Note that in this framework,
lier experiments [32] and numerical simulations [33], but onecanstudystatesincluding—butnotlimitedto—the
may also pave the way for their consideration in future ground-state,planar, ring,or sphericalshell dark soliton
experiments. states,vortexlines, aswellas coplanaror parallelvortex
Our work suggests a number of interesting future re- ringsandhopfions[44], inboth isotropicandanisotropic
search directions, concerning, in particular, related co- traps.
herent structures, their physical relevance and proper- In 2D, and for cylindrical coordinates (ρ,φ,z), the
ties. For instance, this study suggests large scale sta- Laplacian is decomposed as follows:
bility computations for a wide range of pertinent steady ∆ f
stateswithradialsymmetry,includingoneswithdifferent ∇2f =∆Hf + ρG2 , (19)
numbers of radial nodes as well as a number of different
external potentials. Furthermore, it would be relevant where its components read:
to generalize and apply the numerical (and analytical)
1 ∂ ∂f ∂2f
techniquesusedhereintomulticomponentBECsettings. ∆ f = ρ + , (20)
H ρ∂ρ ∂ρ ∂z2
Forexample,onecanstudy“symbioticstates”composed (cid:18) (cid:19)
of sphericalshell dark-brightsolitons,in the spirit ofthe ∂2f
∆ f = . (21)
dark-bright 2D rings of Ref. [45]. On the other hand, G ∂φ2
states emerging robustly from the instabilities of these
spherically(orcylindrically)symmetricones,suchasthe The∆Goperatorhaseigenstates eimφ witheigenvalues
m2, i.e., { }
vortexring(andline)cagesareworthwhiletoinvestigate
−
further in their own right. Efforts along these directions ∆ eimφ = m2eimφ. (22)
G
are currently in progress and will be presented in future −
publications. A stationary state ψ (ρ,z), defined in the domain ρ
0
[0, ) z R, and bearing a topological charge S, sat∈-
∞ × ∈
isfies the following equation:
Acknowledgments
1 S2
∆ ψ + ψ +V(ρ,z)ψ + ψ 2ψ µψ =0, (23)
−2 H 0 2ρ2 0 0 | 0| 0− 0
W.W. acknowledges support from NSF (Grant
Nos. DMR-1151387 and DMR-1208046). P.G.K. grate- As in the 1D case, we construct the linear stability
fullyacknowledgesthesupportofNSF-DMS-1312856,as problemas follows. Let ψ0 a rotationallysymmetric sta-
well as from the US-AFOSR under grant FA950-12-1- tionary state, up to a topological charge S, along the z
0332, and the ERC under FP7, Marie Curie Actions, axis perturbed using the complete basis of eimφ :
{ }
People, International Research Staff Exchange Scheme
(IRSES-605096). R.C.G. gratefully acknowledges the ψ =eiSφ ψ + a (ρ,z,t)eimφ+b∗ (ρ,z,t)e−imφ .
0 m m
support of NSF-DMS-1309035. We thank the Texas " #
m
X(cid:2) (cid:3)
9
Substituting thisexpansionintotheGPEofEq.(1),and by computing each mode m independently and then
upon linearizing and matching the basis expansion on putting them together. In the results based on this
both sides, one can see that the different m modes are methodandpresentedinSec.IV,weusem=0,1,2,...,5.
mutuallyindependent. Anequivalentderivationasin1D It is worth commenting here that the matrix for the 2D
shows that λ are the eigenvalues of the matrix eigenvalue problem with a charge S is less symmetric
than matrices corresponding to the 1D projection and
M M the 2D casewith no charge. Therefore,the full 2Dprob-
M = 11 12 ,
lemwith chargeappearsto demandmore computational
M21 M22!
work. However, it should be noted that, in the 2D case,
theeigenvaluesofthefullmatrixM correspondingtothe
where
projections over the +m and m modes are related via
M = i 1∆ + (m+S)2 +V +2ψ 2 µ , an orthogonal transformation−and, thus, the two set of
11 − −2 H 2ρ2 | 0| − eigenvaluesarecomplexconjugatesofeachother. There-
(cid:16) (cid:17)
M = iψ2, fore, it is only necessary to compute the non-negative m
12 − 0 setand,hence, the chargedandthe unchargedstates ac-
M = iψ∗2,
21 0 tually have similar computational complexity.
M = i 1∆ + (m−S)2 +V +2ψ 2 µ .
22 −2 H 2ρ2 | 0| −
(cid:16) (cid:17)
As before, one can therefore compute the full spectrum
[1] C. J. Pethick and H. Smith, Bose-Einstein Condensa- Kevrekidis, Phys. Rev.Lett. 101 (2008) 130401.
tion inDiluteGases(CambridgeUniversityPress,Cam- [16] G.Theocharis,A.Weller,J.P.Ronzheimer,C.Gross,M.
bridge, 2008). K.Oberthaler,P.G.Kevrekidis,andD.J.Frantzeskakis,
[2] L.P.PitaevskiiandS.Stringari, Bose-Einstein Conden- Phys. Rev.A 81 (2010) 063604.
sation (Oxford University Press, Oxford, 2003). [17] C. Becker, S. Stellmer, P. Soltan-Panahi, S. D¨orscher,
[3] P. G. Kevrekidis, D. J. Frantzeskakis, and R. Carretero- M. Baumert, E.-M. Richter, J. Kronj¨ager, K. Bongs, K.
Gonz´alez (eds.), Emergent Nonlinear Phenomena in Sengstock, Nature Physics 4 (2008) 496–501
Bose-Einstein Condensates. Theory and Experiment [18] S. Stellmer, C. Becker, P. Soltan-Panahi, E.-M. Richter,
(Springer-Verlag, Berlin, 2008); R. Carretero-Gonz´alez, S. D¨orscher, M. Baumert, J. Kronj¨ager, K. Bongs, and
D. J. Frantzeskakis, and P. G. Kevrekidis, Nonlinearity K. Sengstock, Phys. Rev.Lett. 101 (2008) 120406.
21, R139 (2008). [19] B. P. Anderson, P. C. Haljan, C. A. Regal, D. L. Feder,
[4] P. G. Kevrekidis, D. J. Frantzeskakis, and R. Carretero- L.A.Collins,C.W.Clark,andE.A.Cornell,Phys.Rev.
Gonz´alez, The defocusing Nonlinear Schr¨odinger Equa- Lett. 86, 2926 (2001).
tion: From Dark Solitons to Vortices and Vortex Rings [20] P. Engels and C. Atherton, Phys. Rev. Lett. 99, 160405
(SIAM,Philadelphia, 2015). (2007).
[5] K. E. Strecker, G. B. Partridge, A. G. Truscott, and R. [21] I. Shomroni, E. Lahoud, S. Levy and J. Steinhauer, Na-
G. Hulet,Nature 417, 150 (2002). ture Phys.5, 193 (2009).
[6] L.Khaykovich,F.Schreck,G.Ferrari,T.Bourdel,J.Cu- [22] C. Becker, K. Sengstock, P. Schmelcher, R. Carretero-
bizolles, L. D. Carr, Y.Castin, and C. Salomon, Science Gonz´alez,andP.G.Kevrekidis,NewJ.Phys.15,113028
296, 1290 (2002). (2013).
[7] S.L.Cornish,S.T.Thompson,andC.E.Wieman,Phys. [23] S.Donadello,S.Serafini,M.Tylutki,L.P.Pitaevskii,F.
Rev.Lett. 96, 170401 (2006). Dalfovo, G. Lamporesi, and G. Ferrari, Phys. Rev.Lett.
[8] D.J. Frantzeskakis, J. Phys.A 43, 213001 (2010). 113, 065302 (2014).
[9] O.Morsch andM. Oberthaler, Rev.Mod. Phys.78, 179 [24] T.Yefsah,A.T.Sommer,M.J.H.Ku,L.W.Cheuk,W.
(2006) J. Ji, W. S. Bakr, M. W. Zwierlein, Nature 499 (2013)
[10] A. L. Fetter and A.A. Svidzinsky, J. Phys.: Cond. Mat. 426–430.
13, R135 (2001). [25] M.J.H.Ku,W.Ji,B.Mukherjee,E.Guardado-Sanchez,
[11] A.L. Fetter, Rev.Mod. Phys. 81, 647 (2009). L.W.Cheuk,T.Yefsah,andM.W.Zwierlein,Phys.Rev.
[12] S.Komineas,Eur.Phys.J.-Spec.Topics147133(2007). Lett. 113, 065301 (2014).
[13] J. Denschlag, J.E. Simsarian, D.L. Feder, C.W. Clark, [26] G. Theocharis, D. J. Frantzeskakis, P. G. Kevrekidis, B.
L.A. Collins, J. Cubizolles, L. Deng, E.W. Hagley, K. A. Malomed, and Yu. S. Kivshar, Phys. Rev. Lett. 90,
Helmerson,W.P.Reinhardt,S.L.Rolston,B.I.Schneider, 120403 (2003).
and W.D. Phillips, Science 287 (2000) 97–101. [27] L. D. Carr and C. W. Clark Phys. Rev. A 74, 043613
[14] S. Burger, K. Bongs, S. Dettmer, W. Ertmer, K. Seng- (2006).
stock,A.Sanpera,G.V.Shlyapnikov,andM.Lewenstein, [28] G. Herring, L. D. Carr, R. Carretero-Gonz´alez, P. G.
Phys.Rev.Lett. 83 (1999) 5198–5201. Kevrekidis, and D. J. Frantzeskakis, Phys. Rev. A 77,
[15] A. Weller, J.P. Ronzheimer, C. Gross, J. Esteve, M.K. 023625 (2008)
Oberthaler,D.J. Frantzeskakis,G. Theocharis, and P.G. [29] S.Middelkamp,P.G.Kevrekidis,D.J.Frantzeskakis,R.
10
Carretero-Gonz´alez, and P. Schmelcher, Physica D 240, [38] J. Li, D.-S. Wang, Z.-Y. Wu, Y.-M. Yu and W.-M. Liu,
1449 (2011). Phys. Rev.A 86, 023628 (2012).
[30] W. Wang, P. G. Kevrekidis, R. Carretero-Gonz´alez, D. [39] R. N. Bisset, Wenlong Wang, C. Ticknor, R. Carretero-
J. Frantzeskakis,T. J. Kaper, and M. Ma, Phys.Rev.A Gonz´alez, D. J. Frantzeskakis, L. A. Collins, and P. G.
92, 033611 (2015). Kevrekidis, Phys. Rev.A 92, 043601 (2015).
[31] G. Theocharis, P. Schmelcher, M. K. Oberthaler, P. G. [40] D.L.Feder,M.S.Pindzola,L.A.Collins,B.I.Schneider,
Kevrekidis, and D. J. Frantzeskakis, Phys. Rev. A 72, and C. W. Clark Phys. Rev.A 62, 053606 (2000).
023609 (2005). [41] D.S.MorganandT.J.Kaper,PhysicaD192,33(2004).
[32] N.S.Ginsberg,J.Brand,L.V.Hau,Phys.Rev.Lett.94, [42] A. M. Kamchatnov and S. V. Korneev, Phys. Lett. A
040403 (2005). 374, 4625 (2010).
[33] S. Komineas and J. Brand Phys. Rev. Lett. 95, 110401 [43] G.Theocharis,P.G.Kevrekidis,D.J.Frantzeskakis,and
(2005). P. Schmelcher, Phys.Rev. E74, 056608 (2006).
[34] H.Salman, J. Comp. Phys. 258, 185 (2014). [44] R. N. Bisset, W. Wang, C. Ticknor, R. Carretero-
[35] M. R. Hermann and J. A. Fleck Jr., Phys. Rev. A 38, Gonz´alez, D. J. Frantzeskakis, L. A. Collins, and P. G.
6000 (1988). Kevrekidis Phys.Rev.A 92, 063611 (2015).
[36] G. C. Corey, J. Chem. Phys. 97, 4115 (1992). [45] J. Stockhofe, P. G. Kevrekidis, D. J. Frantzeskakis, and
[37] R.Koll´arandR.L.Pego,Appl.Math.Res.Express,2012 P. Schmelcher, J. Phys. B 44, 191003 (2011).
(2012) 1–46.