Double-layer Bose-Einstein condensates: A quantum phase transition in the transverse direction, and reduction to two dimensions
Mateus C. P. dos Santos, Boris A. Malomed, Wesley B. Cardoso
aa r X i v : . [ c ond - m a t . qu a n t - g a s ] S e p Double-layer Bose-Einstein condensates: A quantum phase transition in the transversedirection, and reduction to two dimensions
Mateus C. P. dos Santos, Boris A. Malomed,
2, 3 and Wesley B. Cardoso ∗ Instituto de Física, Universidade Federal de Goiás 74.690-970, Goiânia, Goiás, Brazil Department of Physical Electronics, School of Electrical Engineering, Faculty of Engineering,Tel Aviv University, and Center for Light-Matter Interaction, Tel Aviv University, Tel Aviv 69978, Israel Instituto de Alta Investigación, Universidad de Tarapacá, Casilla 7D, Arica, Chile
We revisit the problem of the reduction of the three-dimensional (3D) dynamics of Bose-Einsteincondensates, under the action of strong confinement in one direction ( z ), to a 2D mean-field equation.We address this problem for the confining potential with a singular term, viz ., V z ( z ) = z + ζ / z ,with constant ζ . A quantum phase transition is induced by the latter term, between the ground state(GS) of the harmonic oscillator and the 3D condensate split in two parallel non-interacting layers,which is a manifestation of the “superselection” effect. A realization of the respective physical settingis proposed, making use of resonant coupling to an optical field, with the resonance detuning modu-lated along z . The reduction of the full 3D Gross-Pitaevskii equation (GPE) to the 2D nonpolynomialSchrödinger equation (NPSE) is based on the factorized ansatz , with the z -dependent multiplier repre-sented by an exact GS solution of the 1D Schrödinger equation with potential V ( z ) . For both repulsiveand attractive signs of the nonlinearity, the 2D NPSE produces GS and vortex states, that are virtuallyindistinguishable from the respective numerical solutions provided by full 3D GPE. In the case of theself-attraction, the threshold for the onset of the collapse, predicted by the 2D NPSE, is also virtuallyidentical to its counterpart obtained from the 3D equation. In the same case, stability and instabilityof vortices with topological charge S =
1, 2, and 3 are considered in detail. Thus, the procedure of thespatial-dimension reduction, 3D → I. INTRODUCTION
Bose-Einstein condensates (BECs) have become a ver-satile platform for realization of various phenomena,such as the production of bright [1–4] and dark [7]solitons, dark-bright complexes [8], vortices [9] andvortex-antivortex dipoles [10–13], persistent flows in thetoroidal geometry [14–16], skyrmions [17], emulation ofgauge fields [18] and spin-orbit coupling [19], quantumNewton’s cradles [20], Anderson localization of mat-ter waves [21, 22], rogue waves [23], quantum droplets(self-trapped states supported by beyond-mean-field in-teractions) [24–33], etc. Further details can be found inreviews, both earlier [34–44] and more recent ones [45–53].Lower-dimensional BECs, i.e., bosonic gases tightlyconfined in one or two transverse direction by a strongpotential, make it possible to study specific phase tran-sitions and collective excitations in quantum settings[54, 55]. In particular, studies of quasi-two-dimensional(quasi-2D) BEC with embedded 2D potentials havedrawn much interest [36, 56]. In this connection, ap-proximations which make it possible to reduce the un-derlying 3D Gross-Pitaevskii equation (GPE) to effec-tive 1D [57–71] and 2D [63, 64, 66, 67, 73–78] equa-tions have been elaborated. In particular, effective low-dimensional equations were developed in Refs. [72],[73], and [74]. In the former work, the adiabatic ap- ∗ [email protected] proximation was applied, using an appropriate analyt-ical expression for the local chemical potential, to elim-inate the transverse dimensions and derive an effective1D equation that governs the axial mean-field dynam-ics of a strongly elongated (cigar-shaped) BEC with re-pulsive interatomic interactions. On the other hand, inRef. [57] effective 1D and 2D time-dependent nonpoly-nomial nonlinear Schrödinger equations (NPSEs) werederived with the help of the variational approximation,which accounts for the structure of the condensate in thetransverse directions. The use of the effective equationswith lower dimensionality is quite relevant, as such sim-plified models may help one to gain deeper insight intothe underlying dynamics, as well as to reduce the costof the computational work. In the experiment, low-dimensional BEC can be created by means of the sametechnique which is used in the 3D setting, i.e., laser cool-ing of atoms in magnetic and/or optical trapping poten-tials [27, 79].In this work, we aim to elaborate a model of BECloaded in a 2D planar harmonic-potential (HO) trapapplied in the plane of ( x , y ) , combined with compet-ing HO trap ∼ z and singular repulsive potential ∼ z acting along the perpendicular axis. This schemeshapes the condensate into a double-pancake configu-ration parallel to the plane of z =
0, as shown in Fig1. Then, using the technique similar to that elaboratedin Ref. [57], we derive an effective 2D equation govern-ing the system’s planar dynamics. Note that, differentlyfrom the attractive potentials ∼ − r , considered inworks [33, 80–84] (see also a brief review in [85]), which FIG. 1. (a) The confining transverse potential, V z = ζ / z + z , from Eq. (3). (b) The 3D isosurface of the local density, | ψ ( x , y , z ) | = g = − λ = ζ = may absorb atoms through the mechanism of the quan-tum collapse (alias “fall onto the center” [86]), the totalnumber of atoms is maintained constant, in the case ofthe repulsive singular potential.The rest of the paper is organized as follows. In thenext section we present the model, including a possi-bility of its physical realization, and derive an effectiveNPSE for it. This section also contains exact analyti-cal results for bound states and a spectrum producedby the singular trapping potential, and for a quantumphase transition between the HO potential and the onewhich includes the singular term. In Sec. III we checkthe accuracy of the effective 2D model, by comparingits predictions with results produced by the full three-dimensional GPE. This is done for static and dynamicalstates alike, including the ground state (GS) and vortexmodes. Stability and instability of vortices with topo-logical charges S =
1, 2, and 3 is also addressed in Sec.III. The paper is concluded by Sec. IV.
II. THE EFFECTIVE TWO-DIMENSIONAL NPSE(NONPOLYNOMIAL SCHRÖDINGER EQUATION) ANDANALYTICAL CONSIDERATIONSA. Basic equations
We start by considering the 3D GPE for atomic BEC,written in the usual scaled form [87, 88]: i ∂ψ∂ t = − ∇ ψ + V ψ + π g | ψ | ψ , (1)where ψ = ψ ( x , y , z , t ) is the mean-field wave functionwith the integral norm set equal to 1, ZZZ | ψ ( x , y , z ) | dxdydz =
1, (2) V ( x , y , z ) is a trapping potential, and g = Na s / a z is thestrength of the two-body interatomic interaction, with N being the number of atoms, while a s and a z are, respec-tively, the s -wave scattering length of atomic collisions and the confinement length of the HO potential acting inthe direction perpendicular to the system’s plane, whichis adopted as the unit of length. Attractive and repul-sive interactions correspond, respectively, to g < g > viz. , V ( x , y , z ) = (cid:18) ζ z + z (cid:19) + U ( x , y ) . (3)The latter term is chosen as the isotropic HO, U ( x , y ) = λ (cid:16) x + y (cid:17) , (4)with strength λ , while term 2 z in the transverse po-tential represents the usual one-dimensional HO trap[87, 88], with the strength normalized with the help ofthe scaling invariance of Eq. (1).The shape of the transverse potential, as defined byEq. (3), is displayed in Fig. 1 (a). As concerns the sin-gular repulsive term in the transverse part of the poten-tial in Eq. (3), with scaled strength ζ , it may represent aspecifically designed physical setting. Indeed, the repul-sive action on cold atoms may be exerted by a nearly-resonant blue-detuned optical field with frequency ω ,close to frequency ω of atomic dipole oscillations (see,e.g., Refs. [89, 90] and references therein), the respectiveinteraction energy being proportional to ( ω − ω ) − .This dependence may be used to induce an effective re-pulsion potential, imposing spatial modulation on ω by means of the Zeeman effect [86] in a spatially inho-mogeneous dc magnetic field, or quadratic Stark - LoSurdo effect [86] in an electrostatic field, cf. Ref. [91],where an inhomogeneous field was used to design spa-tially modulated dipole-dipole repulsion in BEC. In thepresent context, the field should be shaped so as to make ω ( z ) = ω − Ω z , which leads to the singular termin Eq. (4) with ζ ∼ Ω . In particular, in the caseof the Zeeman effect, the necessary spatial maximum orminimum of the magnetic field at z = ζ relevant for the experimental realization is estimatedas ∼
50 ¯m . If the realization is based on the Zeemaneffect, the creation of nonuniform magnetic field withthe respective values of the variation length is possible(see, e.g., Ref. [92]).The use of the transverse potential given by Eq. (3)is relevant for several reasons. First, as shown in detail,below, it offers a relatively simple setting for the realiza-tion of a quantum phase transition, and also gives riseto specific spectrum of excitations, which may be effectsinteresting for the experimental realization. In addition,the potential splits the 3D space in two uncommunicat-ing subspaces, each one admitting its own states, whichthus suggests a possibility to create “quantum chimeras"under the action of this potential. (a) ( x ) x 00.030.060.090.12 -12 -6 0 6 12 (b) ( x ) x 00.020.040.060.08 -16 -8 0 8 16 (c) ( x ) x FIG. 2. Normalized density profile ρ ( x ) in the central cross section (drawn through y =
0) of the 2D GS for the repulsive BECwith nonlinearity strength g = g =
10 (b), and g =
100 (c), in the presence of the in-plane potential (4) with λ = ζ =
1. The corresponding density profiles for the full 3D GPE (1), 2D NPSE (29), 2D cubic NLSE (32), and TFA (35) are plotted bychains of yellow squares, black solid lines, red dashed lines, and blue dashed-dotted lines, respectively.
B. The transverse wave function: analytical results and thequantum phase transition
The separation of the strong transverse and weak pla-nar potentials in Eq. (3) suggests that the reduction 3D →
2D may be facilitated by the adoption of the factor-ized ansatz for the 3D wave function, ψ ( x , y , z , t ) = χ ( z , t ) Φ ( x , y , t ) . (5)First, the substitution of this ansatz in Eq. (1), keepingonly the strong transverse potential and time derivative,leads to the 1D linear Schrödinger equation, i ∂χ∂ t = ˆ H z χ ≡ − ∂ χ∂ z + (cid:18) ζ z + z (cid:19) χ . (6)It is easy to find an exact solution for the GS of this equa-tion, χ ( z , t ) = ( )( α + ) p Γ ( α + ) | z | α exp (cid:16) − i µ ( ) z t − z (cid:17) , (7) α = + q + ζ , (8) µ ( ) z = + q + ζ , (9)where Γ is the Gamma-function, and the constant coef-ficient is determined by the normalization condition, Z + ∞ − ∞ | χ ( z ) | dz =
1. (10)Note that, in particular, α = ζ =
1, and α = ζ =
3. It is worthy to note that the singular term, ∼ ζ / z , in the integral which produces the GS energy, Z + ∞ − ∞ χ ˆ H z χ dx = µ z , (11)with Hamiltonian ˆ H z defined by Eq. (6), converges if thewave function (7) is substituted in Eq. (11). There is another formal solution to Eq. (6), repre-sented by Eq. (7) with α replaced by˜ α = − q + ζ (12)and a nominally lower eigenvalue,˜ µ z = − q + ζ , (13)instead of the values given by Eqs. (8) and (9). How-ever, this solution is physically irrelevant, as the wavefunction is singular at z = ζ ≥ µ z < ζ = | z | replaced by z inEq. (7), produces the first excited state. Furthermore,the actual energy of the singular state, as given by theintegral expression (11), diverges for all ζ >
0, break-ing the equality of the integral expression to µ z , Thus,there is a strong discontinuity ( quantum phase transition )in the spectrum of bound states produced by Hamilto-nian H z in Eq. (6), as there is a jump of the GS, followingthe introduction of an arbitrarily small value of ζ .The transverse potential displayed in Fig. 1(a) seemsas a structure built of two symmetric potential wells sep-arated by a tall barrier. This configuration is used forthe consideration of spontaneous symmetry breaking ofoptical and matter waves in various nonlinear photonicand BEC settings [93]. However, the singular potentialbarrier ∼ z − is so strong that it splits the system intwo non-communicating half-spaces (an effect known as“superselection” [94]), therefore the consideration of Eq.(1) with potential (3) on the entire axis, − ∞ < z < + ∞ ,amounts to solving the same problem on the half-axis,0 ≤ z < ∞ . Indeed, solutions (7) with all values ζ > χ = d χ / dz = z =
0, therefore any two different solutions, built in-dependently at z > z <
0, may be matched atpoint z =
0. This fact implies that, in the one hand, thesymmetry-breaking phenomenology becomes trivial inthe present setting, but, on the other hand, it opens apossibility to construct complexes of completely differ-ent states, which would seem as “quantum chimeras",cf. Refs. [95, 96].In the same vein, it is interesting to compare the sin-gular potential (7) to one featuring a more general sin-gularity, viz ., V z = ζ | z | h + z , (14)with h >
0. It is easy to see that, at h <
2, the expansionof the wave function at | z | → χ ≈ const · e − i µ t (cid:20) + ζ ( − h ) ( − h ) | z | − h (cid:21) , (15)except for the case of h =
1, when Eq. (15) is replaced by χ ≈ const · e − i µ t h + ζ | z | ln ( | z | ) i , (16)The regular dependence of the expansions in Eqs. (15)and (16) on ζ means that the change ζ = → ζ > h > z → χ ≈ const · e − i µ t exp − p ζ ( h − ) | z | ( h − ) /2 ! , (17)which definitely implies that a strong phase transitiontakes place. Thus, the case of h = χ = const · | z | α − z + p + ζ ! exp (cid:16) − i µ ( ) z t − z (cid:17) ,(18) µ ( ) z = + q + ζ , (19)where α is the same as given by Eq. (8). Moreover, itis easy to find an exact full spectrum of eigenvalues forhigher-order excited states, with number n (in the limitof ζ →
0, the spectrum carries over into energy eigen-values of HO with odd numbers, n OH ≡ + n ): µ ( n ) z = ( + n ) + q + ζ , n =
0, 1, 2, ... . (20)Note that this series of eigenvalues is equidistant , similarto the OH spectrum.The formal eigenvalue given by Eq. (13) also gen-erates an infinite series of higher-order ones, ˜ µ ( n ) z = ( + n ) − p + ζ , which are formal counterparts ofOH eigenvalues corresponding to even states, but all therespective wave functions are singular, i.e., unphysical.Eigenvalues ˜ µ ( n ) z are counterparts of those of HO witheven numbers, n OH ≡ n .Lastly, it is also worthy to note that the 2D version ofEq. (6), i.e., i ∂χ ∂ t = − (cid:18) ∂ ∂ r + r ∂∂ r + r ∂ ∂θ (cid:19) χ + (cid:18) ζ z + z (cid:19) χ , (21)where ( r , θ ) are polar coordinates in the 2D plane, pro-duces exact solutions with the azimuthal quantum num-ber, alias vorticity, l =
0, 1, 2, ... ( l = χ ( l ) = const · r η exp (cid:16) − i µ ( l ) t + il θ − r (cid:17) , (22) η = q ζ + l , µ ( l ) = (cid:18)q ζ + l + (cid:19) . (23)Note that the spectrum of the 2D eigenvalues, given byEq. (23), is not equidistant, unlike its 1D counterpart(20). C. Derivation of the two-dimensional equation
The spatial-dimension reduction 3D →
2D proceedsvia the variational approach, which is based on the La-grangian density corresponding to Eq. (1) with potential(3): L = i (cid:18) ψ ∗ ∂ψ∂ t − ψ ∂ψ ∗ ∂ t (cid:19) − |∇ ψ | − (cid:18) ζ z + z (cid:19) | ψ | − U ( x , y ) | ψ | − π g | ψ | . (24)We use the factorized ansatz (5), in which all the timedependence is included in the planar wave function, Φ , and the transverse one is adopted in the form sug-gested by solution (7), except for the phase factor,exp (cid:16) − i µ ( ) z t (cid:17) : ψ = [ Γ ( α + )] − √ σ ! α + | z | α exp (cid:18) − z σ (cid:19) Φ ,(25)where α is defined by Eq. (8), and σ = σ ( x , y , t ) is in-troduced as a variational parameter accounting for pos-sible evolution of the transverse-confinement size. Sta-tionary states generated by ansatz (25) imply the double-pancake shape in the 3D space, with ψ ( z = ) =
0, anda pair of symmetric maxima of density | ψ ( x , y , z ) | at z = ( α /2 ) σ , (26)see Fig. 1 (b) as an illustration. Note also that ansatz (25) is defined so that the transverse factor is subject tothe unitary normalization [cf. Eq. (10)], hence it followsfrom Eq. (2) that the 2D wave function also has its normequal to 1: ZZ | Φ ( x , y , t ) | dxdy =
1. (27)Next, by inserting ansatz (25) in the Lagrangian den-sity (24), performing the integration in the transverse di-rection, and neglecting the spatial derivatives of σ (inthe adiabatic approximation, cf. Ref. [57]), one can de-rive the corresponding effective Lagrangian: L eff = i (cid:18) Φ ∗ ∂ Φ ∂ t − Φ ∂ Φ ∗ ∂ t (cid:19) − |∇ ⊥ Φ | − (cid:18) a + (cid:19) × (cid:18) σ + σ (cid:19) | Φ | − U ( x , y ) | Φ | − π α Γ ( α + ) Γ ( α + ) g | Φ | σ ,(28)where ∇ ⊥ is the gradient operator in Cartesian coor-dinates ( x , y ) . This expression gives rise to the Euler-Lagrange equations: i ∂ Φ ∂ t = − ∇ ⊥ Φ + U ( x , y ) Φ + (cid:18) α + (cid:19) × (cid:18) σ + σ (cid:19) Φ + π α − Γ ( α + ) Γ ( α + ) g | Φ | σ Φ . (29) σ − π α ( α + ) Γ ( α + ) Γ ( α + ) g | Φ | σ − =
0. (30)The 2D nonpolynomial Schrödinger equation (NPSE),Eq. (29), is the main result of the derivation. It deter-mines the density profile in the 2D plane, taking intoregard effects of the transverse BEC structure.First, it is natural to consider the low-density limit ofEq. (29), i.e., with g | Φ | ≪ σ close to 1, as it fol-lows from Eq. (30): σ − ≈ π Γ ( α + ) α + ( α + ) Γ ( α + ) g | Φ | . (31)The small difference of σ from 1, given by Eq. (31), pro-duces no contribution to NPSE (29) in the lowest ap-proximation, hence this equation amounts to the usualnonlinear Schrödinger equation (NLSE) with the cubicterm: i ∂ Φ ∂ t = − ∇ ⊥ Φ + U ( x , y ) Φ + ( α + ) Φ + π α − Γ ( α + ) Γ ( α + ) g | Φ | Φ . (32)In the opposite high-density limit with the repulsivesign of the nonlinearity, g >
0, it is natural to apply the Thomas-Fermi approximation (TFA) to the under-lying three-dimensional GPE (1), neglecting the kinetic-energy term (second derivatives) in it. In the absenceof kinetic energy, one has additional scale invariance,which makes it possible to fix ζ = z = p ζ z ′ , t = t ′ / ζ , ψ = ζ − ψ ′ , λ = p ζλ ′ , g = ζ g ′ .(33)Then, dropping the primes, the TFA may be naturallybased on the following ansatz : ψ ( x , y , z , t ) = √ (cid:18) π (cid:19) z exp (cid:16) − z − i µ t (cid:17) Φ TF ,(34)cf. Eq. (25), where an expression for the 2D wave func-tion is derived by substituting this ansatz in the 3D equa-tion (1), dropping the second derivatives in it, and inte-grating the resultant equation along coordinate z : | Φ TF | = √ √ π g [ µ − − U ( x , y )] , at µ > U ( x , y ) + µ ≤ U ( x , y ) + µ TF , can be ob-tained from the normalization condition, ZZ | Φ TF ( x , y ) | dxdy =
1, (36)as per Eq. (27).
III. NUMERICAL RESULTS
GS solutions of the 2D and 3D equations addressed inthis work were produced by means of the well-knownmethod of the imaginary-time evolution [97]. It wasrealized by means of a split-step scheme, based on theCrank-Nicolson algorithm, with space and time steps ∆ x = ∆ t = g > g < FIG. 3. Chemical potential µ of the GS (ground state), trappedin the in-plane potential (4) with λ = ζ =
1, versusthe self-repulsion strength, g , Displayed are numerical resultsobtained from the full 3D GPE, effective 2D NPSE, 2D cubicNLSE, and TFA, under normalization conditions (2), (27) and(36), respectively. Symbols and curves have the same meaningas in Fig. 2.
468 2 4 6 8 10 (a) (b)
FIG. 4. (a) Mean-squared axial length h z i and (b) chemicalpotential µ of the GS versus the scaled strength ζ of the trans-verse potential barrier. The GS is trapped in the in-plane po-tential (4), with λ = A. GS (ground-state) solutions
1. The repulsive nonlinearity
Here, we compare the 1D lumped density profile pro-duced by the full 3D equation (1) for the stationary GS, ρ ( x ) = Z Z | ψ ( x , y , z , t ) | dydz , (37) with its counterpart obtained from the 2D equations(29), (32), and (35), calculated as ρ ( x ) = Z + ∞ − ∞ | Φ ( x , y , t ) | dy . (38)We start the numerical analysis by considering the re-pulsive nonlinearity ( g >
0) in the presence of the 2D(in-plane) HO trapping potential. Note that the config-uration of the BEC can be defined as double-pancake-shaped if the transverse confinement is much tighterthan the in-plane potential, i.e., λ ≪ α − , see Eq. (26).In Fig. 2 we display in-plane density profiles of theGS in the repulsive condensate ( g >
0) under the ac-tion of potential (4), for three different values of g . Incomparison to the full 3D GPE, the 2D NPSE providesvirtually exact results, while the low-density limit andTFA, based on Eqs. (32) and (35), respectively, producevisible discrepancies. In particular, for the case of rela-tively weak nonlinearity, with g =
1, displayed in Fig. 2(a), at the central point ( x =
0) the error between thenumerically exact value of the density, obtained fromthe full 3D GPE equation, and the 2D approximationsis ≃ ≃ ≃
30% for the 2D NPSE, cubic2D NLSE, and TFA, respectively (naturally, TFA is irrel-evant in the case of weak nonlinearity). In the oppositecase of strong nonlinearity ( g = ≃ ≃ ≃
3% (in this case, the TFA is quite relevant, while thelow-density approximation is not). Note that the errorproduced by the 2D NPSE increases with the increase of g , remaining, nevertheless, fairly small. Similar to thesituation considered in Ref. [57], this happens because ansatz (25), used for the 3D →
2D reduction method, istaken as a solution of Eq. (7) with g =
0, thus gettingless accurate with the increase of g .Another way of evaluating the accuracy of the 2DNPSE is through the calculation of chemical potential µ ,setting ψ ( x , y , t ) = ψ ( x , y ) exp ( − i µ t ) and Φ ( x , y , t ) = Φ ( x , y ) exp ( − i µ t ) in Eqs. (29) and (32), while, as men-tioned above, µ TF is defined by normalization condi-tion (36). Figure 3 shows that the chemical potentialsobtained from the 3D GPE and 2D NPSE always stayvery close, while the low-density approximation andTFA produce discrepancies. Note that the positive slopeof the µ ( g ) dependence is tantamount to d µ / dN >
0, if g is kept constant, while norm N is not fixed byEq. (27) but is allowed to vary. In turn, the latter con-dition ( the anti-Vakhitov-Kolokolov criterion ) is necessaryfor stability of localized states under the action of self-repulsion [100]. In fact, for simple modes, such as GS,this criterion may be sufficient for the stability, which isconfirmed by direct simulations of their perturbed evo-lution (not shown here in detail).We also analyzed the efficiency of the 2D NPSE withrespect to variation of the transverse scaled strength ζ .To this end, we set λ = g =
1, and calculate the (a) ( x ) x 00.060.120.180.240.3 -8 -4 0 4 8 (b) ( x ) x 00.120.240.360.48 -6 -3 0 3 6 (c) ( x ) x FIG. 5. The cross section of 2D density profiles ρ ( x ) for the GS in the attractive condensate with g = − g = − g = − λ = ζ =
1. The results produced by the full 3DGPE, 2D NPSE, and 2D cubic NLSE are shown, severally, by yellow squares, black solid lines, and red dashed lines. mean squared axial length, h z i = ZZ Z α + | z | ( α + ) e − z / σ Γ ( α + ) σ α + | Φ ( x , y , t ) | dxdydz ,(39)based on the ansatz (25) for the 2D NPSE, which we com-pare to h z i = RRR z | ψ ( x , y , z ) | dxdydz , calculated asper the 3D GPE. This result is shown in Fig. 4(a), wherewe see that the effective equation maintains its accuracywith the increase of ζ . Corroborating this result, in Fig.4(b) we present the behavior of the chemical potential µ ( ζ ) , in which the linear dependence µ ( ζ ) , producedby the 2D NPSE, is virtually overlapped with that pro-duced by the 3D GPE.
2. The attractive nonlinearity
Next, we address the case of the self-attraction, g <
0. For this case, in Fig. 5 we plot cross-sections of the2D density profiles, as produced by the full 3D GPE, 2DNPSE, and cubic 2D NLSE (the TFA is irrelevant for theattractive nonlinearity). Similar to the case of g >
0, the2D NPSE predicts the profiles in a virtually exact form,while the cubic 2D NLSE produces a visible discrepancywith the increase of | g | .Furthermore, Eq. (29) with g < | g | > | g c | , as maybe expected in the 2D setting [101]. We have found thevalue of g c for the attractive BECs in the presence ofthe in-plane trapping potential (4), under normalizationconditions (2) and (27), respectively, with λ = g c = − g c = − B. Vortex modes under the action of the repulsive andattractive nonlinearity
Here we address vortex solutions produced by the full3D and effective 2D equations. To this end, the 3D wavefunction is looked for, in the cylindrical coordinates, as ψ ( r , θ , z , t ) = Ψ ( r , z , t ) exp ( iS θ ) , resulting in the follow-ing equation: i ∂ Ψ ∂ t = − (cid:20) ∂ Ψ ∂ r + ∂ Ψ ∂ z + r ∂ Ψ ∂ r (cid:21) + S r Ψ + (cid:18) ζ + z z + λ r (cid:19) Ψ + π g | Ψ | Ψ , (40)where S is integer vorticity. Similarly, substituting Φ ( r , θ , t ) = φ ( r , t ) exp ( iS θ ) in the 2D NPSE equation(29) leads to the radial equation i ∂φ∂ t = − (cid:20) ∂ φ∂ r + r ∂φ∂ r (cid:21) + (cid:18) α + (cid:19) (cid:18) σ + σ (cid:19) φ + λ r φ + S r φ + π α − Γ ( α + ) Γ ( α + ) g | φ | σ φ ,(41)which is combined with Eq. (30).In Fig. 6, we display examples of radial density pro-files of the vortex modes produced by Eqs. (40) and (41),with ρ ( r ) = R + ∞ − ∞ | Ψ ( r , z , t ) | dz and ρ ( r ) = | φ ( r , t ) | , re-spectively. The profiles are presented for both g < g >
0, and for three values of the vorticity, S =
1, 2, 3.The results clearly corroborate the accuracy of the 2DNPSE in describing the vortex states of the 3D GPE.The dimensional reduction method employed here,which is based on ansatz (25), can also be used to pro-duce lumped density profiles of the condensate in theaxial direction, by the integration in the ( x , y ) plane: ρ ( z ) = α + π Γ ( α + ) | z | α Z ∞ exp (cid:18) − z σ (cid:19) | φ ( r , t ) | σ α + rdr .(42)It is relevant to compare this approximate result to itscounterpart, i.e., the integrated density, provided by the (a) ( r ) r (b) ( r ) r 00.00080.00160.00240.00320 4 8 12 16 (c) ( r ) r FIG. 6. Radial density profiles ρ ( r ) for states with vorticity S trapped in the in-plane potential (4) with λ = ζ = S = g = − S = g =
1, and (c) S = g = full 3D GPE as ρ ( z ) = π Z ∞ | Ψ ( r , z , t ) | rdr . (43)The comparison, presented in Fig. 7(a) for typical ax-ial profiles with S =
2, again demonstrates that the 2DNPSE offers virtually exact results, i.e., it can be reliablyused for the full description of the zero-vorticity andvortex modes. In particular, it demonstrates that the ra-dius of the central “hole” in the trapped mode, inducedby the embedded vorticity, increases with the growth of S , as can be clearly seen in Fig. 7(b). This is a generalproperty of solitons with embedded vorticity, which ad-mits an analytical explanation [102]. On the other hand,the axial density profile is very weakly affected by thevalue of S , as shown by Fig. 7(c).Addressing the onset of the collapse in vortex stateswith S =
1, 2 and 3, confined by the in-plane trappingpotential (4) with λ = g c , for the vortices subject to normalizationconditions (2) and (27). These results were produced bythe full 3D GPE, as well as by means of the 2D NPSEand cubic NLSE, which are collected in Table I, wherewe have also included the result for the GS ( S = S . Note that g c strongly increases with the growth of S ,similar to what was observed in other models [103]. C. Stability of the vortex modes with S =
1, 2, 3
To conclude the analysis of the vortex solutions underthe action of the cubic self-attraction, we studied theirstability by means of direct simulations, starting frominputs perturbed by anisotropic deformations, whichmay readily initiate splitting of unstable vortices in non-linear models [103–106]. For this purpose, we used, first,2D NPSE (29) with the in-plane trapping potential (4).The anisotropically deformed initial condition with vor- S ( g c ) − GPE ( g c ) − NPSE ( g c ) cubic 2D − NLSE g c for the onset of the collapse, calculatedfor the ground ( S =
0) and vortex ( S =
1, 2 and 3) states in theself-attractive BEC ( g < λ = ζ =
1. As indicatedin the table, the critical values are calculated as per the full3D-GPE and 2D approximations. ticity S was taken as Φ ( x , y , t = ) = β ( x + iy ) S exp (cid:20) − λ (cid:18) x + y γ (cid:19)(cid:21) ,(44)where β is a constant determined by the normalizationcondition (27), and γ is the parameter of the anisotropicdeformation. Below, we set γ = S = g > − S =
1, ismuch smaller than the modulus of the respective valueat the collapse point, g c ( S = ) = − x and y axes ( eccentricity oscillations [107]), arecaused by the initially imposed anisotropy in Eq. (44).A typical example is shown in Fig. 8(a), where we plota stable vortex profile with S = g = − t =
300 after nine cycles of the quasi-periodic eccen-tricity oscillations. These results are in agreement withthose previously published for 2D NPSE [74] and alsofor cubic 2D NLSE [103–106].In the instability region, i.e., at g < − g = − ( a) ( r ) r 00.160.320.480.64-4 -2 0 2 4 (b) ( z ) z00.160.320.480.64-4 -2 0 2 4 (c) ( z ) z FIG. 7. (a) Lumped axial density profiles ρ ( z ) for the repulsive condensate, with λ = g =
1, and vorticity S =
2, producedby the full 3D GPE and 2D NPSE (yellow squares and black solid lines, respectively. (b) Radial density profiles, ρ ( r ) , for therepulsive condensate with g =
50, under the action of the in-plane trapping potential (4) with λ = ζ =
1, as produced bythe 2D NPSE with normalization (27), for different values of the vorticity: S =
0, 1, 2, and 3. They are depicted, respectively, bythe orange solid, red dashed, blue dotted, and black dashed-dotted lines. (c) Axial density profiles, ρ ( z ) , for the same states as inpanel (b). All the profiles are practically overlapping. -12 -6 0 6 12-12-6 0 6 12 (x 10 -3 )xy (a) 0 4 8 12 16 -12 -6 0 6 12-12-6 0 6 12 (x 10 -3 )xy (b) 0 4 8 12-12 -6 0 6 12-12-6 0 6 12 (x 10 -3 )xy (c) 0 5 10 15 20 25 30 -12 -6 0 6 12-12-6 0 6 12 (x 10 -3 )xy (d) 0 5 10 15 20 25 FIG. 8. Density profiles | Φ ( x , y , t ) | display the evolution ofthe vortex states, with S =
1, which were initially subjectedto the elliptic deformation, as per Eq. (44) with γ = λ = ζ = λ = t = g = −
1. Panels (b), (c) and (d) display snapshots, at t =
0, 60, and140, respectively, of the periodically splitting and recombiningvortex profile, for g = − fission-fusion cycles repeat quasi-periodically. For ex-ample, at g = − τ ff ≃ τ r ≃
80. Thus, the fission-fusion cycles and rotation arenot strictly synchronized, the ratio of the respective pe-riods being τ ff / τ r ≃ | g | ,eventually leading to the onset of the collapse. For in- -12 -6 0 6 12-12-6 0 6 12 (x 10 -3 )xy (a) 0 2 4 6 8 10 -12 -6 0 6 12-12-6 0 6 12 (x 10 -3 )xy (b) 0 2 4 6 8 10 12-12 -6 0 6 12-12-6 0 6 12 (x 10 -3 )xy (c) 0 2 4 6 8 10 12 -12 -6 0 6 12-12-6 0 6 12 (x 10 -3 )xy (d) 0 6 12 18 FIG. 9. The same as in Fig. 8, but for the vortex with S = t =
0, and50, and 100, respectively, for g = − t =
300 for g = − stance, at g = − S =
1, the evolution of the initialprofile (44) ends up with the collapse at t ≃ S = g = − t =
0, 50, 100.In this case, the fission-fusion period is τ ff ≃
31. It isworthy to note that the respective oscillatory motion ofpivots of the two unitary vortices proceeds along the x -0 -12 -6 0 6 12-12-6 0 6 12 (x 10 -3 )xy (a) 0 3 6 9 -12 -6 0 6 12-12-6 0 6 12 (x 10 -3 )xy (b) 0 2 4 6 8 10-12 -6 0 6 12-12-6 0 6 12 (x 10 -3 )xy (c) 0 2 4 6 8 10 -12 -6 0 6 12-12-6 0 6 12 (x 10 -3 )xy (d) 0 4 8 12 16 FIG. 10. The same as in Fig. 9, but for S = axis, keeping y =
0. Simultaneously, the position of themaximum local intensity of the solution rotates with aperiod of τ r ≃
62. Thus, in this cases, the ratio of theperiods is τ ff / τ r ≃ | Φ ( x , y , t ) | is plotted in the ( x , y ) plane for g = − t = S =
3, producing results similar to those re-ported above for S =
2. The triplets are unstable againstfission into a rotating set of three unitary vortices, whosepivots are aligned in the radial direction. For small val-ues of | g | , such as g = − S = S =
2, with thesame values of the fission-fusion and overall-rotationperiods, see Eq. (46). On the other hand, permanentsplitting of the initial triple vortex in a rotating set ofunitary vortices, gradually separating in the radial di-rection, takes place at larger values of | g | . An exampleof the set of three radially separating eddies is plotted inFig. 10(d) for g = − t = S = S =
1, 2, and 3, in Fig. 11 we plot the evolu-tion of ratio h r i / h x i of the spatially averaged 2D radialvariable, r = x + y , and the squared coordinate, h r , x i ≡ ZZ ( r , x ) | φ ( x , y , t ) | dxdy , (47)for S =
1, 2 and 3 (the definition of the average doesnot include a normalization factor, as it cancels in ratio h r i / h x i ). The computation was performed using nu-merical solutions of the effective radial equation (41) in2D Cartesian coordinates.In Fig. 11(a) one observes quasi-periodic oscillationsof h r i / h x i in an unstable state with S = g = − (cid:10) x (cid:11) , when the rotating unitaryvortex is crossing axis x =
0. On the other hand, inFigs. 11(b) and (c) unstable vortices with S ≥ | g | demonstrate practically harmonic oscillationsof the same ratio, (cid:10) R (cid:11) / (cid:10) x (cid:11) , which is explained by thefact that in this case the periods of the fission-fusion cy-cle and rotation are commensurable, see Eq. (46).Finally, as shown above, unstable double and triplevortices, with S = | g | irreversibly splitinto sets of two or three unitary vortices displaced alongthe radial direction, keeping to separate in this direction.Accordingly, the respective curves in Figs. 11(b) and11(c) exhibit quasi-harmonic oscillations with a grow-ing amplitude. Eventually, the separation halts underthe action of the HO trapping potential in Eq. (41). IV. CONCLUSION
The objective of this work is to produce additionalresults for the important problem of the reduction ofthe full 3D dynamics of BEC, loaded in an external po-tential, which imposes strong confinement in one di-rection ( z ), to an effective 2D form. Here, this generalproblem is considered for the specific form of the con-fining potential, given by Eq. (3), which is a combi-nation of the singular repulsive term ζ / z and usualHO (harmonic-oscillator) trap. The singular term splitsthe 3D condensate into a pair of parallel non-interacting“pancakes”, which is an example of the “superselection”phenomenon. A physical realization of this configura-tion is proposed, in terms of a resonant optical field,whose frequency is subjected to appropriate modula-tion in direction z , perpendicular to the “pancakes”. Thereduction of the underlying 3D GPE (Gross-Pitaevskiiequation) to the 2D NPSE (nonpolynomial Schrödinger1 (cid:2)(cid:3) (a) t (cid:8)(cid:9) (cid:10)(cid:11)(cid:12) (cid:13)(cid:14)(cid:15) (b) (cid:16) (cid:17)(cid:18)(cid:19) (cid:20)(cid:21)(cid:22) (cid:25)(cid:26) (cid:27)(cid:28)(cid:29) (cid:30)(cid:31) (c) ! FIG. 11. The ratio of spatially averaged squared coordinates, h r i / h x i [see Eq. (47)], versus time, illustrating the unstableevolution of initial vortex states with (a) S = g = − S = g = − − S =
3. The other parameters are λ = ζ =
1, and γ = equation) is provided by the factorized ansatz , mak-ing use of the fact that the z -dependent potential ad-mits an exact GS (ground-state) solution of the respec-tive Schrödinger equation. The full spectrum of energyeigenvalues in this transverse potential is found in anexact form too. The potential demonstrates a quantumphase transition between the GS (ground state) of HOin the case of ζ = ζ >
0. The resulting two-dimensional NPSE (nonpoly-nomial Schrödinger equation), with the repulsive or at-tractive nonlinearity, produces GS and vortex-state so-lutions, and the threshold for the onset of the collapse,which are virtually identical to their counterparts ob-tained from the numerical solution of the underlying3D GPE. On the other hand, the 2D NLSE with theusual cubic nonlinearity, as well as the TFA (Thomas-Fermi approximation), give rise to conspicuous discrep-ancies, in comparison to the full 3D solution. Thus,the results demonstrate high accuracy of the appropri-ately formulated spatial-dimensionality reduction. Thismethod may be applied to other settings as well.In particular, the existence of stable vorticity stateswith S = S = ACKNOWLEDGMENTS