Regular and chaotic behavior of collective atomic motion in two-component Bose-Einstein condensates
RRegular and chaotic behavior of collective atomic motion in two-componentBose-Einstein condensates
Wei-Can Syu, ∗ Da-Shin Lee, † and Chi-Yong Lin ‡ Department of Physics, National Dong-Hwa University, Hualien, Taiwan, R.O.C. (Dated: 16/06/2020)We theoretically study binary Bose-Einstein condensates trapped in a single-well harmonic po-tential to probe the dynamics of collective atomic motion. The idea is to choose tunable scatteringlengths through Feshbach resonances such that the ground-state wave function for two types ofthe condensates are spatially immiscible where one of the condensates, located at the center of thepotential trap, can be effectively treated as a potential barrier between bilateral condensates of thesecond type of atoms. In the case of small wave function overlap between bilateral condensates,one can parametrize their spatial part of the wave functions in the two-mode approximation to-gether with the time-dependent population imbalance z and the phase difference φ between twowave functions. The condensate in the middle can be approximated by a Gaussian wave functionwith the displacement of the condensate center ξ . As driven by the time-dependent displacement ofthe central condensate, we find the Josephson oscillations of the collective atomic motion betweenbilateral condensates as well as their anharmonic generalization of macroscopic self-trapping effects.In addition, with the increase in the wave function overlap of bilateral condensates by properlychoosing tunable atomic scattering lengths, the chaotic oscillations are found if the system departsfrom the state of a fixed point. The Melnikov approach with a homoclinic solution of the derived z, φ , and ξ equations can successfully justify the existence of chaos. All results are consistent withthe numerical solutions of the full time-dependent Gross-Pitaevskii equations. PACS numbers: 03.75.Lm, 03.75.Fi, 05.30.Jp, 05.45.Ac
I. INTRODUCTION
The first experimental observation of interferencefringes between two freely expanding Bose-Einstein con-densates (BECs) has launched the fascinating possibilityto probe new quantum phenomena on macroscopic scalesrelated to the superfluid nature of the Bose condensates[1]. A chief effect is the Josepshon analog, in which thephase differences between two trapped BECs in a double-well potential can generate Josephson-like current-phaseoscillations via collectively atomic tunneling through thebarrier [2, 3]. Nevertheless, the nonlinearity of tunnelingeffects produces novel anharmonic Josephson oscillationsand macroscopically quantum self-trapping (MQST) ef-fects. The potentially existing such interesting phenom-ena indeed triggers a wide variety of theoretic investiga-tions in the settings of two-component BECs systems inthe double-well potential [4–7], triple-well [8] and the sin-gle BEC system in the triple-well [9], four-well [10] andmultiple-well potentials [11, 12], to cite a few. The ef-fective Josephson dynamics can also be observed in somesystems without an external potential [13]. More impor-tantly, several experiments have successfully observed co-herence oscillations through quantum tunneling of atomsbetween the condensates trapped in a double-well poten-tial [14–16], whereas similar phenomena were also seenfrom attractive to repulsive atomic interactions by ex- ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] perimentally changing atomic scattering lengths throughFeshbach resonances [17]. In addition to the regular co-herent oscillations mentioned above, the studies on thepossibility of the appearance of chaotic motions, by ap-plying an externally time-dependent trap potential, arealso a fascinating task with findings that deserve fur-ther experimental justification [18–21]. Chaotic behav-ior is also studied in the dynamics of of three coupledBECs [22].In the present work, we consider the two-componentBECs system with atoms in two different hyperfine statesin a single-well harmonic trap potential. Such a two-component BEC system provides us an ideal platformfor the study of intriguing phenomena, for example, mim-icking quantum gravitational effects in Ref. [23]. UsingFeshbach resonances to experimentally tune the scatter-ing lengths of atoms allows us to construct the diagramin Fig. 1, showing the typical ground-state wave functionof the binary condensates [24–28]. The effective param-eter characterizing the miscibility or immiscible regimeof binary condensates is mainly determined by the valueof ∆ = g g − g , where g ij denotes the atomic inter-action between i and j atoms [29–31]. The parameterregime for ∆ > < a r X i v : . [ c ond - m a t . qu a n t - g a s ] J un metric bilateral condensates can be prepared by meansof a magnetic field gradient [32–34], allowing the pos-sibility to observe collectively coherence oscillations ofatoms between them. This is an extension of the worksof [2, 3], where the system of a single BEC in an externaldouble-well potential is studied. The collective oscilla-tions of atoms through quantum tunneling between thecondensates centered at each of the potential wells havebeen observed where the full dynamics of the system canbe further simplified, and turned into an analog pendu-lum’s dynamics in the so-called two modes approximation[2, 3, 35, 36]. In our setting, the dynamics of the systemcan also be approximated by analogous coupled pendu-lums dynamics using the variational approach as well astwo modes approximation. One can then reproduce reg-ular oscillations around the stable states studied in Refs.[2, 3] within some parameter regime by considering theeffects from time-averaged trajectories of the central con-densate. Additionally, the collective motion of atoms, ifdeparting from the unstable states, could show chaoticalbehavior, oscillating between bilateral condensates due tothe time-dependent driving force given by the dynamicsof the condensate centered in the middle [18–21]. Thescattering lengths in the parameter regime, which givethe chaotic dynamics of the system, will be useful as areference for experimentalists to observe the effects.To be concrete, consider a mixture of Rb | (cid:105) = | F =2 , m F = − (cid:105) and Rb | (cid:105) = | F = 1 , m F = − (cid:105) byignoring the small mass difference [29]. The scatteringlengths we choose are estimated to be a = 50 a and a = 85 a with a the Bohr radius. The existenceof the Feshbach resonance at 160G permits us to tunethe scattering length a of Rb atoms from 50 a to150 a to be specified in each case later [29]. Let usconsider two-component elongated BECs in quasi-one-dimensional (1D) geometries with the trap potential [37] ω x = 2 π × , ω y,z = 2 π × N = 500, and N = 1000 as an example. Theexperimental realization of the Josephson dynamics bytwo weakly linked BECS in a double-well potential in asingle BEC system is confirmed in Ref. [14]. Both Joseph-son oscillations and nonlinear trapping are observed bytuning the external potential. Such phenomena are alsoexplored in two-component BECs system in a double-wellpotential where apart from their own atomic scatteringprocesses the Rabi coupling between the atoms from eachcomponents of the BECs is introduced. Tuning the scat-tering lengths using the associated Feshbach resonancesand/or the Rabi coupling constant allow us to observe thetransition from Josephson oscillations to nonlinear self-trappings [38]. Our proposal also in the two-componentBECs system but in a single-well potential suggests anovel setup, although involving more rich dynamical as-pects of condensates, the only tunable parameters areatomic scattering lengths again through Feshbach reso-nances to achieve the miscible condensates that needsmore experimental endeavor as compared with [38].Our presentation is organized as follows. In Sec. II, FIG. 1: The ground-state wave function of binary BECs sys-tem in a cigar-shape trap potential is shown for two species Rb | , − (cid:105) state and Rb | , − (cid:105) state with atomic cou-pling constants in miscible/immiscible regimes respectively.The numbers of atoms, say N = 500 , N = 1000, and alsothe potential trap ω x = 2 π × , ω y,z = 2 π × < >
0) below (above)red dashed line determined by ∆ = 0, where ∆ = g g − g .The black solid line separates between symmetric and asym-metric condensate wave functions. we first introduce the model of the binary BECs system,and the corresponding time-dependent Gross-Pitaeviskii(GP) equations for each of condensate wave functions. Inthis two-component BECs system, we study the evolutionof the wave functions within a strong cigar-shape trappotential so that it can effectively treated as a quasi one-dimensional system. We further assume that the con-densates start being displaced from their ground statein the immiscible regime where one of the condensatesis distributed in the trap center while the other conden-sate surrounds the central condensate on its two sides.Then we propose the ansatz of the Gaussian wave func-tions, which can reduce the dynamics of the system intofew degrees of freedom, such as the center of the centralcondensate as well as the population imbalance and therelative phase difference between bilateral condensatesin the two-mode approximation. We later obtain theirequations of motion with the form of coupled oscillatorsequations. In Sec. III, we analyze the stationary-statesolutions and their stability property. In Sec. IV, we ex-amine the time evolution of the system around the sta-ble stationary state, showing the regular oscillatory mo-tion. In Sec. V, the study is focused on when the systemstarts from near unstable stationary-state solutions, ex-hibiting the chaotic behavior by numerical studies. TheMelnikov homoclinic method is adopted for analyticallystudying the existence of chaos. Concluding remarks arein Sec. VI. In Appendix, we provide more detailed deriva-tions or approximations to arrive at coupled oscillatorsequations from the time-dependent GP equations usingthe appropriate ansatz of Gaussian wave functions. II. VARIATIONAL APPROACH ANDANALOGOUS COUPLED PENDULUMSDYNAMICS
We consider the binary BECs for same atoms in twodifferent hyperfine states confined in a strong cigar-shapepotential with the size of the trap L x along the axialdirection,taken in the x direction, which is much largerthan L r along the radial direction. This system can beeffectively treated as the pseudo one-dimensional systemwith the Lagrangian described as [39, 40] L = (cid:88) j =1 , (cid:90) dx (cid:34) i (cid:126) (cid:18) ψ j ∂ψ ∗ j ∂t − ψ ∗ j ∂ψ j ∂t (cid:19) − (cid:32) (cid:126) m (cid:12)(cid:12)(cid:12)(cid:12) ∂ψ j ∂x (cid:12)(cid:12)(cid:12)(cid:12) + V ext | ψ j | + g jj | ψ j | (cid:33) (cid:35) + (cid:90) dx g | ψ | | ψ | , (1)where the external potential in the axial direction is givenby V ext ( x ) = m ω x x /
2. The mass of atoms is m and thecoupling constants between atoms are given in terms ofscattering lengths a ij as g ij = 2 (cid:126) a ij /mL r . The wavefunction ψ j is subject to the normalization (cid:90) d x | ψ j ( x ) | = N j (2)with the numbers of atoms N j in each hyperfine states j .The time-dependent GP equations for each componentof the BEC condensates are given by [41] i (cid:126) ∂ψ ∂t = (cid:18) − (cid:126) m ∂ ∂x + V ext + g | ψ | + g | ψ | (cid:19) ψ , (3) i (cid:126) ∂ψ ∂t = (cid:18) − (cid:126) m ∂ ∂x + V ext + g | ψ | + g | ψ | (cid:19) ψ . (4)The full dynamics of the condensates in a harmonic trappotential can be explored by solving the GP equationsnumerically with given initial wave functions. However,in order to extract their feature analytically, we will pro-pose an ansatz of the wave functions with several relevantparameters. All results from solving the equations ofmotion for the introduced parameters given by the varia-tional approach will be justified by comparing with thosegiven by numerically solving the GP equations. More-over, we will rescale the spatial or temporal variablesinto the dimensionless ones by setting t → t/ (cid:126) ω x and x → x/ (cid:112) (cid:126) /mω x . Henceforth, the dimensionless t and x will be adopted, and their full dimensional expressionswill be recovered, if necessary. FIG. 2: The typical condensate wave functions are plotted forthe stationary-state solutions of the GP equations for scatter-ing length a = 110 a , a = 50 a and a = 85 a that give g /g = 2 . g /g = 1 . < Experimentally, the values of parameters g ij , which arerelated to scattering lengths a ij , can be tuned by Fesh-bach resonances [24, 25, 29]. In fact, we might explorethe evolution of condensate wave functions, in which ψ and ψ are initially prepared in the immiscible regime,determined by the parameter ∆ = g g − g when∆ < ψ centered atthe potential trap center as [39, 40, 42] ψ ( x, t ) = (cid:18) N πw ( t ) (cid:19) / exp (cid:20) − ( x − ξ ( t )) w ( t ) + ixα ( t ) + ix β ( t ) (cid:21) , (5)where four time-dependent variational parameters arethe position of the center ξ ( t ), the width w ( t ), and α ( t )(slope) and β ( t ) ((curvature) / ) of the wave function.We also assume that the ground-state configuration of ψ ( x, t ) is of the two-mode form [2, 3, 35, 36] ψ ( x, t ) = ϕ L ( t ) ψ L ( x ) + ϕ R ( t ) ψ R ( x ) . (6)The involved Gaussian-like spatial functions ψ L ( x ) and ψ R ( x ) together with the initial wave function ψ ( x, t = 0)with width σ will be determined numerically by thestationary-state solutions of the GP equations (3) and(4) with the forms shown in Fig. 2, which is symmet-ric with respect to x = 0. ψ L ( x ) and ψ R ( x ) presum-ably serve as a good basis to express the spatial part ofthe wave functions of bilateral condensates with normal-ization conditions (cid:82) dx | ψ L,R | = 1, and also negligibleoverlap (cid:82) dx ψ L ψ R ≈ x = 0. Thus, the time-dependent part of the wavefunction can be parametrized as ϕ L ( t ) = (cid:112) N L ( t ) e iφ L ( t ) , (7) ϕ R ( t ) = (cid:112) N R ( t ) e iφ R ( t ) (8)with the phase φ L,R ( t ) and the amplitude (cid:112) N L,R ( t )given by the number of atoms in each sides of the con-densates. The total number of atoms ( j = 1) is N = N L + N R . To realize quantum coherence between left andright condensates, and also atomic number oscillations,it finds more convenient to introduce the population im-balance z ( t ) ≡ N L ( t ) − N R ( t ) N L ( t ) + N R ( t ) , (9)and relative phase φ ( t ) ≡ φ R ( t ) − φ L ( t ) . (10)When all atoms locate in the left(right) side of the centralcondensate, z = 1( − g with the wave functions ψ ( x, t )and ψ ( x, t ), which give an additional effect to influencethe dynamics of coherent oscillations due to the wave-function overlap between them. Also, notice that al-though the initial wave function of the central conden-sate is centered at x = 0, which is the minimum ofthe external trapped potential, the later evolution of thewave function will find its stationary state with the cen-ter located at x = ξ , moving around that position withnonzero kinetic energy. Thus, in the presence of ψ ( x, t )shown in Fig. 2 together with the external trapped po-tential, the resulting effective potential experienced bythe wave function ψ ( x, t ) in the GP equation in (3)can be seen in the shape of the double-well potential.However, the dynamics of the central condensate leadsto the effective potential time dependent, contributingits effects to quantum coherence oscillations between bi-lateral condensates. This will be the main purpose ofstudies in Secs. IV and V. This is an extension of thework in Refs. [2, 3] with additional degrees of freedomof the central condensate. When the initial conditionsof the central or bilateral condensates are chosen neartheir stable-state configurations, the dynamics of coher-ence oscillations modifies slightly the behavior found inRefs. [2, 3]. Nevertheless, as their initial conditions arechosen near the unstable-state configurations, the chaos occurs. Although the analogous coupled pendulums dy-namics, which is reduced from the full dynamics of thesystem and later gives essential information on the timeevolution of the condensate wave functions, may producerelatively different trajectories from the solutions of thetime-dependent GP equations, they are still very usefulto provide us analytic analysis on the chaos dynamics byadopting the Melnikov Homoclinic method.Substituting the ansatz of the wave function (5) and(6) into the Lagrangian (1), and carrying out the in-tegration over space, the effective action then becomesthe functional of time-dependent variables α , β , ξ , w , z ,and φ . Their time evolutions follow the Lagrange equa-tions of motion derived from this effective Lagrangian,discussed in detail in the Appendix. When the wave func-tion ψ ( x, t = 0) is slight away from the stationary state,the dynamics of the condensate just undergoes small am-plitude oscillations around that state. We will furtherassume that the displacement ξ and the wave functionwidth, deviated from σ and defined by w = σ + σ , aresmall as compared with σ , namely ξ (cid:28) σ and σ (cid:28) σ .Retaining terms up to linear order in ξ and σ of in-terest, we have found the Lagrange equations for the pa-rameters by the variational approach as˙ z + 2 (cid:18) k + κ σσ (cid:19) (cid:112) − z sin φ = 0 , (11)˙ φ − ∆ E − Λ z − (cid:18) k + κ σσ (cid:19) z √ − z cos φ + 2 (cid:112) N /N ηξ = 0 , (12)¨ ξ + ω ξ ξ − (cid:112) N /N ηz − εσ = 0 , (13)¨ σ + ω σ σ − εξ = 0 . (14)The quantities that appear in the equations above arethe integrals involving the condensate wave functions,namely k = − (cid:90) dx (cid:18) dψ L dx dψ R dx + x ψ L ψ R (cid:19) , (15) k = k − g N √ πσ (cid:90) dxψ L ψ R e − x /σ , (16) κ = g N √ πσ (cid:90) dxψ L ψ R (cid:18) − x σ (cid:19) e − x /σ , (17)∆ E = (cid:90) dx (cid:34) (cid:18) dψ L dx (cid:19) + x ψ L + g N ψ L (cid:35) − (cid:90) dx (cid:34) (cid:18) dψ R dx (cid:19) + x ψ R + g N ψ R (cid:35) , (18)Λ = g N (cid:18)(cid:90) dxψ L + (cid:90) dxψ R (cid:19) , (19) η = − g √ N N √ πσ (cid:90) dx ( ψ L − ψ R ) x e − x /σ , (20) ε = g N √ πσ (cid:90) ∞−∞ dx ( ψ L − ψ R ) (cid:18) xσ − x σ (cid:19) e − x /σ . (21) ω ξ = 1 + g N √ πσ (cid:90) dx ( ψ L + ψ R ) (cid:18) x σ − (cid:19) e − x /σ , (22) ω σ = 1 + 3 σ + g N √ πσ + g N √ πσ (cid:90) dx ( ψ L + ψ R ) (cid:18) − x σ + 4 x σ (cid:19) e − x /σ . (23)The expressions of k (16) and κ (17) are determined bythe wave function overlap between the left and right con-densates, and thus k ∼ κ (cid:28) N > N . Nevertheless, it will beseen that since the frequency of collectively atomic os-cillations ω Jξ ∝ ω x √ k with ω x = 2 π × t ∝ / ( ω x √ k ) isconstrained to be not larger than the typical lifetime ofthis type of the condensates of order ∼
10s [37], giving k > − where N cannot be too large. The variation ofthe Gaussian width induces a time-dependent modifica-tion to the parameter k + κσ/σ in (11) and (12). In thesingle BEC experiment [2, 3], the similar modification to k can be achieved by the time-dependent laser intensityor the magnetic field with the possibility to induce theShapiro effect [43], analog Kaptiza pendulum [44] andeven chaotic motions [18–21]. Here, since σ (cid:28) σ , thetime-dependent modification is relatively small comparedwith k to be ignored in this work. Also, for simplicity,we consider that the spatial part of the wave function ψ ( x, t ) is symmetric with respect to x = 0 shown inFig. 2, thus giving ∆ E = 0 in (12). The coupling of the population imbalance z to the relative phase φ in(12) is given by the strength Λ in (19), a tunable valueby changing the scattering lengths of atoms in bilateralcondensates.It is known that the dynamical equations of pair vari-ables ( z, φ ) in the single BEC system show an analogyto that of a nonrigid pendulum model, whose length de-pends on the angular momentum z [3]. Nevertheless,in our binary BECs system, the pair variables ( z, φ ) arealso coupled to the displacement of the wave function ψ with the coupling constant η in (20) depending on theinterspecies coupling constant g of atoms, and also thewave function overlap between the bilateral and centralcondensates. With the parameters shown in Fig. 1 andthe chosen scattering lengths obeying ∆ <
0, which leadto the spatial parts of the wave functions schematicallyshown in Fig. 1, η is of order one η ∼ .
0, and ε ∼ . σ (0) = ˙ σ (0) = 0, thetime-dependent σ is driven by ξ , giving σ ∼ εξ , whichin turn leads to the term εσ in the equation of ξ (13)of the order of ε ξ with ε (cid:28) η to be safely ignored.From the viewpoint of the wave function ψ , the pres-ence of ψ R and ψ L will make trap potential narrowerand affects slightly the natural vibration frequency in ξ and give a modified frequency ω ξ in (22). Since we justfocus on various types of coherent oscillations betweenbilateral condensates with the dynamical variables z, φ ,apart from their mutual coupling, they are also coupledto the displacement of the central condensate ξ . Ignor-ing the εσ term in (13), the small deviation from thewave-function width σ is decoupled. In fact, in this ap-proximation, once the time-dependence of ξ is found, theequation of σ in (14) can be solved, and then pluggingall solutions to the equations (86) in Appendix can find α and β . Also notice that taking into account the time-dependent σ by involving (14) certainly can improve theagreement with the full numerical results. Thus, afterthe further replacement ξ → (cid:112) N /N ξ , the set of thekey equations can be simplified as˙ z + 2 k (cid:112) − z sin φ = 0 , (24)˙ φ − Λ z − k z √ − z cos φ + ηξ = 0 , (25)¨ ξ + ω ξ ξ − ηz = 0 . (26)These are the main results of this section and we willterm them the coupled pendulums (CP) dynamics. Thecorresponding Hamiltonian then reads as H = Λ2 z − k (cid:112) − z cos φ + 12 ˙ ξ + ω ξ ξ − η z ξ . (27)In our binary BECs case, the above Hamiltonian alsoshows an analogy between the dynamics of coupled pen-dulums and the BEC dynamics as in Ref. [3]. III. EQUILIBRIUM SOLUTIONS ANDSTABILITY ANALYSIS
In the two-component BEC system, the relevant equa-tions of motion for describing quantum coherent oscil-lations between bilateral condensates in the presence ofthe central condensate are effectively described by the dy-namics of two-coupled pendulums. In this case, the sta-tionary states with the vanishing time-derivative termsnow obey 2 k (cid:113) − z sin φ = 0 , (28)Λ z + 2 k z (cid:112) − z cos φ − η ξ = 0 , (29) ω ξ ξ − η z = 0 . (30)As a result, the stationary relative phase is given by φ = 0 , or φ = ± π. (31)Also, Eq. (30) leads to a relation between the popula-tion imbalance z and the displacement of the centralcondensate ξ ξ = ω − ξ ηz . (32)Inserting the relation (32) into Eq. (29) together with φ = 0 or ± π , the population imbalance is found to be z = 0 , or z = ± (cid:115) − k Λ , (33) where Λ eff = Λ − ω − ξ η . The solution z (cid:54) = 0 existsas long as Λ eff / k ≥ eff / k ≤ − eff can be pos-itive or negative, and it is tunable by changing either g or g to vary k (16), Λ (19), and η (20) . Thereare four different types of stationary states, ( z , φ ) =(0 , − I, (0 , ± π ) − II, ( ± (cid:112) − k / Λ , − III and( ± (cid:112) − k / Λ , ± π ) − IV, summarized in Table I. Fur-thermore, the corresponding energy for each stationarystates can be computed through (27) using the relationof stationary-state solutions (32), where the results arelisted in Table II. In the case Λ eff / k >
1, the groundstate is in configuration I where the bilateral condensateshave the same populations z = 0, thus the central con-densate is centered at ξ = 0, consistent with Fig. 1. Asfor Λ eff / k < −
1, the ground state shows the populationimbalance for z (cid:54) = 0 and the wavefunction ψ is centeredat ξ (cid:54) = 0 due to (32) for a given z also seen in Fig. 1.We next study the stability of the system around theequilibrium solutions. To do so, we consider the smallperturbations around ξ , z , and φ defined by ξ ( t ) = ξ + δξ ( t ) , (34) z ( t ) = z + δz ( t ) , (35) φ ( t ) = φ + δφ ( t ) . (36)Substituting (34)–(36) into Eqs. (24)–(26), and using ξ = ω − ξ ηz , one obtains the linearized equations of mo-tion for δξ ( t ) and δz ( t ) as (cid:32) δ ¨ ξ ( t ) δ ¨ z ( t ) (cid:33) = − ω ξ η kη (cid:113) − z cos φ k cos φ (cid:34) ( − Λ + 2Λ z − ω − ξ η z ) (cid:112) − z − k cos φ (cid:35) (cid:32) δξ ( t ) δz ( t ) (cid:33) . (37) TABLE I: Stationary states
Label z φ ω ± I 0 0 (cid:104) ω ξ + ω Jξ ± (cid:113) ( ω ξ − ω Jξ ) + 8 kη (cid:105) / ± π (cid:104) ω ξ + ω Jξ ± (cid:113) ( ω ξ − ω Jξ ) − kη (cid:105) / ± (cid:112) − k / Λ (cid:104) ω ξ + ω Jξ ± (cid:113) ( ω ξ − ω Jξ ) + 8 kη (cid:112) − z (cid:105) / ± (cid:112) − k / Λ ± π (cid:104) ω ξ + ω Jξ ± (cid:113) ( ω ξ − ω Jξ ) − kη (cid:112) − z (cid:105) / TABLE II: Mean field extreme for different Λ eff / k parame-ters stable unstable energyΛ eff / k > E IV > E I − < Λ eff / k < E II > E I Λ eff / k < − E II > E III
Notice that the linear term in δφ vanishes where φ = 0or ± π for the stationary-state configurations is evaluatedin obtaining the linearized equations of motion. Con-sidering the oscillation motion of δξ ( t ) and δz ( t ) with δξ ( t ) , δz ( t ) ∝ exp( iωt ), we obtain the eigenfrequencies ω ± =12 (cid:34) ω ξ + ω Jξ ± (cid:114)(cid:16) ω ξ − ω Jξ (cid:17) + 8 kη (cid:113) − z cos φ (cid:35) . (38)In addition to the nature frequency ω ξ for the centralcondensate, there exists an oscillating frequency ω Jξ = ω J + 2 k cos φ ω − ξ η z (cid:112) − z (39)with ω J = 2 k cos φ (cid:34) k cos φ + Λ(1 − z ) (cid:112) − z (cid:35) (40)which manifests the nature frequency for quantum os-cillations between two bilateral condensates ω J with theeffects from the dynamics of the central condensate [2].However, the true eigenfrequencies are the mixture ofthese two fundamental frequencies ω ξ , ω Jξ as they cou-ple. With the parameters in the immiscible regime, typi-cally for small k where the wave-function overlap betweenbilateral condensates is small, this gives ω ξ (cid:29) ω Jξ . Inthe two-component BECs system, the existing oscilla-tion frequency ω ξ of the central condensate will drive thedynamics of the pair variables ( z, φ ) into the oscillatorymotions with frequencies ω − and ω + where ω + (cid:29) ω − .Later we will manipulate the initial conditions of z (0) , φ (0) as well as ξ (0) , ˙ ξ (0) to effectively switch offthe fast varying mode of frequency ω + , giving the resultspresumably to those in Ref. [2]. Then, the presence ofthe fast varying mode for more general initial conditionswill also be studied to see its effects on quantum coher-ent phenomena ( z, φ ). We will learn that, in the case ofrelatively small k ∼ − with very small wavefunctionoverlap, the time average over the evolution of the fastvarying mode is adopted to find the deviation from theresults in Ref. [2]. Those trajectories are regular motionmoving around the stable states of the system listed inTable II. However, we also probe the parameter regime with relatively large wave-function overlap with relativelylarge value k ∼ − for amplifying the influence of thefast varying mode that can even drive the dynamics of z and φ to the chaotic behavior, if the system starts fromthe unstable state also seen in Table II. In our system,the dynamics of quantum coherence oscillations ( z, φ ) isvery sensitive to the tunneling energy k . For even larger k , the proposed wave-function ansatz fails to describe theevolution of ( z, φ ) where the dynamics of the system canonly be explored by directly solving time-dependent GPequations, which is beyond the scope of the present work. IV. REGULAR MOTIONA. General solutions for small amplitudeoscillations
After studying the stationary-state configurations andexamining their dynamical stability, we step forward todiscuss the time evolution of the system when z and ξ undergo small amplitude oscillations around the station-ary states. Recalling the linearized equations of motion(37), the general solutions will be the superposition oftwo eigenfunctions, δξ ( t ) = A ++ e iω + t + A + − e − iω + t + A − + e iω − t + A −− e − iω − t , (41) δz ( t ) = B ++ e iω + t + B + − e − iω + t + B − + e iω − t + B −− e − iω − t . (42)Substituting Eqs. (41) and (42) into (37), one gets therelations among those coefficients A + , ± = P + B + , ± , (43) A − , ± = P − B − , ± , (44)where P ± depends on the stationary-state solutions ob-tained as P ± = − kη (cid:112) − z cos φ (cid:34) ω ξ − ω Jξ ± (cid:114)(cid:16) ω ξ − ω Jξ (cid:17) + 8 kη (cid:113) − z cos φ (cid:35) . (45)According to the experimental feasibility where an ap-plication of a magnetic gradient trap may cause the dis-placements of condensates [32–34], we then consider theinitial conditions such as δφ (0) = 0, and δ ˙ ξ (0) = 0. Thecorresponding solutions are δξ ( t ) = (cid:18) − P + P − P + − P − δz (0) + P + P + − P − δξ (0) (cid:19) cos ω + t + (cid:18) P + P − P + − P − δz (0) − P − P + − P − δξ (0) (cid:19) cos ω − t, (46) δz ( t ) = (cid:18) − P − P + − P − δz (0) + 1 P + − P − δξ (0) (cid:19) cos ω + t + (cid:18) P + P + − P − δz (0) − P + − P − δξ (0) (cid:19) cos ω − t. (47)In the case that ψ L and ψ R are very spatially separated,resulting in small tunneling energy k in (16), two fre-quencies will have very different values, ω ξ (cid:29) ω Jξ , since ω Jξ ∝ √ k in (39). Then the eigenfrequencies (38) andthe coefficients (45) can be further approximated respec-tively as ω + → ω ξ , ω − → ω Jξ (48)and P + → − ω ξ kη (cid:112) − z cos φ + O ( k ) , P − → ηω ξ + O ( k ) . (49)The general solutions for small k now become δξ ( t ) (cid:39) ω ξ (cid:20) − ω ξ (cid:0) η δz (0) − ω ξ δξ (0) (cid:1) cos ω ξ t + η (cid:18) ω ξ δz (0) + 2 kη (cid:113) − z cos φ δξ (0) (cid:19) cos ω Jξ t (cid:21) , (50) δz ( t ) (cid:39) ω ξ (cid:20) kη (cid:113) − z cos φ (cid:0) η δz (0) − ω ξ δξ (0) (cid:1) cos ω ξ t + ω ξ (cid:18) ω ξ δz (0) + 2 kη (cid:113) − z cos φ δξ (0) (cid:19) cos ω Jξ t (cid:21) . (51)In this binary BEC system, it is important to explorethe effects from the additional mode with frequency ω ξ onthe coherent oscillations between bilateral condensates.To do so, we first choose the initial displacement of thecentral condensate and the initial population imbalanceas − P − δz (0) + δξ (0) = 0 in (46) and (47) where thischoice of the initial conditions leads to the amplitudeof the rapidly oscillatory motion vanishing, thus leavingwith the slowly varying motion only. For the situations of small k , using (49) we obtain a linear relation, δξ (0) = ω − ξ η δz (0) , (52)as in the stationary-state solutions (32). In this case, δξ ( t ) and δz ( t ) undergo slow oscillations and the lin-earized solutions (50) and (51) of the system read as ξ s ( t ) (cid:39) ξ + ω − ξ ηδz (0) cos ω Jξ t , (53) z s ( t ) (cid:39) z + δz (0) cos ω Jξ t. (54)Notice that together with (32) and (52), ξ ( t ) = ω − ξ ηz ( t )and Eqs.(24)–(26) reduce to the same set of the equationsin the single BEC case in Ref. [2] by lettingΛ eff = Λ − η ω − ξ . (55)In addition, the solution (53) shows ˙ ξ s (cid:28) ξ s .Another interesting initial condition is δξ (0) = 0 and inthis case we can activate properly the rapidly oscillatorymotion with the solutions ξ ( t ) (cid:39) ξ s ( t ) + C ξ cos ω ξ t , (56) z ( t ) (cid:39) z s ( t ) + C z cos ω ξ t , (57)where C ξ = − ηω − ξ δz (0) , (58) C z = 2 ω − ξ kη (cid:113) − z cos φ δz (0) (59)and C z (cid:28) C ξ in the regime of small k . Later we willexplore the effects from the rapidly oscillatory mode onthe dynamics of ( z, φ ), in the system. However, thissystem involves two largely separate frequencies, namely ω ξ (cid:29) ω Jξ , for small k with small wave-function overlapbetween bilateral condensates. Although numerical tra-jectories of the CP dynamics is straightforward, we findmore convenient to represent the trajectories in Poincar´emaps (stroboscopic plots at every period, 2 π/ω ξ ). Thesetrajectories in Poincar´e maps can be analytically un-derstood by considering the time-averaged effects overthe evolution of the ω ξ mode in the timescale T where1 /ω ξ (cid:28) T (cid:28) /ω Jξ , with which we construct the corre-sponding effective potential in next section. B. The effective potential
In order to construct an effective potential for the ( z, φ )dynamics, we substitute (24) into (27) to obtain4 k (1 − z ) − ˙ z = (cid:32) Λ2 z + ˙ ξ ω ξ ξ − ηξ z − H (0) (cid:33) , (60)where H (0) is a constant value determined by the initialstate ( z (0) , φ (0) , ξ (0) , ˙ ξ (0)). To get an effective poten-tial for the small amplitude oscillations around z and ξ that includes the effects from the ω ξ mode of the general C , C , and C terms, we substitute (56) and (57) into(60) giving4 k (cid:2) − ( z s + C z cos ω ξ t ) (cid:3) − ( ˙ z s − ω ξ C z sin ω ξ t ) = (cid:20) Λ2 ( z s + C z cos ω ξ t ) + 12 ( ω ξ C ξ sin ω ξ t ) + ω ξ ηω − ξ z s + C ξ cos ω ξ t ) − η ( ηω − ξ z s + C ξ cos ω ξ t ) ( z s + C z cos ω ξ t ) − H (0) (cid:21) , (61)where, since ˙ ξ s (cid:28) ξ s , the term ˙ ξ s is ignored for keep-ing the leading-order effects. Although C z (cid:28) C ξ inthe small k approximation, we still retain the terms of C z , which are small, for seeing how they contribute tothe effective potential. Let us now introduce the time-averaged population imbalance z over the timescale T for1 /ω ξ (cid:28) T (cid:28) /ω Jξ as (cid:104) z s ( t ) (cid:105) ≡ T (cid:90) T dtz s ( t ) = ¯ z ( t ) , (62)where ¯ z ( t ) will be determined self-consistently from theresulting effective potential. It is now straightforward tofind the expression˙¯ z + V eff (¯ z ) = 4 k − H (0) , (63)from which the effective potential is obtained as V eff (¯ z ) = V (¯ z ) + δV (¯ z ) . (64)Among them, V (¯ z ) is the well-known effective potentialfor a single BEC in a double-well trap potential [2] V (¯ z ) = ¯ z (cid:18) k − H (0)Λ eff + Λ z (cid:19) , (65)and δV (¯ z ) is the corrections due to the high-frequency ω ξ mode with the terms of C z and C ξ , given by δV (¯ z )= 12 (cid:20) k + ω ξ − Λ (cid:18) H (0) − Λ eff z (cid:19) + Λ ¯ z (cid:21) C z + ω ξ (cid:18) Λ eff z − H (0) (cid:19) C ξ + η (cid:18) H (0) − Λ eff z (cid:19) C z C ξ − ηω ξ C z C ξ + 18 (cid:0) eff ω ξ + 5 η (cid:1) C z C ξ − η Λ8 C z C ξ + 3Λ C z + ω ξ C ξ . (66) Furthermore, since ξ ( t ) = ω − ξ ηz ( t ) holds true in thecase of C z = C ξ = 0, the whole dynamics of ( z, φ, ξ ) canreduce to that of ( z, φ ) with symmetry of Λ eff → − Λ eff ,and φ → − φ + π as in Ref. [2]. Nonetheless, since C z and C ξ are in general nonzero, the system then does notobey the above-mentioned symmetry. We will verify thisin the later numerical studies. C. Results and Discussions Λ eff / k > (Josephson oscillation, running phase and π -mode self trapping) Consider the numbers of atoms in this binary conden-sates with the atom numbers N = 200 , N = 1000.In Fig. 3, we have chosen scattering lengths a =97 . a , a = 50 a , a = 85 a . As stated previously,we prepare the initial states of the condensates, theGaussian-like spatial functions ψ L ( x ) and ψ R ( x ) togetherwith the initial wave function of the central condensate ψ ( x, t = 0) with width σ determined numerically fromfinding the stationary-state solutions of GP equations (3)and (4) with the forms shown in Fig. 2. The resultingwave functions can in turn give the values of parameters k = 0 . .
91, and η = 1 .
25 via (16), (19), and(20), which lead to Λ eff / k = 68 > φ = 0and z = ξ = 0 shown in Fig. 3(a). The initial con-ditions are chosen as δφ (0) = δ ˙ ξ (0) = 0 together withvarious choices of δz (0) and δξ (0), obeying δξ (0) = ηω − ξ δz (0), with which only the ω Jξ mode becomes rele-vant. Thus, by varying δz (0) and δξ (0) accordingly, theplot shows the regimes of Josephon-like oscillations for¯ z (0) = z + δz (0) < ¯ z c (0), and MQST with the run-ning phase for ¯ z (0) > ¯ z c (0) as in Ref. [2]. The critical¯ z in this case is ¯ z c (0) = 0 .
24, obtained from (69) below.Also, we see the consistency of the solutions (50) and(51) with those of the equations in (24)–(26). In the caseof C z = C ξ = 0, the effective potential (64) reduces tothat in Ref. [2] drawn in Fig. 3(c) of a double-well form.Thus, for ¯ z (0) < ¯ z c (0) the system has large enough ki-netic energy ˙¯ z to move forward and backward betweentwo potential minima around the stationary state z = 0and φ = 0 undergoing the Josephson oscillations. As¯ z (0) increases to ¯ z (0) = ¯ z c (0), the initial kinetic energyof the system drives ¯ z rolling down toward one of thepotential minima and then climbing up the potential hillto reach the state of ¯ z = 0 just with ˙¯ z = 0, obeying thecondition ˙¯ z c (0) + V eff (¯ z c (0)) = V eff (¯ z = 0) . (67)Moreover, when ¯ z (0) > ¯ z c (0), the relatively small kineticenergy constrains the system to move around one of thepotential minima with z (cid:54) = 0, leading to the MQST, also0 FIG. 3: The phase portrait of the population imbalance z andthe phase difference between the bilateral condensates φ aswell as the plots of the effective potential with the parameterΛ = 1 . , η = 1 . , k = 0 . eff / k = 68 . δz (0) = 0 . , δξ (0) = 0 .
31 [blue(gray)], δz (0) = 0 . , δξ (0) = 0 .
30 [red (dark gray)], and δz (0) = 0 . , δξ (0) = 0 .
29 [green (light gray)] where the ω ξ mode is effectively switched; (b) with the initial conditions δz (0) = 0 . , δξ (0) = 0 [blue (gray)], δz (0) = 0 . , δξ (0) = 0[red (dark gray)], and δz (0) = 0 . , δξ (0) = 0 [green (lightgray)] where the the ω ξ mode is effectively activated. Otherinitial conditions for both cases are δφ (0) = 0, ˙ ξ (0) = 0. Theplots (c) and (d) are drawn for the corresponding effectivepotential for (a) and (b), respectively. shown in Fig. 3(c).Considering the initial conditions mentioned above to-gether with the population imbalance ¯ z (0) = z + δz (0)and the displacement of the central condensate ξ (0) = ξ + δξ (0) obeying the relation ξ (0) = ηω − ξ ¯ z (0), H (0)becomes H (0) = Λ eff z (0) − k (cid:112) − ¯ z (0) cos φ (0) . (68)Then, using Eq. (24) at the initial time to replace ˙¯ z c (0)in (67), for a given Λ eff , one can find the value of ¯ z c (0)with the effective potential (60) by setting C z = C ξ = 0,which gives the known result in Ref. [2]. For Λ eff / k > FIG. 4: The evolution of the population imbalance z (a) andthe displacement of the central condensate (b) as a functionof time in the case of the Josephson oscillations are com-pared between the results of the numerical solutions of thetime-dependent GP equations (blue line), and the linearizedequations in Eqs. (24)–(26) (red dashed line) with the sameparameters described in the text and the initial conditions z (0) = − . , φ (0) = ξ (0) = ˙ ξ (0) = 0. if the initial relative phase 0 ≤ | φ (0) | ≤ π/
2, then¯ z c (0) = 4 k Λ (cid:20) Λ eff − k cos φ (0)+ (cid:112) Λ eff (Λ eff − k ) cos φ (0) + 4 k cos φ (0) (cid:21) , (69)and when π/ ≤ | φ (0) | ≤ π we have¯ z c (0) = 4 k Λ (cid:20) Λ eff − k cos φ (0) − (cid:112) Λ eff (Λ eff − k ) cos φ (0) + 4 k cos φ (0) (cid:21) . (70)The critical value of ¯ z c (0) as a function of φ (0) drawsa boundary between the Josephson oscillations and theMQST mode to be discussed in the later numerical stud-ies in Fig. 9. As for a negative Λ eff , the solutions of ¯ z c (0)are given due to the following symmetry Λ eff → − Λ eff and φ → − φ + π . The value of ¯ z c (0) found in Fig. 3(a)and 3(c) can also been achieved from (69) and (70) with φ (0) = 0 for the positive value of ¯ z c (0).In Fig. 3(b), we present the evolutions of z ( t ) and φ ( t ),which include the effects from the ω ξ mode with the sameset parameters discussed previously (see figure captionfor details). To make a comparison with the evolution inFig. 3(a) of the slowly varying mode, it finds more conve-nient to represent the trajectories in Poincar´e maps (stro-boscopic plots at every period T = 2 π/ω ξ ) by collectingthe results from the section of ξ (0) = 0 and ˙ ξ (0) > z (0) < ¯ z c (0)) and therunning phase MQST dynamics ¯ z (0) > ¯ z c (0) are shownwith a transition between them at the shifted criticalvalue ¯ z c (0) = 0 . z c (0) canbe obtained from the conserved quantity H (0) = Λ eff z (0) − k (cid:112) − (¯ z (0) + C z ) cos φ (0)+ (Λ eff + ηω − ξ )2 C z + Λ eff ¯ z (0) C z + ω ξ C ξ − η C z C ξ , (71) where we have substituted the results of (53) and (54)into (27), and evaluated the Hamiltonian at an initialtime. Again, substituting (24) at an initial time into(67) and with the help of effective potential (64), we canobtain the critical value of ¯ z c (0) through the expressionas followsΛ eff = 1¯ z c (0) + 4 C z ¯ z c (0) − C z (cid:40) k (cid:112) − (¯ z c (0) + C z ) cos ( φ (0)) + ηC z C ξ − η ω − ξ C z + (cid:34) (cid:32) k (cid:112) − (¯ z (0) + C z ) cos ( φ (0)) + η C z C ξ − η ω − ξ C z (cid:33) + 16 k ¯ z c (0) (¯ z c (0) + 4 C z ¯ z c (0) − C z )[1 − C z (2¯ z c (0) + C z ) − (1 − (¯ z c (0) + C z ) ) cos ( φ (0))] (cid:35) / (cid:41) , (72)where C z and C ξ are given explicitly in terms of z and δz c (0) via (59) and ¯ z c (0) = z + δz c (0) in (54). For afixed Λ eff , the shifted value of ¯ z c (0) obtained from (72)is consistent with the numerical result in Fig. 3(b).In Fig. 4 we provide numerical comparison on the re-sults of the coupled pendulums equations (24)–(26) andthose by solving the full time-dependent GP equations(3) and (4). The initial wave functions for both conden-sates, as previously said, are prepared from the station-ary solutions of GP equations with the form similar toFig. 2. The wave function of the central condensate canbe approximated by (5) with ξ (0) = α (0) = β (0) = 0giving ˙ ξ (0) = 0. However, for the bilateral condensates,their initial wave functions are prepared for the functionsof x with the relative phase φ (0) = 0 and then tune thewave functions to give z (0) (cid:54) = 0. In Fig. 4, using the sameparameters and the initial conditions we find consistencybetween the solutions from analogous CP dynamics andthe time-dependent GP equations.In Fig. 5, the scattering lengths are slightly changed toprobe the dynamics of z, φ around the stationary state φ = π, z = 0 .
91 ( (cid:54) = 0) obtained from (33) with theparameters Λ = 1 . η = 1 . k = 0 . eff / k = 2 . >
2. We have the π -mode MQST oscil-lations in this case. The initial conditions of δz (0) arechosen being δz (0) (cid:28) z , where the effective potential constructed in(64) can safely be applied. The plots of Figs. 5(c) and5(d) of the respective effective potentials illustrate thesingle-well profile with the effects of turning on or off the ω ξ mode can satisfactorily interpret the dynamics of the population imbalance z shown in Fig. 5(a) and 5(b). Thisindicates that the MQST mode will persist for all rangeof z (0) with − < z (0) <
1. The effects from nonzero C z and C ξ seem changing slightly the turning points of¯ z , as ¯ z oscillates around z . For this Λ eff with φ (0) = π ,¯ z c (0) = 0 is found. This means that for − < z (0) < z c (0) asa function of φ (0) using the effective potential obtainedabove to be shown in Fig. 5. Similar check is done usingthe time-dependent GP equation as in Fig. 4 for the caseof MQST whereas in Fig. 6, large discrepancy is founddue to the fact that the two modes approach in (2) maynot provide a relatively good basis to parametrize thespatial parts of the wave functions of bilateral conden-sates. < Λ eff / k < (Josephson oscillation and π -mode selftrapping) In the interval 1 < Λ eff / k <
2, there exist two solu-tions of ¯ z c (0), if π/ ≤ | φ (0) | ≤ π ,¯ z c (0) = 4 k Λ (cid:20) Λ eff − k cos φ (0) ± (cid:112) Λ eff (Λ eff − k ) cos φ (0) + 4 k cos φ (0) (cid:21) . (73)2 FIG. 5: The phase portrait with the parameters Λ = 1 . η = 1 . k = 0 .
017 (Λ eff / k = 2 . π mode trajecto-ries circle the stationary point z = 0 . , φ = ± π, ξ = 1 . δz (0) =0 . , δξ (0) = 0 .
20 [blue, (gray)], δz (0) = 0 . , δξ (0) = 0 . δz (0) = 0 . , δξ (0) = 0 .
065 [green(light gray)] where the ω ξ mode is effectively switched: (b)with the initial conditions δz (0) = 0 . , δξ (0) = 0 [blue(gray)], δz (0) = 0 . , δξ (0) = 0 [red (dark gray)], and δz (0) = 0 . , δξ (0) = 0 [green (light gray)] where the the ω ξ mode is effectively activated. Other initial conditions forboth cases are δφ (0) = 0, δ ˙ ξ (0) = 0. In bottom panels, theplots (c) and (d) are drawn for the corresponding effectivepotential for (a), and (b) respectively. In Figs. 7 and 8, the parameter Λ eff / k is tuned to thevalue between 1 < Λ eff / k <
2. On the one hand, theJosephson oscillations are seen for φ (0) = 0, and however,the running phase MQST does not exist in this case inFig. 7 where there is no such a ¯ z c (0) found with the single-well potential, shown in (72) and (73). On the otherhand, the transition between the Josephson oscillationsand the MQST can occur in φ (0) = π with the double-well potential, giving the critical value of ¯ z c (0) either by(73) for the positive ¯ z c (0) with the ω Jξ mode only, or by(72) involving the contributions from the ω ξ mode. < Λ eff / k < (Josephson oscillation) As for 0 < Λ eff / k <
1, since there does not exist thesolution of ¯ z c (0) due to the fact that the correspondingeffective potential for both φ (0) = 0 and φ (0) = π showsa single well like the case in Fig. 7, we find the Josephsonoscillations for both relative phase difference cases.To summarize the discussion given above, we show the FIG. 6: The evolution of the population imbalance z (a) andthe displacement of the center condensate (b) as a functionof time in the case of the MQST are compared between theresults of the numerical solutions of the time-dependent GPequations (blue line), and the linearized equations in Eqs.(24), (25), and (26) (red dashed line) with the parameters a = 87 a , a = 50 a , a = 85 a (Λ = 2 . , k = 0 . , η =1 .
4) and the initial conditions are z (0) = 0 . , φ (0) = π, ξ (0) =˙ ξ (0) = 0. phase portraits in Fig. 9 from solving the double pen-dulum equations (24)–(26). The initial conditions wechoose in the plots of the top row are again z (0) = z + δz (0) , φ (0) = φ and ξ (0) = ξ + δξ (0) , ˙ ξ (0) = 0,where δξ (0) = ω − ξ ηδz (0) so that the mode of ω ξ iseffectively switched off. The bottom row figures cor-respond to the initial condition of δξ (0) = 0 insteadwith the solutions that the mode of ω ξ is effectivelyswitched on for comparison. For top row results, inFig. 9(e), there exists a transition between the Josephsonoscillations and the running phase MQST at φ (0) = 0for Λ eff / k >
2, and in Fig. 9(d), the transition oc-curs between the MQST and Josephson oscillations at φ (0) = ± π instead apart from the Josephson oscilla-tions at φ (0) = 0 for 1 < Λ eff / k <
2, and further inFig. 9(c) the Josephson oscillations at φ (0) = 0 , ± π re-spectively occur for 0 < Λ eff / k <
1. For Λ eff / k < C ξ = C z = 0, all results for Λ eff / k > φ → − φ + π . Nevertheless, including the effectsfrom the ω ξ mode, such symmetry mentioned above doesnot exist. V. CHAOTIC DYNAMICS
In previous sections, we mainly focus on variousregimes of the regular motions. In particular, when twofrequencies ω + and ω − have large difference in magni-tude, the dynamics of the rapidly varying mode is aver-aged out giving a mean effect to the mode of slow os-cillations. However, as the difference in magnitude be-tween two frequencies becomes small with a relativelylarger overlap between bilateral condensates, the systemexhibits chaotic phenomena with the scattering lengthsin the regimes marked in Fig. 10. From experimental per-spectives, this can be achieved by increasing the scatter-3 FIG. 7: In top panels, the phase portrait with the parametersof Λ = 1 . η = 1 . η = 0 .
005 (Λ eff / k = 1 .
07) is plottedfor the Josephson oscillation trajectories around the station-ary state: z = 0 φ = 0 , ξ = 0 (a) with the initial conditions, δz (0) = 0 . , δξ (0) = 0 .
38 [blue (gray)], δz (0) = 0 . , δξ (0) =0 .
25 [red (dark gray)], and δz (0) = 0 . , δξ (0) = 0 .
12 [green(light gray)] where the ω ξ mode is effectively switched: (b)with the initial conditions δz (0) = 0 . , δξ (0) = 0 [blue(gray)], δz (0) = 0 . , δξ (0) = 0 [red (dark gray)], and δz (0) = 0 . , δξ (0) = 0 [green (light gray)] where the ω ξ modeis effectively activated. Other initial conditions for both casesare δφ (0) = 0, δ ˙ ξ (0) = 0. In bottom panels, the plots (c) and(d) are drawn for the corresponding effective potential for (a),and (b) respectively. ing length a of the bilateral condensates or decreasingthe inter-species scattering length a , thus resulting inrelatively large tunneling energy k . In Fig. 11, we con-sider the scattering lengths giving Λ eff / k > z = ξ = 0, φ = ± π , a hyperbolic fixed point, whichis a unstable point by changing z , realized from the ef-fective potential in Fig. 3 and also in Table II, but is astable point by changing the phase φ instead. We solvethe time-dependent GP equations using the initial con-dition ψ ( x,
0) with ξ (0) = 0 and ˙ ξ (0) = 0. The bilateralcondensate wave function ψ ( x,
0) is constructed by spa-tial wave functions ψ L and ψ R with the initial populationimbalance z (0) = 0 .
09 and the initial relative phase dif-ference φ (0) = π . The dynamics of z is shown to run be-tween the Josephson oscillations and the running phaseMQST in Fig. 11(b). With the same initial conditionsfor z (0) , φ (0) , ξ (0) and ˙ ξ (0), the numerical results of cou-pled pendulum equations (24), (25), and (26) also showthe similar dynamics from which solutions we can com-pute the so-called Lyapunov exponents for each of thedynamic variables. A chaotic motion can be understood FIG. 8: In top panels, the phase portrait with parametersas in Fig.7 is plotted for the π mode trajectories around thestationary state: z = 0 . , φ = π, ξ = 0 .
47 (a) with theinitial conditions, δz (0) = 0 . , δξ (0) = 0 . δz (0) = 0 . , δξ (0) = 0 .
127 [red (dark gray)], and δz (0) =0 . , δξ (0) = 0 .
064 [green (light gray)] where the ω ξ modeis effectively switched: (b) with the initial conditions δz (0) =0 . , δξ (0) = 0 [blue (gray)], δz (0) = 0 . , δξ (0) = 0 [red(dark gray)], and δz (0) = 0 . , δξ (0) = 0 [green (light gray)]where the the ω ξ mode is effectively activated. Other initialconditions for both cases are δφ (0) = 0, δ ˙ ξ (0) = 0. In bottompanels, the plots (c) and (d) are drawn for the correspondingeffective potential for (a), and (b) respectively. from the fact that, for two initial nearby initial variables q ( II ) (0) → q ( I ) (0), their difference grows exponentially intime t , obeying | q ( II ) ( t ) − q ( I ) ( t ) | = | q ( II ) (0) − q ( I ) (0) | e λt , where the rate λ is the Lyapunov exponent obtained from λ = lim t →∞ lim q ( II ) (0) → q ( I ) (0) t ln (cid:12)(cid:12)(cid:12)(cid:12) q ( II ) ( t ) − q ( I ) ( t ) q ( II ) (0) − q ( I ) (0) (cid:12)(cid:12)(cid:12)(cid:12) . (74)There are four Lyapunov exponents λ z , λ φ , λ ξ , and λ ˙ ξ shown in Fig. 11(c). The chaotic dynamics can be rec-ognized when at least one of them is a nonzero positivevalue at large time, and it is found λ z > ξ ( t ) in (56) and (58)into (27) where again the full dynamics of the systemreduces to that of ( z, φ ) only. Then the correspondingHamiltonian splits into two parts H = H z + H I ( t ) , (75)4 FIG. 9: The phase portraits, choosing different η with fixed parameters of Λ = 1 . , k = 0 .
005 and ω ξ = 1 by solving (24), (25),and (26). The top row figures correspond to the initial conditions: z (0) = z + δz (0) , φ (0) = φ , and ξ (0) = ξ + δξ (0) , ˙ ξ (0) = 0,where δξ (0) = ηω ξ δz (0) away from the various stationary states by changing δz (0) with the solutions that the mode of ω ξ iseffectively switched off. The bottom row figures correspond to the initial condition of δξ (0) = 0 in stead with the solutionsthat the mode of ω ξ is effectively switched on for comparison. The parameters we choose are as follows : (a) Λ eff / k = − . η = 1 . eff / k = − . η = 1 . eff / k = 0, η = 1 . eff / k = 1 . η = 1 . eff / k = 2 . η = 1 . eff / k > z = ξ = 0, φ = ± π seenin Fig. 9(d) (1 < Λ eff / k <
2) and Fig. 9(e) (Λ eff / k > eff / k < − z = ξ = 0, φ = 0 instead, which is also a hyperbolic fixed point in this rangeof Λ eff / k in Fig. 9(a) (Λ eff / k < −
2) and in Fig. 9(b)( − < Λ eff / k < − where H z is the Hamiltonian to account for the dynamics FIG. 11: (a) The phase portrait for the parameters: N =500 , N = 1000 , a = 137 a , a = 50 a , and a = 85 a that result the effective parameters Λ = 7 . , η = 2 .
45, and k = 0 . z (0) = 0 . , φ (0) = π, ξ (0) = ˙ ξ (0) = 0. The red dashedline corresponds to the linearized equation Eqs. (24), (25),and (26) and blue line corresponds to the simulation resultsof real-time GP equations. (c) The corresponding Lyapunovexponents obtained from (74) are λ z , λ ˙ ξ , λ ξ , λ φ from top tobottom. of the slowly varying mode given effectively by H z = Λ eff z − k (cid:112) − z cos φ (76)5 FIG. 12: (a) The phase portrait for the parameters: N =500 , N = 1000 , a = 77 a , a = 50 a , and a = 77 a that result the effective parameters Λ = 4 . , η = 2 .
18, and k = 0 . z (0) (cid:39) , φ (0) (cid:39) π, ξ (0) = ˙ ξ (0) = 0. The red dashedline corresponds to the linearized equation Eqs. (24), (25),and (26) and blue line corresponds to the simulation resultsof real-time GP equations. (c) The corresponding Lyapunovexponents obtained from (74) are λ z , λ ˙ ξ , λ ξ , λ φ from top tobottom. In the case of relatively small difference between the val-ues of ω − and ω + , the dynamics of the fast varying mode ω + can be treated as the time dependent perturbationcontributing to the interaction Hamiltonian as H I ( t ) ≡ ∆ ¯ E ( t ) z ( t ) = η ω − ξ z (0) cos ( ω + t ) z ( t ) . (77)Nevertheless in Refs. [18–21], this time-dependent per-turbation is implemented by introducing an externaldriving force. Here, assuming H I (cid:28) H z , one can writethe solution of z as z = ζ + ˜ ζ where ζ is the solution ofthe unperturbed equation¨ ζ − (Λ eff H z − k ) ζ + Λ ζ = 0 , (78)and ˜ ζ is the first-order correction in z due to the pertur-bation g ( t ) satisfying the equation¨˜ ζ − (Λ eff H z − k )˜ ζ + 32 Λ eff ζ ˜ ζ = g ( t ) , (79) g ( t ) = ∆ ¯ E ( t ) H z + H I Λ eff ζ −
32 Λ eff ∆ ¯ E ( t ) ζ . (80)When Λ eff H z − k >
0, Eq. (78) presents a homoclinicsolution with the double-well potential shown in Fig. 3for Λ eff / k > < Λ eff / k < H z = 2 k . The homoclinic solution might start from aparticular initial condition, and the system will drive z rolling down toward one of the potential minima, andthen climbing up the potential hill of the state with z = 0,a fixed point, just with zero velocity. The solution ζ reads as ζ ( t ) = 2 (cid:113) (Λ eff − k ) / Λ sech (cid:16)(cid:112) Λ eff − k t + D (cid:17) , (81) D = arcsech (cid:32) z (0)2 (cid:112) (Λ eff − k ) / Λ (cid:33) . (82)This homoclinic solution is also a separatrix, which isa boundary between the Josephson oscillations and therunning phase MQST for Λ eff / k >
2, and also theJosephson oscillations and the π phase MQST for 1 < Λ eff / k <
2. One can then construct the Melnikov func-tion from the homogeneous solution of (79), which is˜ ζ = ˙ ζ , and the function g ( t ) as M (0) = (cid:90) ∞−∞ dt ˙ ζ ( t ) g ( t ) . (83)The existence of zero of the Melnikov function shows thechaos [18, 19]. Substituting Eq. (79) into (83) above turnsout to be M (0) = − πω + Λ eff (cid:112) Λ eff − k ∆ ¯ E (0) × (cid:20) H z √ Λ eff − k − (cid:112) Λ eff − k (cid:0) ν (cid:1)(cid:21) × sech (cid:16) πν (cid:17) sin ( D ν ) , (84)where ν = ω + / √ Λ eff − k . For the case of the fixed point z = 0, φ = π , the value of D → ∞ means that withsuch highly rapidly oscillations of the sine function, anyfinite ω + including the frequency given by the fast vary-ing mode will make the Melnikov function zero. Thus,the chaos occurs as shown in Fig. 11 for Λ eff / k > < Λ eff / k <
2, and is shownin Fig. 12. In this case, the dynamics of z runs betweenthe Josephson oscillations and the π mode MQST in-stead, where its behavior can also be analyzed using theMelnikov Homoclinic method discussed above. Moreover,the chaos may appear also in Λ eff / k < − z = 0 , φ = 0 (see Table II) in theregime of the scattering lengths shown in Fig. 10. Noticethat although the semiclassical (mean-field) GP equa-tion can reliably signal the onset of the chaotic motion,it might fail to provide the details of the dynamics on theshort timescales right after entering the chaotic regime,which can be dominated by quantum many-body effects.The exploration of chaos beyond the mean field approx-imation deserves further studies [45, 46]. VI. CONCLUSIONS
In summary, we have proposed a new setting with bi-nary BECs in a single-well trap potential to probe the6dynamics of collective atomic motion. In this setting, Rb atoms | , − (cid:105) and Rb atoms | , − (cid:105) are consid-ered with tunable scattering lengths via Feshbach reso-nances so that the ground-state wave function for twotypes of the condensates are spatially immiscible shownin Fig.1. As such, the condensate of atoms for one ofthe hyperfine states centered at the potential minimumcan be effectively treated as a potential barrier betweenbilateral condensates formed by atoms in the other hy-perfine state. In the case of small wave-function overlapof bilateral condensates, one can parametrize their spatialpart of the wave functions in the two-mode approxima-tion together with time-dependent population of atomsand the phase of each of the wave functions. Besides, thewave function of the condensate in the middle is approx-imated by an ansatz of the Gaussian wave function. Thefull system can be reduced to the dynamics of the imbal-ance population of atoms in bilateral condensates z , aswell as the relative phase difference φ between two wavefunctions together with the time-dependent displacementof the central condensate ξ . For small wavefunction over-lap of bilateral condensates shown in Fig. 10, all sorts ofthe regular trajectories, moving about the stable states,in Refs. [2, 3] can be reproduced. Moreover, the numer-ical results given by solving the equations of z, φ , and ξ are in close agreement with the solutions of the fulltime-dependent GP equations. Nevertheless, with an in-crease in wave-function overlap also shown in Fig. 10,we study the possibility of the appearance of the chaotic oscillations driven by the time-dependent displacementof the central condensate. The application of the Mel-nikov approach with the homoclinic solutions of the z, φ ,and ξ equations successfully predicts the existence of thechaos, which are further justified from solving the fulltime-dependent GP equations. All of the findings in thiswork deserve further experimental investigations usingadvanced techniques for manipulation of atomic conden-sates. Acknowledgments
This work was supported in part by the Ministry ofScience and Technology, Taiwan.
VII. APPENDIX
Section II has discussed a variational approach to thedynamics of binary BEC system. In this appendix weprovide more detailed derivations and approximations toarrive at the equations of the CP dynamics given in (24)–(26). Substituting the ansatz of the ground-state wavefunction (5) and (6) into the Lagrangian (1), and carryingout the integration over space, the effective Lagrangianthen becomes a functional of the time-dependent vari-ables α , β , ξ , w , z , and φ as L = L − N (cid:18) − z ˙ φ + ∆ Ez + Λ2 z − k (cid:112) − z cos φ (cid:19) − N (cid:40) ( ˙ αξ + ˙ βξ + ˙ β w ) + 12 (cid:20) w + 2 β w + ( α + 2 βξ ) (cid:21) + 14 (2 ξ + w ) + g N √ πw (cid:41) − g N N √ πw (cid:90) dx (cid:104) ( ψ L + ψ R ) + z ( ψ L − ψ R ) + 2 (cid:112) − z cos φψ L ψ R (cid:105) exp (cid:20) − ( x − ξ ) w (cid:21) , (85)where the effective parameters are tunneling energy k (15), difference of energies between the wells ∆ E (18)and self-interaction energy Λ (19). Using Lagrangianequations for the parameters α , β , ξ , w , z , and φ , wefirst obtain α = ˙ ξ − ξβ and β = ˙ w w . (86)Then, after inserting Eqs. (86) into (85), the equationsof motion for population imbalance and relative phase between bilateral condensates become˙ z + (cid:18) k − g N √ πw (cid:90) dxψ L ψ R e − ( x − ξ )2 w (cid:19) × (cid:112) − z sin φ = 0 , (87)˙ φ − ∆ E − Λ z − (cid:18) k − g N √ πw (cid:90) dxψ L ψ R e − ( x − ξ )2 w (cid:19) z √ − z cos φ − g N √ πw (cid:90) dx ( ψ L − ψ R ) e − ( x − ξ )2 w = 0 . (88)7For the center condensate, the equations of motion of ξ and w are obtained as¨ ξ + ξ + g N √ πw (cid:90) dx (cid:2) ( ψ L + ψ R ) + z ( ψ L − ψ R )+2 (cid:112) − z cos φψ L ψ R (cid:105) ( x − ξ ) e − ( x − ξ )2 w = 0 , (89)¨ w + w − w − g N √ πw + g N √ πw (cid:90) dx (cid:104) ( ψ L + ψ R )+ z ( ψ L − ψ R ) + 2 (cid:112) − z cos φψ L ψ R (cid:105) × (cid:20) x − ξ ) w − (cid:21) e − ( x − ξ )2 w = 0 . (90)It can be understood that the presence of bilateral con-densates contribute a time-dependent deformation forthe central condensate by coupling to population imbal- ance z ( t ) in the integrand of Eqs. (89) and (90). Weintroduce Gaussian functions (see Fig. 2) ψ L ( x ) = (cid:18) πλ (cid:19) / e − ( x + ζ )22 λ , (91) ψ R ( x ) = (cid:18) πλ (cid:19) / e − ( x − ζ )22 λ , (92)and substitute them into Eqs. (87)–(90). In our cases,the displacement ξ and the width defined as w = σ + σ with σ determined initially by solving time-independentGP equation for finding the ground state solution and σ driven by ξ , satisfy the conditions ξ (cid:28) (cid:113) λ + σ and σ (cid:28) (cid:113) λ + σ . (93)In the case of λ ∼ σ , we have ξ (cid:28) (cid:112) λ + σ ∼ σ and σ (cid:28) (cid:112) λ + σ ∼ σ . Back to Eqs. (87)–(90), it is allowedto expand the equations in terms of small σ/σ and ξ/σ as˙ z + (cid:34) k − g N √ πσ (cid:18) − σσ + O (cid:18) σ σ (cid:19)(cid:19) (cid:90) dxψ L ψ R e − ( x − ξ )2 σ (1 − σ/σ ) (cid:35) (cid:112) − z sin φ = 0 , (94)where the integral above can be further expanded as (cid:20) − σσ + O (cid:18) σ σ (cid:19)(cid:21) (cid:90) dxψ L ψ R e (cid:20) − (cid:16) xσ (cid:17) +2 (cid:18) xξσ (cid:19) + O (cid:18) ξ σ (cid:19)(cid:21)(cid:18) − σσ + O (cid:18) σ σ (cid:19)(cid:19) = (cid:90) dxψ L ψ R e − x /σ (cid:20) − σσ (cid:18) − x σ (cid:19)(cid:21) + (cid:90) dxψ L ψ R e − x /σ (cid:18) xξσ (cid:19) + O (cid:18) ξ σ , σ σ (cid:19) . (95)The term of order ξ/σ vanishes due to the odd function in the integrand. Therefore Eq. (87) can be further simplifiedby keeping the terms of order σ/σ as˙ z + (cid:20) k − g N √ πσ (cid:90) dxψ L ψ R e − x /σ + 2 g N √ πσ (cid:90) dxψ L ψ R e − x /σ (cid:18) − x σ (cid:19) (cid:18) σσ (cid:19)(cid:21) (cid:112) − z sin φ = 0 , (96)and reduces to (11) with the definition of the constants in (16) and (17). Following the same procedure to approximateEq. (88), we have˙ φ − ∆ E − Λ z − (cid:20) k − g N √ πσ (cid:90) dxψ L ψ R e − x /σ + 2 g N √ πσ (cid:90) dxψ L ψ R e − x /σ (cid:18) − x σ (cid:19) (cid:18) σσ (cid:19)(cid:21) × z √ − z cos φ − g N √ πσ (cid:90) dx ( ψ L − ψ R ) x e − x σ ξ = 0 , (97)which leads to (12).8Now we turn to linearize Eq. (89), which is¨ ξ + ξ + g N √ πσ (cid:20) − σσ + O (cid:18) σ σ (cid:19)(cid:21) (cid:90) dx (cid:104) ( ψ L + ψ R ) + z ( ψ L − ψ R ) + 2 (cid:112) − z cos φψ L ψ R (cid:105) × ( x − ξ ) e − x /σ (cid:20) − x σ + 2 (cid:18) xσ (cid:19) (cid:18) ξσ (cid:19) + 2 (cid:18) x σ (cid:19) (cid:18) σσ (cid:19) + O (cid:18) ξ σ , σ σ (cid:19)(cid:21) = 0 . (98)Considering the vanishing of the integral due to the odd function in the integrand, we conclude¨ ξ + (cid:20) g N √ πσ (cid:90) dx ( ψ L + ψ R ) (cid:18) x σ − (cid:19) e − x /σ (cid:21) ξ + (cid:20) g N √ πσ (cid:90) dx ( ψ L − ψ R ) xe − x /σ (cid:21) z + (cid:20) g N √ πσ (cid:90) dx ( ψ L − ψ R ) (cid:18) x σ − xσ (cid:19) e − x /σ (cid:21) σ = 0 , (99)that gives (13). Finally, we can linearize Eq. (90) as¨ σ + (cid:20) σ + g N √ πσ + g N √ πσ (cid:90) ∞−∞ dx ( ψ L + ψ R ) (cid:18) − x σ + 4 x σ (cid:19) e − x /σ (cid:21) σ − (cid:20) g N √ πσ (cid:90) ∞−∞ dx ( ψ L − ψ R ) (cid:18) xσ − x σ (cid:19) e − x /σ (cid:21) ξ = 0 , (100)giving (14). [1] M. R. Andrews, C. G. Townsend, H.-J. Miesner, D. S.Durfee, D. M. Kurn, and W. Ketterle, Science , 637(1997).[2] A. Smerzi, S. Fantoni, S. Giovanazzi, and S. R. Shenoy,Phys. Rev. Lett. , 4950 (1997).[3] S. Raghavan, A. Smerzi, S. Fantoni, and S. R. Shenoy,Phys. Rev. A , 620 (1999).[4] B. Juli´a-D´ıaz, M. Guilleumas, M. Lewenstein, A. Polls,and A. Sanpera, Phys. Rev. A , 023616 (2009).[5] B. Sun and M. S. Pindzola, Phys. Rev. A , 033616(2009).[6] I. I. Satija, R. Balakrishnan, P. Naudus, J. Heward,M. Edwards, and C. W. Clark, Phys. Rev. A , 033616(2009).[7] M. Mel´e-Messeguer, B. Juli´a-D´ıaz, M. Guilleumas,A. Polls, and A. Sanpera, New Journal of Physics ,033012 (2011).[8] Andrea Richaud and Vittorio Penna, New Journal ofPhysics , 105008 (2018).[9] B. Liu, L.-B. Fu, S.-P. Yang, and J. Liu, Phys. Rev. A , 033601 (2007).[10] S. De Liberato and C. J. Foot, Phys. Rev. A , 035602(2006).[11] F. S. Cataliotti, S. Burger, C. Fort, P. Maddaloni, F. Mi-nardi, A. Trombettoni, A. Smerzi, and M. Inguscio, Sci-ence , 843 (2001).[12] M. Nigro, P. Capuzzi, H. M. Cataldo, and D. M. Jezek,Phys. Rev. A , 013626 (2018).[13] M. Abad, M. Guilleumas, R. Mayol, M. Pi, and D. M.Jezek, Phys. Rev. A , 035601 (2011).[14] M. Albiez, R. Gati, J. F¨olling, S. Hunsmann, M. Cris- tiani, and M. K. Oberthaler, Phys. Rev. Lett. , 010402(2005).[15] L. J. LeBlanc, A. B. Bardon, J. McKeever, M. H. T.Extavour, D. Jervis, J. H. Thywissen, F. Piazza, andA. Smerzi, Phys. Rev. Lett. , 025302 (2011).[16] M. Pigneur, T. Berrada, M. Bonneau, T. Schumm,E. Demler, and J. Schmiedmayer, Phys. Rev. Lett. ,173601 (2018).[17] G. Spagnolli, G. Semeghini, L. Masi, G. Ferioli,A. Trenkwalder, S. Coop, M. Landini, L. Pezz`e, G. Mod-ugno, M. Inguscio, A. Smerzi, and M. Fattori, Phys.Rev. Lett. , 230403 (2017).[18] F. K. Abdullaev and R. A. Kraenkel, Phys. Rev. A ,023613 (2000).[19] C. Lee, W. Hai, L. Shi, X. Zhu, and K. Gao, Phys. Rev.A , 053604 (2001).[20] H. Jiang, H. Susanto, T. M. Benson, and K. A. Cliffe,Phys. Rev. A , 013828 (2014).[21] J. Tomkoviˇc, W. Muessel, H. Strobel, S. L¨ock,P. Schlagheck, R. Ketzmerick, and M. K. Oberthaler,Phys. Rev. A , 011602 (2017).[22] R. Franzosi and V. Penna Phys. Rev. E , 046227(2003).[23] W.-C. Syu, D.-S. Lee, and C.-Y. Lin, Phys. Rev. D ,104011 (2019).[24] G. Thalhammer, G. Barontini, L. De Sarlo, J. Catani,F. Minardi, and M. Inguscio, Phys. Rev. Lett. ,210402 (2008).[25] S. Tojo, Y. Taguchi, Y. Masuyama, T. Hayashi, H. Saito,and T. Hirano, Phys. Rev. A , 033609 (2010).[26] J. Catani, L. De Sarlo, G. Barontini, F. Minardi, and M. Inguscio, Phys. Rev. A , 011603 (2008).[27] L. Wacker, N. B. Jørgensen, D. Birkmose, R. Horchani,W. Ertmer, C. Klempt, N. Winter, J. Sherson, and J. J.Arlt, Phys. Rev. A , 053602 (2015).[28] S. A. Moses, J. P. Covey, M. T. Miecnikowski, B. Yan,B. Gadway, J. Ye, and D. S. Jin, Science , 659 (2015).[29] S. B. Papp, J. M. Pino, and C. E. Wieman, Phys. Rev.Lett. , 040402 (2008).[30] F. Riboli and M. Modugno, Phys. Rev. A , 063614(2002).[31] K. Sasaki, N. Suzuki, D. Akamatsu, and H. Saito, Phys.Rev. A , 063611 (2009).[32] D. J. McCarron, H. W. Cho, D. L. Jenkin, M. P.K¨oppinger, and S. L. Cornish, Phys. Rev. A , 011603(2011).[33] Y. Eto, M. Kunimi, H. Tokita, H. Saito, and T. Hirano,Phys. Rev. A , 013611 (2015).[34] Y. Eto, M. Takahashi, K. Nabeta, R. Okada, M. Kunimi,H. Saito, and T. Hirano, Phys. Rev. A , 033615 (2016).[35] Alessia Burchianti, Chiara Fort, and Michele Modugno,Phys. Rev. A , 023627 (2017).[36] D. Ananikian and T. Bergeman, Phys. Rev. A , 013604(2006). [37] Becker, C., Stellmer, S., Soltan-Panahi, P. et al, NaturePhys , 496501 (2008)[38] T. Zibold, E. Nicklas, C. Gross, and M. K. Oberthaler,Phys. Rev. Lett. , 204101 (2010).[39] V. M. P´erez-Garc´ıa, H. Michinel, J. I. Cirac, M. Lewen-stein, and P. Zoller, Phys. Rev. Lett. , 5320 (1996).[40] L. Salasnich, A. Parola, and L. Reatto, Phys. Rev. A , 043614 (2002).[41] C. J. Pethick and H. Smith, Bose-Einstein Condensationin Dilute Gases , 2nd ed. (Cambridge University Press,Cambridge, 2008).[42] C.-Y. Lin, E. J. V. de Passos, and D.-S. Lee, Phys. Rev.A , 055603 (2000).[43] J. Grond, T. Betz, U. Hohenester, N. J. Mauser,J. Schmiedmayer, and T. Schumm, New Journal ofPhysics , 065026 (2011).[44] E. Boukobza, M. G. Moore, D. Cohen, and A. Vardi,Phys. Rev. Lett. , 240402 (2010).[45] G. P. Berman, A. Smerzi, and A. R. Bishop, Phys. Rev.Lett.88