Phase Separation and Dynamics of Two-component Bose-Einstein Condensates
aa r X i v : . [ c ond - m a t . qu a n t - g a s ] J un Phase Separation and Dynamics of Two-component Bose-Einstein Condensates
R. Navarro , R. Carretero-Gonz´alez , ∗ , P.G. Kevrekidis Nonlinear Dynamical Systems Group † , Computational Science Research Center,Department of Mathematics and Statistics, San Diego State University, San Diego CA, 92182-7720, USA Department of Mathematics and Statistics, University of Massachusetts, Amherst, MA 01003-4515 USA (Dated: Submitted to
Phys. Rev. A , 2009)We study the interactions between two atomic species in a binary Bose-Einstein condensate torevisit the conditions for miscibility, oscillatory dynamics between the species, steady state solutionsand their stability. By employing a variational approach for a quasi one-dimensional, two-atomicspecies, condensate we obtain equations of motion for the parameters of each species: amplitude,width, position and phase. A further simplification leads to a reduction of the dynamics into a simpleclassical Newtonian system where components oscillate in an effective potential with a frequencythat depends on the harmonic trap strength and the interspecies coupling parameter. We developexplicit conditions for miscibility that can be interpreted as a phase diagram that depends on theharmonic trap’s strength and the interspecies species coupling parameter. We numerically illustratethe bifurcation scenario whereby non-topological, phase-separated states of increasing complexityemerge out of a symmetric state, as the interspecies coupling is increased. The symmetry-breakingdynamical evolution of some of these states is numerically monitored and the associated asymmetricstates are also explored.
I. INTRODUCTION
Over the past 15 years, the study of Bose-Einstein con-densates (BECs) has gained a tremendous momentum,stemming from an intense and wide variety of theoreti-cal, as well as experimental studies that have now beensummarized in a number of books [1, 2]. One of the par-ticularly intriguing aspects of the system is its effectivenonlinearity stemming from a mean-field representationof the inter-atomic interactions. This, in turn, has ledto a wide range of developments in the area of nonlin-ear matter waves in BECs [3] and the drawing of naturalparallels between this field and that of nonlinear optics,where similar nonlinear Schr¨odinger (NLS) types of mod-els arise [4].One of the particularly interesting aspects of investi-gation of BECs within the realm of NLS (often referredto in the BEC context as Gross-Pitaevskii) equations isbased on the examination of multi-component systems.Starting from the early work on ground-state solutions[5, 6] and small-amplitude excitations [7] of the order pa-rameters, numerous investigations have been focused onthe study of two hyperfine states or two different atomicspecies that can be condensed and confined concurrently.More specifically, a few among the numerous topics in-vestigated involved the structure and dynamics of binaryBECs [8, 9, 10], the formation of domain walls betweenimmiscible species [8, 11], bound states of dark-dark [12],or dark-bright [13, 14], or coupled-vortex [15], or evenspatially periodic states [16]. The early experimental ef-forts produced such binary mixtures of different hyper-fine states of Rb [17] and of Na [19], but also of mixed ∗ URL: http://rohan.sdsu.edu/ ∼ rcarrete/ † URL: http://nlds.sdsu.edu condensates [18]. Efforts were later made to create two-component BECs with different atomic species, such as K– Rb [20] and Li–
Cs [21], among others. Re-cently the interest in multi-component BECs has beenrenewed by more detailed and more controlled experi-mental results illustrating the effects of phase separation[22, 23, 24], which have, in turn, motivated further the-oretical studies in the subject [25, 26]; see also [27] for arecent review.Although in the present work, we will focus on two-component condensates, it is relevant to note in passingthe increasingly growing interest in three-component, so-called, spinor condensates [28]. Among the numerousthemes of investigation within the latter context we men-tion spin domains [29], polarized states [30], spin textures[31], as well as multi-component (vectorial) solitons ofbright [32, 33, 34, 35], dark [36], gap [37], and bright-dark [38] types.Our aim in the present manuscript is to revisit thetheme of binary condensates in quasi-one-dimensionalBEC settings, in an attempt to offer additional both an-alytical and numerical insights on the phenomenology ofphase separation. Our manuscript is organized as fol-lows. In Sec. II, a Gaussian trial function is used ina variational approach to obtain six first-order ordinarydifferential equations (ODEs) for the time evolution ofthe parameters of the two-component ansatz: position,amplitude, width, phase, wave number and chirp. InSec. III A, the fixed points of the system of ODEs areobtained to yield the equilibrium position, amplitude,and width of the two species. Bifurcation diagrams ofthe ansatz’ parameters are produced as a function of theinterspecies coupling strength. In Sec. III B, phase dia-grams are produced and analytical conditions for mis-cibility are obtained relating the interspecies couplingwith the system’s chemical potential and parabolic trapstrength. The dynamics of this system is compared toresults obtained by numerical integration of the Gross-Pitaevskii equation in Sec. III C. In Sec. III D the sys-tem of ODEs is further reduced, upon suitable approx-imations, to a classical Newtonian system; the latter issimple to analyze and instructive with respect to the in-terpretation of the fundamental interactions driving thesystem’s dynamics. In Sec. IV, we numerically analyzein a systematic way the existence and stability of higherexcited phase-separated states as a function of the in-terspecies interaction. The dynamical instability evolu-tion of the latter class of states is examined numericallyin Sec. V. Motivated by the numerical experiments ofSec. V, in Sec. VI we study the existence of asymmetricstates when the chemical potentials of the two compo-nents differ. Finally, in Sec. VII we summarize our re-sults and present some interesting directions for futureinvestigation.
II. VARIATIONAL MODELA. Coupled equations
The Gross-Pitaevskii (GP) equation, which is a variantof the NLS equation accounting for the potential confin-ing the atomic species, governs the dynamics of bosonicparticles near absolute zero temperatures [1, 2]. In thecontext considered herein (related to the case of Rbwhich is common in relevant experiments [17, 23]), thetwo hyperfine states of the same atom are described bya set of coupled GP equations [9] i ~ ∂ψ ∂t = (cid:18) − ~ m ∇ + V + g | ψ | + g | ψ | (cid:19) ψ ,i ~ ∂ψ ∂t = (cid:18) − ~ m ∇ + V + g | ψ | + g | ψ | (cid:19) ψ ,V j = 12 (cid:0) ω jx x + ω jy y + ω jz z (cid:1) , j = 1 , . (1)Here, g jk = 4 π ~ a jk /m are the self coupling interactionparameters of the first species for j = 1, k = 1, for thesecond species j = 2, k = 2 and for j = 1, k = 2 or j = 2, k = 1 are the cross species coupling parameters. Theseinteraction strengths depend on the scattering lengthsbetween same ( a and a ) and different species ( a = a ). The external harmonic trapping potential for eachatomic species, V j = V j ( r ) depends on the radial distance r from the center of the trap, while the atomic densityis given by the square modulus of ψ j = ψ j ( r , t ). Theatomic mass is denoted by m .We will make the customary assumption that the ex-ternal trap’s effect on each species is the same, making V ( r ) = V ( r ). More importantly, in the interest of an-alytical tractability of our results, we will also assumethat the self interactions for each species is the same, g = g . Both assumptions are very good approxima-tions of the physical reality, although they are not exact;see e.g., the relevant discussion of [23]. In a highly anisotropic trap, where the frequency ofthe longitudinal component of the trap is much smallerthan the transverse components ω x ≪ ω y = ω z , an effec-tive one-dimensional (1D) system of partial differentialequations (PDEs) can be obtained [3] i ∂u ∂t = (cid:18) − ∂ ∂x + Ω x + | u | + g | u | (cid:19) u , (2) i ∂u ∂t = (cid:18) − ∂ ∂x + Ω x + | u | + g | u | (cid:19) u , (3)where time, space, and wave function have been rescaledto reduce the system’s parameters to just two (Ω and g ).In this 1D reduction, the chemical potential correspondsto µ D , ψ j → u j p g jj /µ D , x → x √ mµ D / ~ is the longi-tudinal distance from the center of the trap, t → tµ D / ~ is the rescaled time, g = g /g = g /g is the rescaledspecies interaction term. In the literature, the conditionof miscibility ∆ = ( g g − g g ) /g = g − ~ ω x /µ D . B. Ansatz and Euler-Lagrange Equations
To develop a variational model, we substitute in therescaled two-component Lagrangian L = Z ∞−∞ ( L + L + L + L ) dx, where L j = E j + i (cid:18) u j ∂ u ∗ j ∂t − u ∗ j ∂u j ∂t (cid:19) ,E j = 12 (cid:12)(cid:12)(cid:12)(cid:12) ∂u j ∂x (cid:12)(cid:12)(cid:12)(cid:12) + V ( x ) | u j | + 12 | u j | ,L = L = 12 g | u | | u | , the Gaussian ansatz of the form u ( x, t ) = Ae − ( x − B )22 W e i ( C + Dx + Ex ) , (4) u ( x, t ) = Ae − ( x + B )22 W e i ( C − Dx + Ex ) . (5)Starred variables indicate complex conjugation. The pa-rameters A , B , C , D , E , and W are assumed to be timedependent and they represent the amplitude, position,phase, wave number, chirp, and width of the Gaussianansatz, respectively. For large Ω, the steady state solu-tion very closely resembles two Gaussian functions sep-arated by a distance 2 B , with constant rotation of thephase C = µt , where µ is the chemical potential (with-out loss of generality we take µ = 1). Upon interactionbetween the species and with the harmonic trap, accel-eration induces an inhomogeneity in the carrier wave,known as chirp, that accordingly affects the phase of eachspecies. It is important to note here that the ansatz isinvariant upon transposition of space and atomic speciescomponent x → − x and u → u . This allows us toreduce the number of parameters in the system to sixinstead of having twelve, which would take into accountindependent variation of the parameters in u and u .However, this simplification comes at a certain cost as,in particular, it is not possible to monitor asymmetric(between the two components) states within this ansatz;the latter type of states will be partially explored withinthe dynamics of the species in Sec. V.When the Lagrangian is evaluated for the proposedansatz, a spatially-averaged effective Lagrangian is ob-tained L = −√ πA W (cid:20) Ω B + 2 E W + ( D + 2 BE ) + A √ ge − B W ) + 12 W + Ω W dCdt + 2 B dDdt + dEdt (2 B + W ) (cid:21) , (6)and the equations of motion for the parameters are ob-tained through the corresponding Euler-Lagrange equa-tions ∂L∂p j − ddt (cid:18) ∂L∂ ˙ p j (cid:19) = 0 , (7)where the parameter p j , j = 1 , , . . . , A , B , C , D , E , W and ˙ p j = dp j /dt . The Euler-Lagrange equations for A , B , and W evaluate to ∂L/∂A = 0, ∂L/∂B = 0 and ∂L/∂W = 0since the second term in the Euler-Lagrange equation(7) is zero for these. This results in equations that in-volve linear combinations of dC/dt , dD/dt and dE/dt that when solved, give Eqs. (10), (11) and (12). Therest of the Euler-Lagrange equations give equations thatinvolve linear combinations of dA/dt , dB/dt and dW/dt that can be solved to give Eqs. (8), (9) and (13). Thefollowing equations are the result of solving the Euler-Lagrange equations for the time derivatives of all the pa- rameters of our ansatz: dAdt = − AE, (8) dBdt = D + 2 BE, (9) dCdt = B W − D − W + √ A W (2 B − W ) + √ A g W e − B W (8 B + 2 B W + 5 W ) , (10) dDdt = √ A Bg W e − B W (4 B + W ) −√ A b W − BW − DE, (11) dEdt = √ A g W e − B W ( − B + W ) + √ A W + 12 W − E − Ω , (12) dWdt = 2 EW. (13)As described below, these ODEs reflect fairly accu-rately the principal dynamical features of the system.In particular, they capture the oscillations of the twospecies past each other when the equilibrium separation B is zero and the oscillations about their correspondingphase-separated equilibrium position when the two com-ponents are phase-separated. III. STEADY STATE SOLUTIONSA. Phase Bifurcations
The steady state of Eqs. (8)–(13) is obtained by setting dA/dt = dB/dt = dC/dt = dE/dt = dD/dt = dE/dt = dW/dt = 0 witch leads immediately to the steady statesolution E ∗ = D ∗ = 0 and C ∗ = µ . When the equilibriumseparation between the species is zero, the equilibriumamplitude and width reduce to B ∗ = 0 , (14) A ∗ = 2 √ (cid:16) µ − p + 4 µ (cid:17) g ) , (15) W ∗ = (cid:16) µ + p + 4 µ (cid:17) . (16)When the equilibrium separation is nonzero, the result-ing equilibrium amplitude and width are given by the −4 −2 0 2 40.10.20.30.40.50.6 x u and u (a) PDE ODE−4 −2 0 2 40.10.20.30.40.5 x u and u (b) PDE ODE FIG. 1: (Color online) (a) Steady state solution, u and u , forthe mixed state when B = 0. (b) Steady state solution for theseparated state when B = 0. The solid (blue) line representsthe steady state of the GP equations and the dashed (red)line is the steady state solution of the system of ODEs. Here,Ω = 0 . µ = 1 g = 1 for (a) and g = 20 for (b). transcendental relationsΩ − √ A ∗ gW ∗ e − B ∗ W ∗ = 0 , (17) µ − W ∗ − W ∗ (cid:16) √ A ∗ g + Ω W ∗ (cid:17) = 0 , (18) µ + 34 W ∗ −
54 Ω (cid:0) W ∗ + 2 B ∗ (cid:1) = 0 . (19)Figure 1 shows that the steady state solution of thefull GP model closely matches the steady state from theODEs, indicating that the ansatz successfully capturesthe relevant PDE behavior. For large values of Ω, thesteady state solution of the GP deviates from the Gaus-sian shape, resembling an inverted parabola, often re-ferred to as the Thomas-Fermi approximation [1, 2].Detailed bifurcation diagrams for the amplitude, po-sition, and width of each species can be obtained by −1−0.8−0.6−0.4−0.200.20.40.60.81 B ACCBBD(a)
Unstable (ODE)Stable (ODE)Unstable (PDE)Stable (PDE) A ACBD(b)
Unstable (ODE)Stable (ODE)Unstable (PDE)Stable (PDE) W ACBD(c)
Unstable (ODE)Stable (ODE)Unstable (PDE)Stable (PDE)
FIG. 2: (Color online) Equilibrium (a) separation, (b) am-plitude, and (c) width for the two condensed species, as afunction of the interspecies coupling strength g for Ω = 0 . µ = 1. The (red) thin solid (stable) and dashed (un-stable) lines correspond to the steady states for the reducedODE model [Eqs. (8)-(13)], while the (blue) thick solid (sta-ble) and dots (unstable) correspond to the steady state fromthe full GP model [Eqs. (2)-(3)]. solving Eqs. (17)–(19) for A , B , and W for the phaseseparated state and Eqs. (14)–(16) for the mixed state.The steady state solution reveals a pitchfork bifurcationfor the position of each condensate as the interspeciescoupling strength is increased as can be seen in Fig. 2a.Interestingly, the steady state of the GP equations pro-duces a supercritical pitchfork bifurcation at point D ofFig. 2. More specifically, for small values of g , a sta-ble mixed phase can be identified; as g is increased, themixed phase becomes unstable and the phase separatedstate becomes stable. On the other hand, the system ofODEs also predicts a pitchfork bifurcation, however theapproximate nature of the ansatz results in the identifi-cation of the bifurcation as a subcritical one (point A) oc-curring in the vicinity of a symmetric pair of saddle nodebifurcations (points B). Both bifurcation diagrams agreewith each other away from the transition between phases.In the vicinity of the transition point, clearly, the natureof the ansatz is insufficient to capture the fine details ofthe bifurcation structure (thus inaccurately suggesting aphenomenology involving bistability, and hysteretic be-havior of the system). It should be noted that the phaseseparation from the GP model (point D) lies near thesaddle node bifurcation point (point B) and the subcriti-cal bifurcation point (point A) from the system of ODEs.At small values of the harmonic trap strength, the truepoint of phase transition is closer to the subcritical bi-furcation point and at larger values of the harmonic trapstrength, the true point of phase transition lies closer tothe saddle node bifurcation point. B. Phase Separation
Using the zero separation amplitude and width, pointA in Fig. 2, an expression for the onset of phase separa-tion can be approximated in terms of the critical inter-species interaction g cr = 6 µ + 3 p + 4 µ µ − p + 4 µ . (20)Despite the deviation of point A in Fig. 2, from the rel-evant point D of the corresponding PDE, the analyticalexpression offers valuable insight on the dependence ofthe critical interspecies interaction for phase separationon parameters of the trap (in particular, its frequency)and those of the condensate (in particular, its chemicalpotential); see also Fig. 3 for a detailed comparison ofthe ODE and PDE bifurcation points. More specifically,the equation predicts that when the harmonic trap’s fre-quency approaches zero, Ω →
0, phase separation occurswhen g cr → g − → Ω cr = 4 √ µ/ ≈ . µ phase separation occurs at g cr → ∞ . This behavior canbe qualitatively understood since tighter (larger Ω) trapstend to “squeeze” both components together, thus frus-trating the system’s tendency towards phase separation.The prediction of the system of ODEs for the loca-tion of the bifurcation point agrees well with the resultsfrom the GP model for small values of the harmonic trapstrength. Recall that the phase separation for the ODEmodel is located at the subcritical pitchfork bifurcationpoint A in Fig. 2. However, as can be observed from Ω g Stable separated state Stable mixed state
Supercritical Pt. D (PDE) Subcritical Pt. A (ODE) Saddle−Node Pt. B (ODE)
FIG. 3: (Color online) The (blue) solid line represents theboundary of zero species separation from the GP equations(point D in Fig. 2), where states that lie to the right of thisline are mixed and values to the left are phase-separated. The(red) dashed line represents the boundary of zero species sep-aration for the system of ODEs given by Eq. (20) (point A inFig. 2). The (black) dotted line shows the saddle-node point(point B in Fig. 2) obtained from the system of ODEs.FIG. 4: Oscillations for a mixed state from direct numericalintegration of the GP model (2)-(3) for g = 1 ( g < g cr ), µ = 1and Ω = 0 .
6. Lighter gray corresponds to one the density ofone component and darker gray to the other component.
Fig. 3, a better approximation for the full system’s phaseseparation (point D), in the case of large trapping fre-quencies, can be given by the saddle-node bifurcationpoint B for the ODE reduced system for Ω values closeto Ω cr . A (a)−101 B (b)00.51 C + µ t (c)−0.500.5 D (d)−0.200.2 E (e) W t(f) FIG. 5: (Color online) (a) Amplitude, (b) position, (c) phase,(d) velocity, (e) chirp, and (f) width corresponding to the mixed oscillating state show in Fig. 4. Solid lines representresults from direct numerical simulations of the GP equationwhile dashed lines depict the results for the ODE reduction(9)–(13).FIG. 6: Same as in Fig. 4 for the time evolution of oscillatingspecies about their phase separated configuration when g =20 ( g > g cr ). A (a)−101 B (b)−2.5−2−1.5−1−0.5 C + µ t (c)−0.500.5 D (d)−0.200.2 E (e) W t(f) FIG. 7: (Color online) Same as in Fig. 5 for the time evolutionof oscillating species about their phase separated configura-tion as shown in Fig. 6.
C. Dynamics of the Reduced System
For relatively small values of the interspecies coupling g (i.e., to the left of point D in Fig. 2) the two species donot separate and thus, when given opposite direction ve-locities from the mixed state, they will oscillate througheach other as depicted in Fig. 4. This case is analyzed inFig. 5 where the two components oscillate through eachother about their common equilibrium separation of zeroand the prediction from the system of ODEs (red dashedlines) agrees very well with the direct integration of theGP equations (blue solid lines). Because of conservationof mass, the amplitude and width oscillate out of phase:the amplitude is maximized and width is minimized whenthe acceleration of the species is maximized; on the otherhand, the width is maximized and the amplitude is min-imized when the velocity is maximized. The velocity ofeach species [see Eq. (9)] has two components, one thatdepends on the wave number D and another that dependson the product of chirp and position 2 EB . If there is nochirp, the wave number is the velocity of the condensatespecies. In Fig. 5c and 7c, a factor of µt has been addedto show the deviation of the phase from the steady statevalue.On the other hand, if we assume a well-separated stateas the PDEs initial condition, (right of point D in Fig. 2),it is possible that the two components entertain oscil- −5 0 500.511.522.533.544.5 B U e ff g=0g=2.6g=6g=20 FIG. 8: (Color online) The effective potential, U eff , for severalvalues of the species coupling parameter: g = 0, g = 2 . g = 6, and g = 20. Here, Ω = 0 . µ = 1. lations about their phase-separated steady state as de-picted in Fig. 6. As illustrated in Fig. 7, the two phaseseparated components collide against each other as theyoscillate (but do not go through each other) about anonzero position. The prediction from the system of or-dinary differential equations once again agrees very wellwith the numerical integration of the GP equation forthe position of the two species, but differs somewhat forthe other parameters of motion. This occurs because thetime evolution of the solution of Eqs. (2) and (3) in thiscase, due to the oscillation and interaction, deviates froma Gaussian waveform and the corresponding variationalprediction begins to lose accuracy. D. Newtonian Reduction
To develop a more tractable model for the dynamics, aclassical Newtonian system for the motion of the centerof each species is desirable. Taking a time derivative ofEq. (9) and substitution of Eqs. (9), (11), and (12) yields d Bdt = − Ω − √ A ∗ gW ∗ e − B W ∗ ! B, (21)where we simplified the dynamics by replacing the timedependent variations of A and W by their respective equi-librium values A ∗ and W ∗ . This simplification is justifiedby the fact that the oscillations in A and W are relativelyweak as it can be observed from Figs. 5a, 5f, 7a, and7f. These phase separated oscillations contain two funda-mental physical features: the external trapping potentialand an exponential repulsive interaction that dependson the cross species coupling parameter g . IntegratingEq. (21) with respect to B yields a Newtonian equation A m p li t ude v (a)System of ODEsNewtonian ODEPDE T v (b)System of ODEsNewtonian ODEPDE FIG. 9: (Color online) (a) Amplitude of oscillation and (b)period of oscillation, T , as a function of increasing initial ve-locity, v . The (blue) solid line represents results from thecoupled GP equations, the (red) dashed line represents re-sults from the system of ODEs, and the black dotted linerepresents the results from the Newtonian reduction. Here, g = 1, µ = 1 and Ω = 0 . of motion under the effective potential U eff = Ω B + √ A ∗ g e − B W ∗ . (22)This reduced dynamics gives an effective double well po-tential for a fully phase separated state ( g > g cr ) anda nonlinear single well potential for the mixed states( g < g cr ). It is important to note that in the expres-sion of Eq. (22), A ∗ and B ∗ vary as a function of g . Thisdependence has been incorporated in Fig. 8 which showsthe effective potential for increasing values of g , where ittransitions from a single well potential to a double wellpotential. For Ω = 0 . µ = 1, the reduced Newtonianmodel predicts that at g = 2 . T (a) ODEPDEEigenvalues T (b) ODEPDEEigenvalues FIG. 10: (Color online) (a) The period of oscillation, T , isshown as a function of g for an initial velocity of v = 0 . g for an initialvelocity of v = 0 .
1. The (blue) solid line represents resultsfrom the coupled GP equations, the (red) dashed line rep-resents results from system of ODEs, while the black dotsrepresent the results from the evaluation of the eigenvaluesof the Jacobian of the system of ODEs at equilibrium. Here,Ω = 0 . µ = 1. of two particle-like excited states (i.e., dark solitary mat-ter waves) within the same species and has been found tobe extremely successful in comparisons with experimen-tal results [39]. In that case, as well, the fundamentalcharacteristics were the parabolic trapping and the ex-ponential tail-tail interaction between the waves.Figure 9 depicts both the (a) amplitude and (b) periodof oscillations predicted by the Newtonian reduction incomparison with the corresponding PDE findings. It isseen that the GP equation and the ODE system agreevery well for a wide range of values of initial veloci-ties. This figure also shows the nonlinearity of oscillationswhere small amplitude oscillations have a period of 18.3and large oscillations yield a period of T ≈ .
5, whichcorresponds to the harmonic trap’s period T → π/ Ω. Figure 10a shows that the period of oscillation from thesystem of ODEs matches that of the GP model for smallvalues of g . As g → . g → .
1, the sys-tem of ODEs begins to phase-separate. Since these aresmall oscillations, the eigenvalues from the Jacobian ofthe system of ODEs at equilibrium matches very wellthe oscillations of the system of ODEs. Figure 10b showssimilar results to Fig. 10a but for larger oscillations. Wecan see that for larger oscillations the ODEs’ period moreclosely matches that of the GP for smaller g , while forlarger values of the interspecies strength, the deviationbecomes more significant. Furthermore, for larger val-ues of g , the period obtained from the eigenvalues alsodeviates from the period from the ODEs. For large oscil-lations, the two species do not effectively interact (sincethey go through each other too rapidly to “feel” eachother) and the period is roughly independent of g as pre-dicted by Eq. (21) when | B | ≪
1, yielding simple har-monic oscillations. In this case, the GP’s period is closeto that predicted by the ODEs.
IV. EXCITED STATES
For sufficiently weak traps Ω < .
5, as the interspeciescoupling strength is increased for g >
1, new, high ordermixed states emerge. These higher excited states cor-respond to alternating bands dominated successively byeach of two species. As Ω is decreased or alternatively g is increased, the number of alternating bands increaseswithin the solution profile and the population imbalancewithin each band less pronounced. Each solution branchis found by using parameter continuation on the param-eter g using a Newton fixed point iteration to find thestationary state. The stability for each computed pro-file was obtained by computing the eigenvalues of thelinearized dynamics at the fixed points using standardtechniques [3]. We now describe the series of bifurca-tions that occur as the interspecies coupling g increasesas depicted in Fig. 11: • For g < .
08, the only state that exists is the mixedstate and it is stable. This threshold is close to thetraditional miscibility condition ∆ = g − • At g = 1 .
08 (point A), the mixed states becomesunstable, through the previously described super-critical pitchfork bifurcation, leading to the emer-gence of a stable phase-separated state where onecomponent moves to the left and the other onemoves to the right. • At g = 1 .
24 (point B), a second supercritical pitch-fork occurs, rendering the symmetric (across com-ponents) solution more unstable and giving rise toa state where one species has a single hump and theother one has a double hump. We call this state a1-2 hump configuration. ∆ h gA B C D E FIG. 11: (Color online) Degree of phase separation ∆ h = u ( x c ) − u ( x c ), where x c is the location of the closest densitymaximum (irrespective of component) to the trap center, as a function of the interspecies coupling parameter for variousexcited states for g ∈ [0 . , .
5] and Ω = 0 .
2. Stable (unstable) solutions branches are depicted with black solid (dashed) lines.Typical solutions for each branch are depicted in the surrounding insets. λ r A B C D E
FIG. 12: (Color online) Real part of largest eigenvalue of themixed state as a function of the interspecies coupling param-eter for the various excited states depicted in Fig. 11 with g ∈ [0 . , . . • At g = 1 .
63 (point C), the third supercritical pitch-fork bifurcation of the series arises, leading thistime to a situation where one species moves to theleft and the other one moves to the right, both forming a double hump (a 2-2 hump configuration). • At g = 2 .
22 (point D), a double hump with a triplehump state arises (a 2-3 hump configuration). • At g = 3 . g is increased, higher “modulationwavenumbers” become unstable, leading to the “quan-tized” (associated with the quantization of wavenum-bers in the effectively finite box) cascade of supercriticalpitchfork bifurcations and associated further destabiliza-tions of the symmetric state. The degree of phase sep-aration in Fig. 11 between these bands is computed as∆ h = u ( x c ) − u ( x c ), where x c is the location of theclosest density maximum, irrespective of component, tothe trap center.The emergence of high order states can be inferredby observing the real part of the eigenvalues as the in-terspecies coupling is varied. The eigenvalues collide0with the imaginary axis as new states emerge. Eigen-value analysis shows that all excited states are unstablewith the exception of the single-hump phase-separatedstate (which results from the first pitchfork bifurcation,namely the 1-1 hump state). These excited states are,however, less unstable than the mixed state from whichthey arise. The latter, as shown in Fig. 12, becomes pro-gressively more unstable, as expected, as further multi-hump branches arise. Each of the bifurcation points, forwhich the same designation as in Fig. 11 is used, corre-sponds to a further pair of real eigenvalues appearing forthe mixed branch. g λ r BCD E 1−2 hump2−2 hump2−3 hump3−3 hump
FIG. 13: (Color online) Real part of the largest eigenvaluesof the principal phase-separated states as a function of theinterspecies coupling parameter, g , for Ω = 0 . As indicated by Fig. 11, the 1-1 hump phase separatedstate is stable even for large values of g . In fact, forΩ ≪ g ≫
1, the two components repel each otherso strongly that the center of magnetic trap becomesa domain wall, where each species abruptly transitionsfrom near zero atomic density to maximum atomic den-sity. Apparently, this 1-1 hump phase-separated stateis the only stable state of the system after the mixedstate loses its stability past the bifurcation point A (seeFig. 11). Nonetheless, as depicted in Fig. 13, the addi-tional multi-humped excited states are successively cre-ated as g is increased (bifurcation points B, C, D, andE in Fig. 11) and are progressively less unstable as g increases. Therefore, for sufficiently large g , the instabil-ity of higher excited multi-humped states might be weakenough for these states to be observable within experi-mentally accessible times. It is interesting to note thatsince these excited multi-humped states emanate fromthe (already unstable) mixed state, they feature a rela-tively strong instability close to their bifurcation point. (b) −5 0 500.20.4 x | u , ( x ) | t= 0t=40 (c) −5 0 500.20.4 x | u , ( x ) | t= 0t=145 (e) −5 0 500.20.4 x | u , ( x ) | t=0t=10000 FIG. 14: (Color online) Dynamics of the mixed state for g = 5and Ω = 0 .
2. (a) Top and bottom subpanels correspond tothe evolution of the densities for the first and second compo-nents, respectively, after applying an initial spatially randomperturbation of size 1 × − . Panels (b) and (c) depict snap-shots of the initial density and at the times indicated in panel(a) by the white vertical dashed lines. (d) Long term dy-namics showing that the mixed state eventually approaches aseparated state. (e) Snapshots of the densities at t = 0 and t = 10 ,
000 (see white vertical dashed line in panel (d)).
V. DYNAMICS OF UNSTABLE STATES
In this section we present the dynamics of the differ-ent unstable states that were identified above. We startby analyzing the destabilization of the mixed state fora value of g to the right of the bifurcation point A (seeFig. 11). In fact, for all the dynamical destabilizationresults presented in this section we chose a value of g = 5and Ω = 0 . (b) −5 0 500.20.40.60.8 x | u , ( x ) | t= 0t=758 (c)0 1000 2000 3000 4000 500010 t || | u i ( t ) | − | u i ( ) | || FIG. 15: (Color online) Dynamics of the unstable 1-2 humpstate for g = 5 and Ω = 0 .
2. (a) Same as in Fig. 14.a. (b)Snapshots of the densities at t = 0 and t = 758 (see whitevertical dashed line in panel (a)). (c) Growth of the norm ofthe perturbation vs. time (in semi-log plot). Circles (crosses)correspond to the perturbation of the first (second) compo-nent and the red solid line is the instability growth obtainedfrom the eigenvalue computation (max( ℜ ( λ ))) = 0 . E (see Fig. 11). Therefore, in this regime we have co-existence of several unstable multi-hump solutions andthe stable phase-separated 1-1 hump state. Other pa-rameter regions (not shown here) gave qualitatively thesame results.In Fig. 14 we show the destabilization of the mixedstate. As it can be observed from the figure, the mixedstate suffers an initial modulational instability that be-comes apparent for t >
35 that seeds a highly perturbed2-3 hump solution [see panel (b)]. Since this new solu-tion is also unstable for the chosen parameter values, itis rapidly converted into a relatively long lived 1-2 humpstate [see panel (c)]. Nonetheless, as it is clear from thelong term dynamics presented in panel (d), the 1-2 humpstate, being unstable, eventually “decays” to the sepa-rated (1-1 hump) state. We use here the term “decaying”in quotes since our system is conservative (no dissipation)and thus there is no real decay.In Fig. 15 we show the destabilization of the unstable1-2 hump state. As it is obvious from the figure, al-though the instability eigenvalue is sizeable [ λ r = 0 . (b) −5 0 500.20.40.60.8 x | u , ( x ) | t= 0t=4500 FIG. 16: (Color online) Dynamics of the unstable 2-2 humpstate for g = 5 and Ω = 0 . (b) −5 0 500.20.4 x | u , ( x ) | t= 0t=149 (d) −5 0 500.20.4 x | u , ( x ) | t=0t=50000 FIG. 17: (Color online) Dynamics of the unstable 2-3 humpstate for g = 5 and Ω = 0 . t > , (b) −5 0 500.20.4 x | u , ( x ) | t= 0t=74 (d) −5 0 500.20.4 x | u , ( x ) | t=0t=15000 FIG. 18: (Color online) Dynamics of the unstable 3-3 humpstate for g = 5 and Ω = 0 . t > , for a short period of time (becoming slightly asymmetric)and then comes back close to the 1-2 hump (symmetric)steady state configuration. This indicates that, for thisparameter combination, the 1-2 hump state is a saddlefixed point. Thus, the orbit remains close to the steady1-2 hump state and it is eventually “kicked out” along theunstable manifold. Then, it performs a quick excursionand “comes back” through the stable manifold.In Fig. 16 we depict the destabilization dynamics of the2-2 hump state. As it is evident from the figure, whencompared to Fig. 15, the destabilization happens earlier( t ≈ t ≈ t ≈
90 (earlier than its 1-2 and 2-2 hump counterparts given its larger instability growthrate [see Fig.13]). The dynamics goes through a transient1-2 state and eventually settles into a highly perturbedseparated state that is preserved thereafter as was thecase also for the mixed state dynamics (see Fig. 14) andthe 2-2 state (see Fig. 16) presented above.Finally, in Fig. 18 we present the dynamical destabi-lization for the 3-3 hump state. Again, this state desta-bilizes even faster ( t = 40) than the previous waveformsbecause of its higher instability eigenvalue (see Fig. 13).In this case, the 3-3 hump state first destabilizes into ahighly perturbed 2-2 state. Since the 2-2 hump state isunstable, the ensuing dynamics results into a separated1-2 hump state after t >
100 which, in turn, eventuallysettles to a highly perturbed separated state. (a)0 25 50 75−0.500.5 t−50000 R e ( u , ) (b)0 25 50 75−0.500.5 t−100000 R e ( u , ) FIG. 19: (Color online) Dynamics of the real part of the wavefunctions for the 2-3 hump state of Fig. 17 for g = 5 andΩ = 0 .
2. Panels (a) and (b) correspond, respectively, to timescentered about t = 50 ,
000 and t = 100 , VI. EXISTENCE AND STABILITY OFASYMMETRIC STATES
It is interesting to observe in the above computationsthat the transient 1-2 hump state emanating from thedestabilization of the mixed state (see Figs. 14) or fromhigher order states (see Figs. 17 and 18) has a relativelylong life span. However, being this state also unstable,it eventually tends to a perturbed separated (i.e., 1-1hump) state. The process of converting a 1-2 state intothe separated state involves the effective shift of massfrom one of the two humps (in the component with twohumps) to the other one until all the mass is “swallowedup” by one hump resulting in a 1-1 hump. For exam-ple, as it can be observed in panel (d) of Fig. 17, at t = 50 ,
000 the right hump of the first component has3 µ , (a) 0.6 0.8 1 µ , (b) 0 25,000 50,000 75,000 100,000 0.6 0.8 1 t µ , (c) FIG. 20: (Color online) Local chemical potential at the loca-tion of maximum density for each component. Blue solid linecorresponds to the first component while the green dashed linecorresponds to the second component. The different panelscorrespond to the evolution of the local chemical potentials forthe steady states: (a) mixed state (cf. Fig. 14), (b) 2-3 humpstate (cf. Fig. 17), and (c) 3-3 hump states (cf. Fig. 18). more mass than its left hump. Since the chemical po-tential (i.e., rotation frequency of the wave function inthe complex plane) is closely related to the mass of thecondensate, it is possible to follow the local change inmass by following the local chemical potential. In Fig. 19we depict the oscillations of the real part of the wavefunctions at the location of the maximal density (seesquares and circles in the figure). We then fit a sinu-soidal curve (see solid lines in the figure) through thedata to obtain the oscillation frequency. As it is obviousfrom the figure, although both components started withthe 2-3 hump steady state with the same chemical po-tential µ = µ = µ = 1, the two components oscillate atdifferent rates. In fact, at t = 50 ,
000 the local chemicalpotentials for each species is, respectively, µ = 0 .
776 and µ = 0 . t = 100 ,
000 we have µ = 0 .
688 and µ = 0 . µ = µ = 1but rapidly drop when the initial instability of the steadystate develops for the three cases. Then, the chemical po-tentials slowly drift until they acquire values very close toeach other at around the time when the dynamics settlesto a perturbed phase separated state.The fact that during the dynamical evolution of theunstable higher order states the local chemical poten-tials vary naturally prompts the question of existence and stability of steady states with different chemical poten-tials between the species. Up to this point all the steadystates we analyzed supposed that both components hadthe same chemical potential. Namely, the evolution foreach components is u ( x, t ) = w ( x ) e − iµ t u ( x, t ) = w ( x ) e − iµ t , where the steady state profiles are w , ( x ) and µ = µ =1. We now relax this symmetric assumption and look forsteady states with different chemical potentials betweenthe components ( µ = µ ). A full bifurcation diagram asa function of µ , for fixed µ = 1, is depicted in Fig. 21.In this diagram we use in the vertical axis the differencein the variance of the steady state distribution betweenthe two components. The different bifurcation branchescorrespond to the particular asymmetric states depictedin the surrounding insets. This diagram provides an ac-count of how the different solutions bifurcate from eachother. It is very interesting to follow the bifurcation pathof all the states labeled from A to T during which all thehigher order states are browsed continuously. Due tosymmetry, a similar scenario is present when one keeps µ = 1 fixed and varies µ (results not shown here).The stability in the bifurcation diagram depicted inFig. 21 is denoted, as before, with a solid line for stablestates and a dashed line for unstable states. As it can benoticed from the figures, as it was the case for symmetricstates, all the asymmetric states are unstable except themixed state. This is due to the fact that the interspeciescoupling for these diagrams was chosen as g = 5, i.e.,high enough so that the mixed state in unstable. In fact,in Fig. 22 we show the largest real part of the eigenvaluesfor the different asymmetric states depicted in the bifur-cation diagram in Fig. 21. As it was observed for thesymmetric case, the higher the order of the asymmetricstate the more unstable it becomes. It is crucial to notethat the 1-2 hump state seems to get stabilized for µ = 1and µ > .
12 (see top panel of Fig. 22). However, af-ter close inspection, see the bottom panel in Fig. 22, thereal part of the eigenvalue for the 1-2 hump state never vanishes but instead becomes extremely small (between10 − and 10 − ) until the branch disappears (at µ ≈ . µ = 1). This very weak instability explains the factwhy the transient 1-2 hump state appears to be sustainedfor very long times in the evolutions depicted in Figs. 14,17, and 18. VII. CONCLUSIONS
In the present work, we have analyzed the emergenceof non-topological, phase-separated states in the immis-cible regime out of mixed ones in the miscible regimes,as a natural parameter of the system (namely the inter-species interaction strength) was varied. Our analysis4
A B C D E FGHIJKLMNOPQRST 0.6 1 1.4 1.8 −0.1 0 µ V a r( u ) − V a r( u ) ABC,K,RDE HGF J IL N 0M PQST
FIG. 21: (Color online) Bifurcation diagram of asymmetric states for constant µ = 1 as a function of µ . The vertical axiscorresponds to the difference of the variance of the steady state configurations for both components for g = 5 and Ω = 0 . µ λ r u u u u u −8 −6 −4 −2 µ λ r u FIG. 22: (Color online) Real part of largest eigenvalue forthe asymmetric mixed states as a function of µ with µ = 1, g = 5 and Ω = 0 .
2. The bottom panel corresponds to theeigenvalue for the 1-2 hump state plotted in logarithmic scale. was presented for the case of magnetically trapped two-component Bose-Einstein condensates i.e., two incoher-ently coupled NLS equations with a parabolic potential.Using a variational approach, we obtain a reduced modeldescribing the statics, stability and dynamics of each con-densate cloud. We are able to elucidate the miscibilityboundary for the interspecies coupling parameter as thestrength of the magnetic trap (and/or the chemical po-tential of the system) is varied. The approach is alsocapable of accurately capturing the spatial oscillationsof the clouds about the stable stationary states (bothfor mixed and for phase-separated states). In particular,for relatively small interspecies coupling, the two con-densed clouds do not phase-separate (mixed state) giv-ing rise, if the BEC clouds are initially displaced withrespect to each other, to oscillations through each other.On the other hand, for relatively large interspecies cou-pling, the two clouds form a stable phase-separated statewhich can entertain oscillations about the equilibriumseparation between the components. A further dynam-ical reduction allows to understand this behavior moreintuitively based on an effective potential that undergoesa bifurcation from a single well (mixed state) to a dou-ble well (phase-separated state) form as the interspeciescoupling parameter is increased. We also describe the bi-furcation scenario of higher order phase-separated states,5as the interspecies coupling parameter is increased. Weobserve that several (interwoven between the two compo-nents) bands of density modulations progressively ariseout of the mixed states giving rise to higher excitedstates. Among all the phase-separated states, only thefirst excited one with one hump in each component isfound to be dynamically stable for all values of the in-terspecies interaction strength past its bifurcation point.We furthered our analysis by studying the existence andstability of asymmetric states for which the chemical po-tentials for each species is different. We found, similarto the symmetric case of equal chemical potentials, thatthe only stable steady state for high enough interspeciescoupling is the separated state. Nonetheless, we foundthat in some regimes, the (asymmetric variant of the)state with one band in one component surrounded by twobands of the other component (the 1-2 hump state) canhave an extremely weak instability and thus facilitatingits potential observability in numerical experiments.There are numerous directions along which the presentwork can be naturally extended. For example, within theone-dimensional context, it is natural to seek to relaxthe simplifying assumption of the intraspecies scatteringlengths (and by extension the self-interaction coefficientsof the two components, g and g ) being equal. How-ever, this extension will unfortunately have to come at the expense of an ansatz with different amplitude, width,etc. parameters between the components, which will ren-der the intuitive and explicit analytical results obtainedherein much more tedious.Another natural extension is to try to generalize theideas presented herein into higher-dimensional or largernumber of component (i.e., spinor) settings. Especiallyin the former case, more complicated waveforms suchas crosses and propellers have been predicted in two-component condensates in two-dimensions [40] and itwould be particularly interesting to examine whetherthese, as well as more complicated multi-hump wave-forms emerge systematically from the correspondingmixed state via two-dimensional generalizations of thebifurcations presented herein. Such studies are presentlyin progress and will be presented in future publications. Acknowledgments
The authors acknowledge support from NSF-DMS-0806762. PGK also acknowledges support from NSF-DMS-0349023 (NSF-CAREER) and from the Alexandervon Humboldt Foundation. [1] C. J. Pethick and H. S. Smith,
Bose-Einstein Condensa-tion in Dilute Gases , Cambridge University Press, Cam-bridge, 2002.[2] L.P. Pitaevskii and S. Stringari,
Bose-Einstein Conden-sation , Oxford University Press (Oxford, 2003).[3] P.G. Kevrekidis, D.J. Frantzeskakis, and R. Carretero-Gonz´alez (eds.),
Emergent nonlinear phenomena inBose-Einstein condensates. Theory and experiment (Springer-Verlag, Berlin, 2008).[4] Yu.S. Kivshar and G.P. Agrawal,
Optical solitons: fromfibers to photonic crystals , Academic Press (San Diego,2003).[5] T.-L. Ho and V.B. Shenoy, Phys. Rev. Lett. , 3276(1996); H. Pu and N.P. Bigelow, Phys. Rev. Lett. ,1130 (1998).[6] B.D. Esry et al. , Phys. Rev. Lett. , 3594 (1997).[7] Th. Busch et al. , Phys. Rev. A , 2978 (1997); R.Graham and D. Walls, Phys. Rev. A , 484 (1998); H.Pu and N.P. Bigelow, Phys. Rev. Lett. , 1134 (1998);B.D. Esry and C.H. Greene, Phys. Rev. A , 1265(1998).[8] M. Trippenbach, K. Goral, K. Rzazewski, B. Malomed,and Y.B. Band, J. Phys. B , 4017 (2000).[9] I.M. Merhasin, B. A. Malomed, R. Driben, J. Phys. B:At. Mol. Opt. Phys. , 877-892 (2005).[10] I.M. Merhasin, B.A. Malomed and R. Driben, Phys.Scripta T116 , 18 (2006).[11] S. Coen and M. Haelterman, Phys. Rev. Lett. , 140401(2001).[12] P. ¨Ohberg and L. Santos, Phys. Rev. Lett. , 2918(2001). [13] Th. Busch and J.R. Anglin, Phys. Rev. Lett. , 010401(2001).[14] C. Becker, S. Stellmer, S. Soltan-Panahi, S. Dorcher, M.Baumert, E.-M. Richter, J. Kronjager, K. Bongs and K.Sengstock, Nature Phys. , 496 (2008).[15] K. Kasamatsu, M. Tsubota and M. Ueda, Int. J. Mod.Phys. B , 1835 (2005).[16] B. Deconinck, J.N. Kutz, M.S. Patterson and B.W.Warner, J. Phys. A , 5431 (2003).[17] C.J. Myatt, E.A. Burt, R.W. Ghrist, E.A. Cornell, andC.E. Wieman Phys. Rev. Lett. , 586 (1997).[18] D.S. Hall, M.R. Matthews, J.R. Ensher, C.E. Wieman,and E.A. Cornell, Phys. Rev. Lett. , 1539 (1998).[19] D.M. Stamper-Kurn, M. R. Andrews, A. P. Chikkatur,S. Inouye, H.-J. Miesner, J. Stenger, and W. Ketterle,Phys. Rev. Lett. , 2027 (1998).[20] G. Modugno, G. Ferrari, G. Roati, R.J. Brecha, A. Si-moni, M. Inguscio, Science , 1320 (2001).[21] M. Mudrich, S. Kraft, K. Singer, R. Grimm, A. Mosk,and M. Weidem¨uller, Phys. Rev. Lett. , 253001 (2002).[22] V. Schweikhard, I. Coddington, P. Engels, S. Tung andE.A. Cornell, Phys. Rev. Lett. , 210403 (2004).[23] K.M. Mertes, J.W. Merrill, R. Carretero-Gonz´alez, D.J.Frantzeskakis, P.G. Kevrekidis and D.S. Hall, Phys. Rev.Lett. , 190402 (2007).[24] S.B. Papp, J.M. Pino and C.E. Wieman, Phys. Rev. Lett. , 040402 (2008).[25] K. Kasamatsu, and M. Tsubota, Phys. Rev. A , 013617(2006).[26] G. Karali, P.G. Kevrekidis and N.K. Efremidis, J. Phys.A: Math. Theor , 045206 (2009). [27] R. Carretero-Gonz´alez, D.J. Frantzeskakis, and P.G.Kevrekidis, Nonlinearity , R139 (2008).[28] M.-S. Chang, C.D. Hamley, M.D. Barrett, J.A. Sauer,K.M. Fortier, W. Zhang, L. You, and M.S. Chapman,Phys. Rev. Lett. , 140403 (2004).[29] J. Stenger, S. Inouye, D.M. Stamper-Kurn, H.-J. Mies-ner, A.P. Chikkatur, and W. Ketterle. Nature , 345(1998).[30] H.E. Nistazakis, D.J. Frantzeskakis, P.G. Kevrekidis,B.A. Malomed, R. Carretero-Gonz´alez, and A.R. Bishop,Phys. Rev. A , 063603 (2007).[31] A.E. Leanhardt, Y. Shin, D. Kielpinski, D.E. Pritchard,and W. Ketterle, Phys. Rev. Lett. , 140403 (2003).[32] J. Ieda, T. Miyakawa, and M. Wadati, Phys. Rev. Lett. , 194102 (2004).[33] J. Ieda, T. Miyakawa, and M. Wadati, J. Phys. Soc. Jpn. , 2996 (2004). [34] L. Li, Z. Li, B.A. Malomed, D. Mihalache, and W.M. Liu,Phys. Rev. A , 033611 (2005).[35] W. Zhang, ¨O.E. M¨ustecaplioglu, and L. You, Phys. Rev.A , 043601 (2007).[36] M. Uchiyama, J. Ieda, and M. Wadati, J. Phys. Soc. Jpn. , 064002 (2006).[37] B.J. Dabrowska-W¨uster, E.A. Ostrovskaya, T.J. Alexan-der, and Y.S. Kivshar, Phys. Rev. A , 023617 (2007).[38] H.E. Nistazakis, D.J. Frantzeskakis, P.G. Kevrekidis,B.A. Malomed, and R. Carretero-Gonz´alez, Phys. Rev.A , 033612 (2008).[39] A. Weller, J.P. Ronzheimer, C. Gross, J. Esteve, M.K.Oberthaler, D.J. Frantzeskakis, G. Theocharis, and P.G.Kevrekidis, Phys. Rev. Lett. , 130401 (2008).[40] B.A. Malomed, H.E. Nistazakis, D.J. Frantzeskakis, andP.G. Kevrekidis, Phys. Rev. A70