PHYSICAL REVIEW A 89, 013630 (2014) Demixing and symmetry breaking in binary dipolar Bose-Einstein-condensate solitons S. K. Adhikari* Instituto de Fı́sica Teórica, UNESP-Universidade Estadual Paulista, 01.140-070 São Paulo, São Paulo, Brazil (Received 16 October 2013; published 31 January 2014) We demonstrate a fully demixed (separated) robust and stable bright binary dipolar Bose-Einstein-condensate soliton in a quasi-one-dimensional setting formed due to dipolar interactions for repulsive contact interactions. For large repulsive interspecies contact interaction, the first species may spatially separate from the second species, thus, forming a demixed configuration, which can be spatially symmetric or symmetry broken. In the spatially symmetric case, one of the species occupies the central region, whereas, the other species separates into two equal parts and stays predominantly out of this central region. In the symmetry-broken case, the two species stay side by side. Stability phase diagrams for the binary solitons are obtained. The results are illustrated with realistic values of parameters in the binary 164Dy-168Er and 164Dy-162Dy mixtures. The demixed solitons are really soliton molecules formed of two types of atoms. A proposal for creating dipolar solitons in experiments is also presented. DOI: 10.1103/PhysRevA.89.013630 PACS number(s): 03.75.Hh, 03.75.Mn, 03.75.Kk, 03.75.Lm I. INTRODUCTION A quasi-one-dimensional (quasi-1D) bright soliton and a soliton train were created and were investigated experimentally in Bose-Einstein condensates (BECs) of 7Li [1] and 85Rb atoms [2]. These quasi-1D solitons appear for attractive contact interaction in an axially free BEC under a radial harmonic trap [3]. However, due to collapse instability, in three dimensions (3D), such solitons are fragile and can accommodate only a small number of atoms. The recent observation of the BECs of 164Dy [4,5], 168Er [6], and 52Cr [7,8] atoms with large magnetic dipole moments has opened new directions of research in BEC solitons. In dipolar BECs, in addition to the quasi-1D solitons [9,10], one could also have quasi-two-dimensional (quasi-2D) solitons [11] and vortex solitons [12]. Moreover, dipolar BEC solitons can be stabilized when the harmonic trap(s) is (are) replaced by periodic optical-lattice trap(s) in quasi-1D [13] and quasi- 2D [12] configurations. More interestingly, one can have dipolar BEC solitons for fully repulsive contact interaction [9]. Hence, these solitons stabilized solely by long-range dipolar interaction are less vulnerable to collapse in 3D for a large number of atoms due to the repulsive contact interaction. This creates a new scenario for robust solitons of large number of atoms stabilized by short-range repulsion and long-range inter- and intraspecies dipolar attractions. One can study the interplay between the long-range anisotropic dipolar interactions and the variable contact interactions using a Feshbach resonance [14] in a binary dipolar BEC soliton. In a quasi-1D binary dipolar mixture [15,16] free to move along the polarization z direction, for large repulsive interspecies contact interaction, we demonstrate fully demixed (separated) bright dipolar binary solitons with a minimum of spatial overlap between the two species stabilized by long- range dipolar interactions for repulsive intraspecies contact interactions. The demixed binary soliton could be spatially symmetric with one of the species occupying the central region *adhikari@ift.unesp.br; http://www.ift.unesp.br/users/adhikari (|z| < z0) and the other species staying symmetrically in both peripheral regions (|z| > z0). For the same set of parameters (number of atoms and dipolar and contact interactions), the binary soliton could also stay in a stable but symmetry-broken spatially asymmetric shape where the two species lie side by side along the z direction. The demixed quasi-1D dipolar binary solitons are re- ally soliton molecules. The binding between two (spatial- symmetry-broken) or three (spatially symmetric) solitons formed of two types of dipolar atoms is provided by the long-range interspecies dipolar interaction, which is attractive along the axial polarization direction. Previous consideration of dipolar soliton molecules was limited to two component solitons of the same species of atoms [9]. We consider the numerical solution of a three-dimensional mean-field model [15,16] of bright binary dipolar BEC solitons. We illustrate our findings using realistic values of inter- and intraspecies contact and dipolar interactions in the binary dipolar 164Dy-168Er and 164Dy-162Dy mixtures. The repulsive interspecies contact interaction plays a crucial role for the appearance of the demixed binary dipolar solitons. The domain of demixed solitons is illustrated in phase diagrams involving a critical number of atoms and the interspecies scattering length. The profiles of the binary soliton are displayed in 3D isodensity plots of the two components. The stability of the solitons is demonstrated by considering a time evolution of the binary soliton when the inter- or intraspecies contact interaction is perturbed with a small time-dependent oscillating part. A steady propagation of the binary soliton under such a perturbation over a long period of time establishes the stability. As these solitons are stable and robust, they can be created and can be studied in a laboratory, and we suggest a way to achieve this goal. First, a highly cigar-shaped dipolar BEC with an appropriate number of atoms has to be created. Then, the axial trap has to be removed slowly, when the axially trapped quasi-1D dipolar BEC will turn into an axially free quasi-1D dipolar BEC soliton. The viability of this procedure is demonstrated by real-time simulation of the mean-field model. As the dipolar soliton is an eigenstate of the mean-field model, 1050-2947/2014/89(1)/013630(10) 013630-1 ©2014 American Physical Society http://dx.doi.org/10.1103/PhysRevA.89.013630 S. K. ADHIKARI PHYSICAL REVIEW A 89, 013630 (2014) a very slow turnoff of the axial trap will always lead to a stable soliton. However, a reasonably quick turnoff is found to lead to a stable soliton as demonstrated in this paper. A nondipolar soliton is highly fragile and cannot be created in this fashion. In Sec. II, the time-dependent 3D mean-field model for the binary dipolar BEC soliton is presented. The results of the numerical calculation are shown in Sec. III. The domain of a stable binary soliton is illustrated in stability phase diagrams showing the maximum number of atoms versus interspecies scattering length for fixed values of inter- and intraspecies dipolar interactions and intraspecies scattering lengths. The domain of mixed and demixed binary solitons are identified, and the demixed solitons are further classified by spatially symmetric and symmetry-broken types. The 3D isodensity profiles of the symmetric and symmetry-broken binary solitons are contrasted for the same set of parameters. The viability of preparing these solitons experimentally by relaxing the axial trap on a quasi-1D binary mixture is demonstrated by real-time simulation of the mean-field model. A similar demonstration for a single-component dipolar quasi-1D soliton was also performed. Finally, in Sec. IV, a brief summary of our findings is presented. II. MEAN-FIELD MODEL FOR A BINARY MIXTURE The extension of the mean-field Gross-Pitaevskii (GP) equation to binary dipolar boson-boson [15,16] and boson- fermion [17] mixtures are well established, and, for the sake of completeness, we make a brief summary of the same, appropriate for this study. We consider a bright binary dipolar BEC soliton with the mass, number of atoms, magnetic dipole moment, and scattering length for the two species i = 1,2, given by mi, Ni, μi , and ai , respectively. The intra- (Vi) and interspecies (V12) interactions for two atoms at r and r′ are Vi(R) = μ0μ 2 i 4π 1 − 3 cos2 θ |R|3 + 4π� 2ai mi δ(R), (1) V12(R) = μ0μ1μ2 4π 1 − 3 cos2 θ |R|3 + 2π� 2a12 mR δ(R), (2) respectively, where a12 is the intraspecies scattering length, μ0 is the permeability of free space, θ is the angle made by the vector R with the polarization z direction, R = r − r′, and mR = m1m2/(m1 + m2) is the reduced mass. With these interactions, the coupled GP equations for the binary dipolar BEC can be written as [15,16,18] i� ∂φ1(r,t) ∂t = [ − � 2 2m1 ∇2 + 1 2 m1ω 2 1ρ 2 + 4π� 2 m1 a1N1|φ1(r,t)|2 + 2π� 2 mR a12N2|φ2(r,t)|2 + N1 μ0 μ2 1 4π ∫ Vdd(R)|φ1(r′,t)|2dr′ +N2 μ0μ1μ2 4π ∫ Vdd(R)|φ2(r′,t)|2dr′ ] φ1(r,t), (3) i� ∂φ2(r,t) ∂t = [ − � 2 2m2 ∇2 + 1 2 m2ω 2 2ρ 2 + 4π� 2 m2 a2N2|φ2(r,t)|2 + 2π� 2 mR a12N1|φ1(r,t)|2 + N2 μ0μ 2 2 4π ∫ Vdd(R)|φ2(r′,t)|2dr′ +N1 μ0 μ1μ2 4π ∫ Vdd(R)|φ1(r′,t)|2dr′ ] φ2(r,t), (4) Vdd(R) = 1 − 3 cos2 θ R3 , ρ2 = x2 + y2, i = √−1, (5) with normalization ∫ dr|φi(r,t)|2 = 1. Here ωi’s are the angular frequencies of the traps acting on the two species in the transverse ρ direction. The binary soliton is free to move in the axial z direction. The intra- and interspecies dipolar interactions are usually expressed in terms of the dipolar lengths a (i) dd and a (12) dd , defined by a (i) dd = μ0μ 2 i mi 12π�2 , a (12) dd = μ0μ1μ2mR 6π�2 . (6) We express the strengths of the dipolar interactions by these lengths and transform Eqs. (3) and (4) into the following dimensionless form [15]: i ∂φ1(r,t) ∂t = [ −∇2 2 + 1 2 ρ2 + g1|φ1|2 + g12|φ2|2 + g (1) dd ∫ Vdd(R)|φ1(r′,t)|2dr′ + g (12) dd ∫ Vdd(R)|φ2(r′,t)|2dr′ ] φ1(r,t), (7) i ∂φ2(r,t) ∂t = [ −m12 ∇2 2 + 1 2 mωρ2 + g2|φ2|2 + g21|φ1|2 + g (2) dd ∫ Vdd(R)|φ2(r′,t)|2dr′ + g (21) dd ∫ Vdd(R)|φ1(r′,t)|2dr′ ] φ2(r,t), (8) where mω = ω2 2/(m12ω 2 1), m12 = m1/m2, g1 = 4πa1N1, g2 = 4πa2N2m12, g12 = 2πm1a12N2/mR, g21 = 2πm1a12N1/mR, g (2) dd = 3N2a (2) dd m12, g (1) dd = 3N1a (1) dd , g (12) dd = 3N2a (12) dd m1/(2mR), and g (21) dd = 3N1a (12) dd m1/(2mR). In Eqs. (7) and (8), length is expressed in units of oscillator length l0 = √ �/(m1ω1), energy is expressed in units of oscillator energy �ω1, probability density |φi |2 is expressed in units of l−3 0 , and time is expressed in units of t0 = 1/ω1. 013630-2 DEMIXING AND SYMMETRY BREAKING IN BINARY . . . PHYSICAL REVIEW A 89, 013630 (2014) The dimensionless GP equation for a single-component dipolar quasi-1D soliton is [9] i ∂φ(r,t) ∂t = [ −∇2 2 + 1 2 ρ2 + 4πaN |φ(r,t)|2 + 3addN ∫ Vdd(R)|φ(r′,t)|2dr′ ] φ(r,t), (9) where N is the number of atoms, a is the scattering length, and add is the dipolar length. III. NUMERICAL RESULTS We demonstrate stable demixed spatially symmetric and asymmetric binary dipolar solitons for realistic values of atom numbers and interaction parameters in the 164Dy-168Er and 164Dy-162Dy mixtures. The 164Dy and 168Er atoms have the largest magnetic moments of all the dipolar atoms used in BEC experiments. For the 164Dy-168Er mixture, the 164Dy atoms are labeled i = 1, and the 168Er atoms are labeled i = 2. The magnetic moment of a 164Dy atom is μ1 = 10μB [5] and of a 168Er atom is μ2 = 7μB [6] with μB as the Bohr magneton leading to the dipolar lengths a (1) dd ≡ μ0μ 2 1m1/(12π� 2) ≈ 132.7a0, a (2) dd ≡ μ0μ 2 2m2/(12π� 2) ≈ 66.6a0, and a (12) dd ≡ μ0μ1μ2mR/(6π� 2) ≈ 94.0a0 where a0 is the Bohr radius. For the 164Dy-162Dy mixture, the 164Dy atoms are labeled i = 1, and the 162Dy atoms are labeled i = 2. The magnetic moment of a 162Dy atom is μ2 = 10μB . Consequently, the dipolar lengths are as follows: a (1) dd ≈ 132.7a0, a (2) dd ≈ 131.0a0, and a (12) dd ≈ 131.9a0. The fundamental constants used in the evaluation of these dipolar lengths are as follows: μB = 9.274 01 × 10−24 A m2, a0 = 5.292 × 10−11 m, μ0 = 4π × 10−7 N/A2, � = 1.054 57 × 10−34 m−2 kg−1 s, and 1 amu = 1.660 54 × 10−27 kg. The dipolar interaction in 164Dy atoms is roughly double that in 168Er atoms and about eight times larger than that in 52Cr atoms with a dipolar length add ≈ 15a0 [7]. We consider the trap frequencies of ω1 = ω2 = 2π × 61 Hz so that the length scale is l0 ≡ √ �/m1ω1 ≈ 1 μm and the time scale is t0 ≡ ω−1 1 ≈ 2.6 ms. To the best of our knowledge, the experimental values of none of the scattering lengths considered in this paper are known, except for the intraspecies scattering length of 168Er atoms where preliminary cross-dimensional thermalization measurements point to a scattering length between 150a0 and 200a0 [6]. Nevertheless, scattering lengths can be adjusted experimentally using a Feshbach resonance [14]. We solve the 3D Eqs. (7) and (8) by the split-step Crank-Nicolson discretization scheme using both real- and imaginary-time propagation in 3D Cartesian coordinates in- dependent of the underlying trap symmetry using a space step of 0.1–0.2 and a time step of 0.001–0.005 [19]. The dipolar potential term is treated by Fourier transformation in momentum space using a convolution theorem in the usual fashion [20]. We study demixing in binary solitons supported by long- range dipolar interactions for repulsive interspecies and intraspecies contact interactions. The repulsive contact inter- actions make the dipolar solitons less vulnerable to collapse for a larger number of atoms compared to the nondipolar solitons. 102 103 104 105 0 20 40 60 80 100 120 140 cr it Collapse Stable Stable Collapse Stable U nb ou nd U nb ou nd Dy ( ) 102 103 104 105 0 20 40 60 80 100 120 140 cr it Collapse Stable Stable Collapse Stable U nb ou nd U nb ou nd Dy ( ) Dy ( ) 102 103 104 105 0 20 40 60 80 100 120 140 cr it Collapse Stable Stable Collapse Stable U nb ou nd U nb ou nd Dy ( ) Dy ( ) Er ( ) 102 103 104 105 0 20 40 60 80 100 120 140 cr it Collapse Stable Stable Collapse Stable U nb ou nd U nb ou nd Dy ( ) Dy ( ) Er ( ) Er ( ) 102 103 104 105 0 20 40 60 80 100 120 140 cr it Collapse Stable Stable Collapse Stable U nb ou nd U nb ou nd 102 103 104 105 0 20 40 60 80 100 120 140 cr it Collapse Stable Stable Collapse Stable U nb ou nd U nb ou nd FIG. 1. (Color online) Stability phase diagram illustrating the critical number of atoms in a single-component quasi-1D dipolar BEC soliton of 164Dy or 168Er atoms from numerical calculation (n) and variational approximation (v). The system is repulsive and unbound for a � add. Stable quasi-1D solitons appear for a � add and the number of atoms N below a critical number Ncrit. The oscillator length is l0 = 1 μm. However, to have a net attraction in the system, atoms with large magnetic dipole moments, such as 164Dy and 168Er are considered. In the quasi-1D shape with a radial harmonic trap, the dipolar interaction is attractive in the axial z direction. The net repulsion in the transverse direction is balanced by the external harmonic trap. A. Single-component 164Dy and 168Er solitons To have a feeling about the maximum number of atoms in a single-component quasi-1D dipolar soliton, first we solve Eq. (9) for different values of the scattering length a. We find that, for interaction parameters of 164Dy and 168Er atoms, the solitons are stable up to a critical maximum number of atoms, beyond which the system collapses. In Fig. 1, we plot this critical number Ncrit versus a from numerical simulation using Eq. (9) and from a Gaussian variational approximation to it as developed in Ref. [9]. The variational approximation leads to overbinding and, hence, can accommodate a larger number of atoms in a stable soliton as can be seen in Fig. 1. We find that a stable soliton is possible for a � add and for a number of atoms below this critical number. The critical number of atoms increases with the increase in contact repulsion as a → add, which is counterintuitive. The solitons are bound by long-range dipolar interaction, and an increase in contact repulsion gives more stability against collapse for a fixed dipolar interaction strength. In this phase diagram, three regions are shown: stable, collapse, and unbound. In the unbound region (a � add), contact repulsion dominates over dipolar attraction, and the soliton cannot be bound. In the collapse region, the opposite happens, and the soliton collapses due to excess of attractive dipolar interaction along the axial z direction. In the stable region, there is a balance between attraction and repulsion, and a stable soliton can be formed. So, in our study of binary solitons, we will take a � add. In this fashion, a binary soliton of a large number of atoms could be created, and this would be of greater experimental interest. In this study, we take for the 164Dy atoms, a1 = 120a0, and for the 168Er atoms, we 013630-3 S. K. ADHIKARI PHYSICAL REVIEW A 89, 013630 (2014) FIG. 2. (Color online) Integrated 1D density |φ(z,t)|2 =∫ dx ∫ dy|φ(r,t)|2 of 164Dy atoms during real-time propagation when the axial trap of angular frequency ω = 2π × 2 Hz on a quasi-1D dipolar BEC of 8000 164Dy atoms is removed linearly for 20 > t/t0 > 0 and the resultant axially free soliton is propagated in real time. Parameters used in the simulation: add = 132.7a0, a = 120a0, ωρ = 2π × 61 Hz, and l0 = 1 μm. take a2 = 60a0. From Fig. 1, we see that, for these values of scattering length, a stable bright soliton can accommodate a large number of atoms of the two species. As the dipolar solitons are stable and robust, it is relatively easy to make these solitons of a large number of atoms. The nondipolar solitons are usually tiny and fragile containing a small number of atoms [1,2]. We illustrate a method for creating a dipolar soliton of 8000 164Dy atoms by real-time simulation of the GP equation (9) with the scattering length a = 120a0. By imaginary-time propagation, we create a bound quasi-1D BEC in the trap V (r) = (ρ2 + λ2z2)/2 with λ = 0.03. This corresponds to taking a weak axial trap along the z axis of angular frequency ωz ≈ 2π × 2 Hz compared to the strong radial trap in the x-y plane of angular frequency ωρ = 2π × 61 Hz. The 3D profile of the BEC as obtained in imaginary-time propagation is used as the initial state in the real-time routine. During real-time propagation, from t/t0 = 0 to 20, the axial trap λ2z2/2 gradually (linearly) is reduced to zero so that, for t/t0 > 20, the axially free quasi-1D soliton condition is realized. We continue the real-time propagation for t/t0 > 20 for a reasonably long interval of time. Long sustained propagation of the central peak establishes the soliton nature of the axially free dipolar BEC. The result of this simulation is presented in Fig. 2 where we show the integrated 1D density |φ(z,t)|2 = ∫ dx ∫ dy|φ(r,t)|2 during real-time propagation. Upon relaxation of the axial trap reasonably quickly, the soliton expands a little in the beginning before attaining the final shape as can be seen in Fig. 2. We note that the solitonic nature of the axially free dipolar BEC is evident in Fig. 2. This approach can be used in a laboratory to create dipolar solitons. Fragile nondipolar solitons supported by con- tact attraction only cannot be easily prepared in this fashion. B. Binary 164Dy-168Er soliton Scattering lengths play important roles in the preparation of binary solitons and can be experimentally controlled indepen- dently by magnetic [14] and optical [21] Feshbach resonance techniques. For 164Dy atoms, we take a1 = 120a0, and for 168Er atoms, we take a2 = 60a0. However, the interspecies scattering length a12 plays a crucial role in demixing in binary solitons and will be considered as a variable. After some experimentation in the numerical simulation, we realized that a reasonably large interspecies scattering length a12 (�70a0) and the number of atoms of the species below a critical value are needed for the creation of a demixed binary dipolar soliton. For much larger values of the interspecies scattering length (a12 � 120a0), the repulsion could be sufficiently strong so that solitons of the two species could not be bound to form a stable binary dipolar soliton. For larger values of the number of atoms, the system collapses due to a strong net dipolar interaction even for repulsive intraspecies interaction. The demixing takes place due to interspecies contact repulsion. For smaller values of interspecies scattering length, the repulsion is not sufficient for demixing, and a mixed or overlapping configuration of a binary dipolar BEC soliton is realized [13]. A spatially symmetric binary soliton is obtained by solving Eqs. (7) and (8) with imaginary-time propagation with the fol- lowing spatially symmetric Gaussian input for both functions φi : φi(r) = π−3/4 wρ √ wz exp [ − ρ2 2w2 ρ − z2 2w2 z ] , (10) where w’s are the respective widths. An initial guess of widths close to their final converged values facilitates convergence. The input (10) and the final binary soliton are spatially symmetric around z = 0. To obtain a symmetry-broken binary soliton, a symmetry-broken initial input is needed. Now, we display the typical profile of a spatially symmetric demixed binary dipolar BEC soliton in Fig. 3 for 2000 164Dy atoms and 5000 168Er atoms for the interspecies scattering length a12 = a(Dy-Er) = 90a0 and for large intraspecies scat- tering lengths: a1 = a(Dy) = 120a0 and a2 = a(Er) = 60a0. The species containing the larger number of atoms (168Er) stays at the center, and the species containing the smaller number of FIG. 3. (Color online) Three-dimensional isodensity contour (a) |φ1|2 of 164Dy, (b) |φ2|2 of 168Er, and (c) (|φ1|2 + |φ2|2) of the 164Dy-168Er mixture for the binary soliton of 2000 164Dy and 5000 168Er atoms with a(Dy) = 120a0, a(Er) = 60a0, and a(Dy-Er) = 90a0. The dimensionless lengths x, y, and z are in units of l0(≡ 1 μm). The density on contour is 0.002. 013630-4 DEMIXING AND SYMMETRY BREAKING IN BINARY . . . PHYSICAL REVIEW A 89, 013630 (2014) 102 103 104 105 60 70 80 90 100 110 120 cr it( D y) (Dy-Er)/ 0 (Er) = 1000 collapse 102 103 104 105 60 70 80 90 100 110 120 cr it( D y) (Dy-Er)/ 0 (Er) = 1000 collapse 102 103 104 105 60 70 80 90 100 110 120 cr it( D y) (Dy-Er)/ 0 (Er) = 1000 collapse 102 103 104 105 60 70 80 90 100 110 120 cr it( D y) (Dy-Er)/ 0 (Er) = 1000 collapse demixed demixed mixed (a) variational 102 103 104 105 60 70 80 90 100 110 120 cr it( D y) (Dy-Er)/ 0 (Er) = 5000 collapse 102 103 104 105 60 70 80 90 100 110 120 cr it( D y) (Dy-Er)/ 0 (Er) = 5000 collapse 102 103 104 105 60 70 80 90 100 110 120 cr it( D y) (Dy-Er)/ 0 (Er) = 5000 collapse 102 103 104 105 60 70 80 90 100 110 120 cr it( D y) (Dy-Er)/ 0 (Er) = 5000 collapse demixed demixedmixed(b) variational FIG. 4. (Color online) Stability phase diagram for a binary 164Dy-168Er soliton for (a) N (Er) = 1000 and (b) N (Er) = 5000 atoms. The variational mixed-collapse border is shown by solid circles. The parameters used in the simulation are l0 = 1 μm, a(Dy) = 120a0, and a(Er) = 60a0. A stable demixed binary soliton is observed in the darker (pink and green) regions to the right, whereas, a stable mixed binary soliton appears in the lighter (yellow) region to the left. For a spatially symmetric demixed binary soliton, in the upper darker (pink) region, the 168Er atoms stay out, and in the lower darker (green) region, the 164Dy atoms stay out of the central region (z = 0). atoms (164Dy) breaks into two equal parts, leaves the central region occupied by the 168Er atoms, and stays symmetrically on two sides of the 168Er soliton with a minimum of interspecies overlap. Note that the “3D densities” plotted in this paper are the norms of the respective wave functions |φi |2, whereas, the atom densities are Ni |φi |2. Next, we show the stability phase diagrams for the appear- ance of the demixed binary dipolar soliton in the 164Dy-168Er mixture for different numbers of atoms and variable inter- species scattering length a12. In Figs. 4(a) and 4(b), we show the critical number of 164Dy atoms N1 = Ncrit(Dy) versus a12 for the numbers of 168Er atoms N2 = N (Er) = 1000,5000, respectively. The binary system collapses for N (Dy) larger than these critical values. For small values of a12, the binary soliton is mixed, and it gets demixed for large values of a12 due to strong interspecies contact repulsion. A Gaussian variational approximation to Eqs. (7) and (8) was developed in Ref. [16] for a mixed binary soliton. This approximation, using a mixed configuration of the two solitons, can be used to calculate the mixed-collapse boundary. The variational results for this boundary are also shown in Fig. 4. As the variational approximation overbinds the binary soliton, it shows a larger stable region than the full numerical solution. For a12 � 80a0, the binary solitons are spatially symmetric and mixed, and the two species of atoms lie on top of each other. However, for larger a12’s, the contact repulsion between the two species of atoms increases, and the system gets demixed. There are two 102 103 104 105 60 70 80 90 100 110 120 cr it( E r) (Dy-Er)/ 0 (Dy) = 1000 collapse 102 103 104 105 60 70 80 90 100 110 120 cr it( E r) (Dy-Er)/ 0 (Dy) = 1000 collapse 102 103 104 105 60 70 80 90 100 110 120 cr it( E r) (Dy-Er)/ 0 (Dy) = 1000 collapse 102 103 104 105 60 70 80 90 100 110 120 cr it( E r) (Dy-Er)/ 0 (Dy) = 1000 collapse demixed demixed mixed (a) variational 102 103 104 105 60 70 80 90 100 110 120 (Dy) = 5000 collapse 102 103 104 105 60 70 80 90 100 110 120 (Dy) = 5000 collapse 102 103 104 105 60 70 80 90 100 110 120 (Dy) = 5000 collapse 102 103 104 105 60 70 80 90 100 110 120 (Dy) = 5000 collapse 102 103 104 105 60 70 80 90 100 110 120 (Dy) = 5000 collapse 102 103 104 105 60 70 80 90 100 110 120 (Dy) = 5000 collapse demixed mixed demixed (b) variational cr it( E r) (Dy-Er)/ 0 cr it( E r) (Dy-Er)/ 0 cr it( E r) (Dy-Er)/ 0 cr it( E r) (Dy-Er)/ 0 FIG. 5. (Color online) The same as in Fig. 4 for (a) N (Dy) = 1000 and (b) N (Dy) = 5000 atoms. ways this demixing takes place. The demixed binary soliton can either be spatially symmetric around z = 0, or it can be spatially asymmetric around z = 0. The two configurations are possible for identical parameters of the binary system, e.g., number of atoms of the two components and all dipolar and contact interaction strengths. However, there is a region around N1 ≈ N2 for these larger values of a12 where the solitons continue in a mixed configuration. On both sides of this region, demixing can take place. In the spatially symmetric configuration, one of the components (the one containing the larger number of atoms) continues in the spatially symmetric state at z = 0. The other component (the one containing the smaller number of atoms) divides into two equal pieces, separates from each other, leaves the central region occupied by the first component around z = 0, moves in opposite directions, and stays for |z| > z0 in a spatially symmetric configuration as shown in Fig. 3. In the spatial-symmetry- broken configuration, both components remain as single pieces and move away from each other and eventually lie on both sides of z = 0 in a spatially asymmetric configuration. In Figs. 5(a) and 5(b), we show the critical number of 168Er atoms N2 = Ncrit(Er) versus a12 for the numbers of 164Dy atoms N1 = N (Dy) = 1000,5000, respectively. A similar scenario emerges in Fig. 5 as in Fig. 4. In Fig. 3, we presented a spatially symmetric demixed binary dipolar soliton with 2000 164Dy and 5000 168Er atoms. In this case, the 168Er atoms lie in the central region, whereas, 164Dy atoms move away from the central region. The situation changes with a larger number of 164Dy atoms. This is illustrated Fig. 6 with 5000 164Dy atoms and 1000 168Er atoms. In this case, in the isodensity contour of the binary soliton, the roles of the 164Dy and 168Er atoms are changed compared to that in Fig. 3. Now, contrary to that in Fig. 3, the 164Dy atoms occupy the central region and the 168Er atoms occupy the peripheral region. In Figs. 4 and 5 in the darker pink region in the spatially 013630-5 S. K. ADHIKARI PHYSICAL REVIEW A 89, 013630 (2014) FIG. 6. (Color online) Three-dimensional isodensity contour (a) |φ2|2 of 168Er, (b) |φ1|2 of 164Dy, and (c) (|φ1|2 + |φ2|2) of the 164Dy-168Er mixture for the binary soliton of 5000 164Dy and 1000 168Er atoms with a(Dy) = 120a0, a(Er) = 60a0, and a(Dy-Er) = 90a0. The dimensionless lengths x, y, and z are in units of l0(≡ 1 μm). The density on the contour is 0.002. symmetric configuration, the 164Dy atoms stay at the center and the 168Er atoms stay in the outer region, and the opposite happens in the darker green region. The binary soliton of Fig. 3 lies in the darker green area of Fig. 4(b), and the binary soliton of Fig. 6 lies in the darker pink area of Fig. 5(b). It is important to establish the stability of these spatially symmetric solitons under a small variation in the parameters of the initial state. To achieve this, we consider the final state of the spatially symmetric binary soliton as obtained from the imaginary-time routine and use it as the initial state in the real-time program. First, we perform the numerical simulation with the binary soliton of Fig. 3 with the time-dependent interspecies scattering length a12 = 90a0 + 10a0 cos(t/t0), whereas, the initial state was obtained with a12 = 90a0. All other contact and dipolar interaction strengths are maintained at their initial values. This perturbation should be sufficient to destroy an unstable soliton. The system executes small oscillation around the stable position. Otherwise, no noticeable change is observed after 100 units of real-time propagation. Our results are illustrated in Figs. 7 and 8. In Fig. 7, we illustrate the contour plot of integrated 1D density |φi(z,t)|2 =∫ ∫ dx dy|φ(x,y,z,t)|2 during real-time propagation. The vi- sual profile of the solitons in Fig. 7 maintains a constant width and does not show a sign of instability. To have a more quantitative view, in Fig. 8 , we plot the initial and final 1D densities |φi(z,t)|2 at times t/t0 = 0 and 100 during real-time propagation. This clearly demonstrates the stability of the binary dipolar soliton under small oscillation. The binary dipolar solitons studied above were spatially symmetric. One can also have symmetry-broken spatially asymmetric binary solitons for the same parameters used as in the case of spatially symmetric binary solitons. However, to generate these solitons by imaginary-time propagation, one should consider spatially asymmetric initial states. In particular, we consider the following spatial-symmetry-broken FIG. 7. (Color online) Contour plot of integrated 1D densities |φi(z,t)|2 = ∫ dx ∫ dy|φ(r,t)|2 of 164Dy and 168Er atoms when the binary soliton of Fig. 3 is propagated in real time with the time-dependent interspecies scattering length a(Dy-Er) = 90a0 + 10a0 cos(t/t0) with all other parameters maintained at their time- independent initial values. (asymmetric around z = 0) initial states with a spatial displace- ment between the two profiles in imaginary-time propagation: φ1(r) = π−3/4 wρ √ wz exp [ − ρ2 2w2 ρ − (z − z0)2 2w2 z ] , (11) φ2(r) = π−3/4 wρ √ wz exp [ − ρ2 2w2 ρ − (z + z0)2 2w2 z ] , (12) in place of the initial states (10). It is important to mention that the imaginary-time propagation converges to the lowest- energy state maintaining the symmetry of the initial state. For example, in the 1D harmonic-oscillator problem, the imaginary-time propagation with a spatially symmetric initial state will converge to the ground state of the problem, whereas, with a spatially antisymmetric initial state, it will converge to the spatially antisymmetric first excited state. With the initial states (11) and (12), we obtain the symmetry-broken spatially asymmetric binary solitons with the same parameters as used in the spatially symmetric binary solitons. We tried some other forms of the symmetry-broken initial state and found that the solution always converges to the same final state. After some experimentation, for the same set of parameters of the problem, Dy (Fi) Er (Fi) Dy (In) Er (In) FIG. 8. (Color online) Initial (In) at t = 0 and final (Fi) at t = 100 values of integrated 1D densities |φi(z)|2 = ∫ dx ∫ dy|φ(r)|2 of 164Dy and 168Er atoms in the binary 164Dy-168Er soliton of Fig. 3 during real-time propagation of Fig. 7. 013630-6 DEMIXING AND SYMMETRY BREAKING IN BINARY . . . PHYSICAL REVIEW A 89, 013630 (2014) FIG. 9. (Color online) The same as in Fig. 3 for the symmetry- broken binary soliton of 2000 164Dy and 5000 168Er atoms. we could locate only two definite demixed final states: (a) the symmetric states of Figs. 3 and 6 and (b) the symmetry-broken states to be presented next. We perform imaginary-time propagation with the initial states (11) and (12) using the same parameters as in Fig. 3. The result was insensitive to the value of the parameter z0 in Eqs. (11) and (12), and we take z0 = 10 in this study. The 3D isodensity contours of the resultant soliton profiles in this case for 2000 164Dy and 5000 168Er atoms are shown in Fig. 9. The converged binary dipolar soliton is clearly spatially asymmetric. The second species with a larger number of 168Er atoms has a longer spatial extension than the 164Dy atoms. Then, we consider imaginary-time propagation using the same parameters as in Fig. 6. The 3D isodensity contours of the converged state are shown in Fig. 10. In this case with 5000 164Dy and 1000 168Er atoms, the 164Dy BEC clearly occupies a larger region of space than the 168Er atoms. FIG. 10. (Color online) The same as in Fig. 3 for the symmetry- broken binary soliton of 5000 164Dy and 1000 168Er atoms. FIG. 11. (Color online) The same as in Fig. 7 for the symmetry- broken binary soliton of Fig. 10. In the case of the demixed solitons, it is interesting to ask which of the two configurations—symmetric or asymmetric— has the lower energy and, hence, will be favorable experi- mentally. To this end, we calculated the respective energies numerically in several cases and consistently found that the spatially asymmetric configuration has slightly lower energy. However, this will not be of phenomenological interest as the energy difference is very small, being on the order of 0.1% of the total energy of either state. Now, we investigate the stability of these symmetry- broken binary dipolar solitons. To this end, we consider the converged solutions of the imaginary-time routine and use these as the initial states in the real-time routine with the oscillating interspecies scattering length a12 = a(Dy-Er) = 90a0 + 10a0 cos(t/t0) as in the case of oscillation presented in Fig. 7. During real-time propagation, the system executes small oscillation, and the contour plot of the integrated 1D densities is illustrated in Fig. 11. There is no visible change in the widths of the condensates during real-time propagation, confirming stability. We also compared the 1D and 3D profiles of the initial and final states (not presented here) and found that they are practically identical. The demixed binary solitons, as presented in Figs. 3, 6, 9, and 10, are ideal examples of soliton molecules of different species where the participating solitons maintain their identities with a minimum of overlap. The binding between the two solitons comes from the long-range dipolar interaction. In nondipolar systems, such binding can only come from short-range interspecies attraction, and a composite soliton molecule can possibly appear only with complete overlap between participating solitons [16,22]. The demixed binary dipolar soliton or the binary soliton molecule can be prepared and can be observed experimentally. We illustrate the possibility of experimental realization of the soliton shown in Figs. 6 and 9 by realistic numerical simulation using the binary GP Eqs. (7) and (8). In case of the spatially symmetric binary soliton of Fig. 6, we prepare a spatially symmetric binary quasi-1D trapped dipolar BEC of 5000 164Dy and 1000 168Er atoms with the same parameters as in Fig. 6 by imaginary-time simulation but with traps V (r) = (ρ2/2 + λ2z2/2) and V (r) = mω(ρ2/2 + λ2z2/2) with λ = 0.03, respectively, on 164Dy and 168Er atoms in place of the axially free traps V (r) = ρ2/2 and V (r) = mωρ2/2. For both 164Dy and 168Er atoms, these traps correspond to one of 013630-7 S. K. ADHIKARI PHYSICAL REVIEW A 89, 013630 (2014) FIG. 12. (Color online) Integrated 1D density |φi(z,t)|2 =∫ dx ∫ dy|φi(r,t)|2 of 164Dy (red) and 168Er (green) atoms during real-time propagation when the axial trap of angular frequency ωz = 2π × 2 Hz on a (a) quasi-1D spatially symmetric binary dipolar BEC of 5000 164Dy and 1000 168Er atoms and a (b) quasi-1D spatially asymmetric binary dipolar BEC of 2000 164Dy and 5000 168Er atoms are removed linearly for 20 > t/t0 > 0 and the resultant soliton is propagated in real time. Parameters used in the simulation: a(Dy) = 120a0, a(Er) = 60a0, a(Dy-Er) = 90a0, ωρ = 2π × 61 Hz, and l0 = 1 μm. axial angular frequency ωz = 2π × 2 Hz and of radial angular frequency ωρ = 2π × 61 Hz. This binary bound dipolar BEC is then used as the initial state in the real-time program. During real-time propagation, from t/t0 = 0 to 20, the axial traps λ2z2/2 are gradually (linearly) reduced to zero so that, for t/t0 > 20, the axially free quasi-1D binary dipolar soliton emerges. The simulation is continued for a long interval of time, and a steady propagation of the binary demixed dipolar soliton is established. The outcome of this simulation is shown in Fig. 12(a) where we plot the integrated 1D density |φi(z,t)|2 of the two component solitons as obtained from real-time propagation. A similar simulation for the spatially asymmetric binary soliton of Fig. 9 is shown in Fig. 12(b). After being released from the axial trap for t/t0 > 20, the solitons expand a little and then oscillate around their appropriate sizes. The stability of the demixed binary dipolar solitons is evident in the simulations illustrated in Figs. 12(a) and 12(b). 102 103 104 105 60 70 80 90 100 110 120 130 140 cr it( 16 4 D y) 162Dy-164Dy)/ 0 (162Dy) = 1000 collapse 102 103 104 105 60 70 80 90 100 110 120 130 140 cr it( 16 4 D y) 162Dy-164Dy)/ 0 (162Dy) = 1000 collapse 102 103 104 105 60 70 80 90 100 110 120 130 140 cr it 16 4 D y) 162Dy-164Dy)/ 0 (162Dy) = 1000 collapse 102 103 104 105 60 70 80 90 100 110 120 130 140 cr it 16 4 D y) 162Dy-164Dy)/ 0 (162Dy) = 1000 collapse demixed demixed mixed variational FIG. 13. (Color online) The same as in Fig. 4 for a binary 164Dy-162Dy soliton for N (162Dy) = 1000 atoms. The variational mixed-collapse boundary is indicated by solid circles. For a spatially symmetric demixed binary soliton, in the upper darker (pink) region, the 162Dy atoms stay out, and in the lower darker (green) region, the 164Dy atoms stay out of the central region. Parame- ters used in the simulation: l0 = 1 μm and a(164Dy) = a(164Dy) = 120a0. C. Binary 164Dy-162Dy soliton Here, we consider yet another type of binary dipolar soliton, e.g., the 164Dy-162Dy soliton. This is particularly interesting as Lev and his collaborators are studying this binary mixture in the laboratory at the Stanford University [23]. In order to permit a large number of atoms in the solitons, we consider a large value for the scattering lengths, e.g., a(162Dy) = a(164Dy) = 120a0, viz. Fig. 1. The interspecies scattering length is considered as a variable. First, we consider the stability phase plot for this binary soliton. In Fig. 13, we show the critical number of 164Dy atoms in the binary soliton for 1000 162Dy atoms. Again, mixed and demixed binary solitons appear. As the mass and dipolar lengths are almost the same for the two isotopes, the binary plot is quasisymmetric. The plot for 1000 164Dy atoms in the binary soliton will lead to a practically identical plot as that in Fig. 13 with the role of the two isotopes interchanged. This system is highly dipolar compared to the 164Dy-168Er mixture and, hence, can accommodate a smaller number of solitons as can be seen by comparing Fig. 13 with Figs. 4 and 5 with 1000 atoms of one of the components. Compared to the 164Dy-168Er mixture, demixed binary solitons in 164Dy-162Dy only appear for a larger value of interspecies scattering length to compensate for larger interspecies dipolar attraction. In Figs. 4 and 5, demixing appears for interspecies scattering length a12 � 100a0, whereas, in Fig. 13, demixing appears for a12 � 120a0. In this case, one also has mixed, spatially symmetric demixed, and spatial symmetry-broken demixed binary soli- tons. Without repeating a detailed discussion of different types of solitons, in Fig. 14, we show the isodensity contour of a typical spatially symmetric binary 164Dy-162Dy soliton for 1000 162Dy atoms and 5000 164Dy atoms for the interspecies scattering length a(162Dy-164Dy) = 130a0. The component 162Dy, with a smaller number of atoms, stays out of the central region, whereas, the component 164Dy, with a larger number of atoms, stays at the center. 013630-8 DEMIXING AND SYMMETRY BREAKING IN BINARY . . . PHYSICAL REVIEW A 89, 013630 (2014) FIG. 14. (Color online) Three-dimensional isodensity contour (a) |φ1|2 of 162Dy, (b) |φ2|2 of 164Dy, and (c) (|φ1|2 + |φ2|2) of the 164Dy-162Dy mixture for the binary soliton of 1000 162Dy and 5000 164Dy atoms with a(162Dy) = a(164Dy) = 120a0 and a(162Dy-164Dy) = 130a0. The dimensionless lengths x, y, and z are in units of l0(≡1 μm). The density on contour is 0.002. IV. SUMMARY AND DISCUSSION Using the numerical solution of a set of coupled 3D mean-field GP equations, we demonstrate the existence of demixed dipolar binary 164Dy-168Er and 164Dy-162Dy solitons stabilized by inter- and intraspecies dipolar interactions in the presence of repulsive inter- and intraspecies contact interactions. The domain of the appearance of the binary soliton is highlighted in stability diagrams of a number of atoms in the two components and interspecies scattering length a12 for fixed dipolar and intraspecies contact interactions. The binary soliton is stable for a maximum number of atoms beyond which it collapses. For small interspecies interaction a12, the binary soliton is mixed, and for large a12, it is demixed. The mixed-collapse boundary was also calculated using a Gaussian variational approximation [16] and was found to be in reasonable agreement with the numerical result. We found two types of demixing in this case: (a) spatially symmetric and (b) spatial-symmetry-broken demixed binary solitons. In the spatially symmetric case, the species with the smaller number of atoms breaks up into two equal parts, and the parts leave the central region occupied by the species with a larger number of atoms and consolidate on both sides of the BEC component with the larger number of atoms as shown in Figs. 3, 6, and 14. In the spatially asymmetric case, the two solitons move away from each other and finally stabilize side by side as shown in Figs. 9 and 10. The two species of the demixed soliton have a minimum of overlap and are stabilized by dipolar interaction in the presence of reasonably strong contact interaction. Such demixed solitons are not possible in the absence of dipolar interaction. We also tested the stability of the demixed solitons by considering real-time propagation with oscillating scattering length with small amplitude. A stable binary soliton should then execute small oscillations. The converged solution of the binary soliton in the imaginary-time propagation is used as the initial state of real-time propagation. The solitons are found to exhibit small oscillation over a long time, thus, establishing their stability. The solitons considered in this paper are stabilized by long- range dipolar attraction and short-range contact repulsion. The dipolar interaction is attractive in the axial direction, and the dipolar repulsion in the transverse direction is compensated by a harmonic trap. Hence, unlike normal BEC solitons stabilized by short-range contact attraction, the present dipolar BECs will be more immune to collapse due to short-range repulsion and can easily accommodate 10 000 atoms of the binary 164Dy-168Er mixture as can be seen from Figs. 4, 5, and 13. The dipolar (Nadd) and contact (Na) nonlinear strengths for such solitons could be typically ∼50 for 5000 164Dy atoms. A nondipolar soliton is created only by contact attraction without any repulsion and is fragile due to collapse instability. In a stable nondipolar soliton, the maximum contact interaction strength is N |a| ∼ 0.6 [1–3]. Thus, the dipolar soliton with larger interaction strengths can accommodate a much larger number of atoms than the nondipolar solitons and should be of great experimental interest. A way of experimentally realizing these demixed dipolar solitons is suggested, and the viability of this scheme is demonstrated by real-time simulation. First, a quasi-1D demixed binary BEC is to be formed with a weak trap in the axial z direction. Then, the axial trap is to be removed linearly in a reasonably short time, when the demixed binary BEC turns into a demixed binary soliton. With the present experimental techniques, such binary dipolar BECs can be observed, and the conclusions of this study can be verified. ACKNOWLEDGMENTS We thank FAPESP and CNPq (Brazil) for partial support. [1] K. E. Strecker, G. B. Partridge, A. G. Truscott, and R. G. Hulet, Nature (London) 417, 150 (2002); L. Khaykovich, F. Schreck, G. Ferrari, T. Bourdel, J. Cubizolles, L. D. Carr, Y. Castin, and C. Salomon, Science 256, 1290 (2002). [2] S. L. Cornish, S. T. Thompson, and C. E. Wieman, Phys. Rev. Lett. 96, 170401 (2006). [3] V. M. Perez-Garcia, H. Michinel, and H. Herrero, Phys. Rev. A 57, 3837 (1998). [4] M. Lu, S. H. Youn, and B. L. Lev, Phys. Rev. Lett. 104, 063001 (2010); J. J. McClelland and J. L. Hanssen, ibid. 96, 143005 (2006); S. H. Youn, M. W. Lu, U. Ray, and B. V. Lev, Phys. Rev. A 82, 043425 (2010). [5] M. Lu, N. Q. Burdick, S. H. Youn, and B. L. Lev, Phys. Rev. Lett. 107, 190401 (2011). [6] K. Aikawa, A. Frisch, M. Mark, S. Baier, A. Rietzler, R. Grimm, and F. Ferlaino, Phys. Rev. Lett. 108, 210401 (2012). [7] T. Lahaye et al., Nature (London) 448, 672 (2007); A. Griesmaier, J. Stuhler, T. Koch, M. Fattori, T. Pfau, and S. Giovanazzi, Phys. Rev. Lett. 97, 250402 (2006). 013630-9 http://dx.doi.org/10.1038/nature747 http://dx.doi.org/10.1038/nature747 http://dx.doi.org/10.1038/nature747 http://dx.doi.org/10.1038/nature747 http://dx.doi.org/10.1126/science.1071021 http://dx.doi.org/10.1126/science.1071021 http://dx.doi.org/10.1126/science.1071021 http://dx.doi.org/10.1126/science.1071021 http://dx.doi.org/10.1103/PhysRevLett.96.170401 http://dx.doi.org/10.1103/PhysRevLett.96.170401 http://dx.doi.org/10.1103/PhysRevLett.96.170401 http://dx.doi.org/10.1103/PhysRevLett.96.170401 http://dx.doi.org/10.1103/PhysRevA.57.3837 http://dx.doi.org/10.1103/PhysRevA.57.3837 http://dx.doi.org/10.1103/PhysRevA.57.3837 http://dx.doi.org/10.1103/PhysRevA.57.3837 http://dx.doi.org/10.1103/PhysRevLett.104.063001 http://dx.doi.org/10.1103/PhysRevLett.104.063001 http://dx.doi.org/10.1103/PhysRevLett.104.063001 http://dx.doi.org/10.1103/PhysRevLett.104.063001 http://dx.doi.org/10.1103/PhysRevLett.96.143005 http://dx.doi.org/10.1103/PhysRevLett.96.143005 http://dx.doi.org/10.1103/PhysRevLett.96.143005 http://dx.doi.org/10.1103/PhysRevLett.96.143005 http://dx.doi.org/10.1103/PhysRevA.82.043425 http://dx.doi.org/10.1103/PhysRevA.82.043425 http://dx.doi.org/10.1103/PhysRevA.82.043425 http://dx.doi.org/10.1103/PhysRevA.82.043425 http://dx.doi.org/10.1103/PhysRevLett.107.190401 http://dx.doi.org/10.1103/PhysRevLett.107.190401 http://dx.doi.org/10.1103/PhysRevLett.107.190401 http://dx.doi.org/10.1103/PhysRevLett.107.190401 http://dx.doi.org/10.1103/PhysRevLett.108.210401 http://dx.doi.org/10.1103/PhysRevLett.108.210401 http://dx.doi.org/10.1103/PhysRevLett.108.210401 http://dx.doi.org/10.1103/PhysRevLett.108.210401 http://dx.doi.org/10.1038/nature06036 http://dx.doi.org/10.1038/nature06036 http://dx.doi.org/10.1038/nature06036 http://dx.doi.org/10.1038/nature06036 http://dx.doi.org/10.1103/PhysRevLett.97.250402 http://dx.doi.org/10.1103/PhysRevLett.97.250402 http://dx.doi.org/10.1103/PhysRevLett.97.250402 http://dx.doi.org/10.1103/PhysRevLett.97.250402 S. K. ADHIKARI PHYSICAL REVIEW A 89, 013630 (2014) [8] J. Stuhler, A. Griesmaier, T. Koch, M. Fattori, T. Pfau, S. Giovanazzi, P. Pedri, and L. Santos, Phys. Rev. Lett. 95, 150406 (2005); K. Goral, K. Rzazewski, and T. Pfau, Phys. Rev. A 61, 051601 (2000); T. Koch et al., Nat. Phys. 4, 218 (2008); T. Lahaye et al., Rep. Prog. Phys. 72, 126401 (2009). [9] L. E. Young-S, P. Muruganandam, and S. K. Adhikari, J. Phys. B 44, 101001 (2011). [10] S. Giovanazzi, D. O’Dell, and G. Kurizki, J. Phys. B 34, 4757 (2001). [11] R. Nath, P. Pedri, and L. Santos, Phys. Rev. Lett. 102, 050401 (2009); P. Pedri and L. Santos, ibid. 95, 200404 (2005); I. I. Tikhonenkov, B. A. Malomed, and A. Vardi, ibid. 100, 090406 (2008); P. Köberle, D. Zajec, G. Wunner, and B. A. Malomed, Phys. Rev. A 85, 023630 (2012); R. Eichler, D. Zajec, P. Köberle, J. Main, and G. Wunner, ibid. 86, 053611 (2012); R. Eichler, J. Main, and G. Wunner, ibid. 83, 053604 (2011); A.-X. Zhang and J.-K. Xue, ibid. 82, 013606 (2010); V. M. Lashkin, ibid. 75, 043607 (2007). [12] S. K. Adhikari and P. Muruganandam, J. Phys. B 45, 045301 (2012); I. Tikhonenkov, B. A. Malomed, and A. Vardi, Phys. Rev. A 78, 043614 (2008). [13] S. K. Adhikari and P. Muruganandam, Phys. Lett. A 376, 2200 (2012). [14] S. Inouye et al., Nature (London) 392, 151 (1998). [15] Luis E. Young-S. and S. K. Adhikari, Phys. Rev. A 86, 063611 (2012); 87, 013618 (2013). [16] S. K. Adhikari and Luis E. Young-S., J. Phys. B 47, 015302 (2014). [17] S. K. Adhikari, Phys. Rev. A 88, 043603 (2013). [18] R. M. Wilson, C. Ticknor, J. L. Bohn, and E. Timmermans, Phys. Rev. A 86, 033606 (2012); H. Saito, Y. Kawaguchi, and M. Ueda, Phys. Rev. Lett. 102, 230403 (2009). [19] P. Muruganandam and S. K. Adhikari, Comput. Phys. Commun. 180, 1888 (2009); D. Vudragovic, I. Vidanovic, A. Balaz, P. Muruganandam, and S. K. Adhikari, ibid. 183, 2021 (2012). [20] K. Goral and L. Santos, Phys. Rev. A 66, 023613 (2002); S. Yi and L. You, ibid. 63, 053607 (2001). [21] S. Blatt, T. L. Nicholson, B. J. Bloom, J. R. Williams, J. W. Thomsen, P. S. Julienne, and J. Ye, Phys. Rev. Lett. 107, 073202 (2011). [22] S. K. Adhikari, Phys. Lett. A 346, 179 (2005). [23] B. L. Lev (private communication). 013630-10 http://dx.doi.org/10.1103/PhysRevLett.95.150406 http://dx.doi.org/10.1103/PhysRevLett.95.150406 http://dx.doi.org/10.1103/PhysRevLett.95.150406 http://dx.doi.org/10.1103/PhysRevLett.95.150406 http://dx.doi.org/10.1103/PhysRevA.61.051601 http://dx.doi.org/10.1103/PhysRevA.61.051601 http://dx.doi.org/10.1103/PhysRevA.61.051601 http://dx.doi.org/10.1103/PhysRevA.61.051601 http://dx.doi.org/10.1038/nphys887 http://dx.doi.org/10.1038/nphys887 http://dx.doi.org/10.1038/nphys887 http://dx.doi.org/10.1038/nphys887 http://dx.doi.org/10.1088/0034-4885/72/12/126401 http://dx.doi.org/10.1088/0034-4885/72/12/126401 http://dx.doi.org/10.1088/0034-4885/72/12/126401 http://dx.doi.org/10.1088/0034-4885/72/12/126401 http://dx.doi.org/10.1088/0953-4075/44/10/101001 http://dx.doi.org/10.1088/0953-4075/44/10/101001 http://dx.doi.org/10.1088/0953-4075/44/10/101001 http://dx.doi.org/10.1088/0953-4075/44/10/101001 http://dx.doi.org/10.1088/0953-4075/34/23/318 http://dx.doi.org/10.1088/0953-4075/34/23/318 http://dx.doi.org/10.1088/0953-4075/34/23/318 http://dx.doi.org/10.1088/0953-4075/34/23/318 http://dx.doi.org/10.1103/PhysRevLett.102.050401 http://dx.doi.org/10.1103/PhysRevLett.102.050401 http://dx.doi.org/10.1103/PhysRevLett.102.050401 http://dx.doi.org/10.1103/PhysRevLett.102.050401 http://dx.doi.org/10.1103/PhysRevLett.95.200404 http://dx.doi.org/10.1103/PhysRevLett.95.200404 http://dx.doi.org/10.1103/PhysRevLett.95.200404 http://dx.doi.org/10.1103/PhysRevLett.95.200404 http://dx.doi.org/10.1103/PhysRevLett.100.090406 http://dx.doi.org/10.1103/PhysRevLett.100.090406 http://dx.doi.org/10.1103/PhysRevLett.100.090406 http://dx.doi.org/10.1103/PhysRevLett.100.090406 http://dx.doi.org/10.1103/PhysRevA.85.023630 http://dx.doi.org/10.1103/PhysRevA.85.023630 http://dx.doi.org/10.1103/PhysRevA.85.023630 http://dx.doi.org/10.1103/PhysRevA.85.023630 http://dx.doi.org/10.1103/PhysRevA.86.053611 http://dx.doi.org/10.1103/PhysRevA.86.053611 http://dx.doi.org/10.1103/PhysRevA.86.053611 http://dx.doi.org/10.1103/PhysRevA.86.053611 http://dx.doi.org/10.1103/PhysRevA.83.053604 http://dx.doi.org/10.1103/PhysRevA.83.053604 http://dx.doi.org/10.1103/PhysRevA.83.053604 http://dx.doi.org/10.1103/PhysRevA.83.053604 http://dx.doi.org/10.1103/PhysRevA.82.013606 http://dx.doi.org/10.1103/PhysRevA.82.013606 http://dx.doi.org/10.1103/PhysRevA.82.013606 http://dx.doi.org/10.1103/PhysRevA.82.013606 http://dx.doi.org/10.1103/PhysRevA.75.043607 http://dx.doi.org/10.1103/PhysRevA.75.043607 http://dx.doi.org/10.1103/PhysRevA.75.043607 http://dx.doi.org/10.1103/PhysRevA.75.043607 http://dx.doi.org/10.1088/0953-4075/45/4/045301 http://dx.doi.org/10.1088/0953-4075/45/4/045301 http://dx.doi.org/10.1088/0953-4075/45/4/045301 http://dx.doi.org/10.1088/0953-4075/45/4/045301 http://dx.doi.org/10.1103/PhysRevA.78.043614 http://dx.doi.org/10.1103/PhysRevA.78.043614 http://dx.doi.org/10.1103/PhysRevA.78.043614 http://dx.doi.org/10.1103/PhysRevA.78.043614 http://dx.doi.org/10.1016/j.physleta.2012.05.030 http://dx.doi.org/10.1016/j.physleta.2012.05.030 http://dx.doi.org/10.1016/j.physleta.2012.05.030 http://dx.doi.org/10.1016/j.physleta.2012.05.030 http://dx.doi.org/10.1038/32354 http://dx.doi.org/10.1038/32354 http://dx.doi.org/10.1038/32354 http://dx.doi.org/10.1038/32354 http://dx.doi.org/10.1103/PhysRevA.86.063611 http://dx.doi.org/10.1103/PhysRevA.86.063611 http://dx.doi.org/10.1103/PhysRevA.86.063611 http://dx.doi.org/10.1103/PhysRevA.86.063611 http://dx.doi.org/10.1103/PhysRevA.87.013618 http://dx.doi.org/10.1103/PhysRevA.87.013618 http://dx.doi.org/10.1103/PhysRevA.87.013618 http://dx.doi.org/10.1088/0953-4075/47/1/015302 http://dx.doi.org/10.1088/0953-4075/47/1/015302 http://dx.doi.org/10.1088/0953-4075/47/1/015302 http://dx.doi.org/10.1088/0953-4075/47/1/015302 http://dx.doi.org/10.1103/PhysRevA.88.043603 http://dx.doi.org/10.1103/PhysRevA.88.043603 http://dx.doi.org/10.1103/PhysRevA.88.043603 http://dx.doi.org/10.1103/PhysRevA.88.043603 http://dx.doi.org/10.1103/PhysRevA.86.033606 http://dx.doi.org/10.1103/PhysRevA.86.033606 http://dx.doi.org/10.1103/PhysRevA.86.033606 http://dx.doi.org/10.1103/PhysRevA.86.033606 http://dx.doi.org/10.1103/PhysRevLett.102.230403 http://dx.doi.org/10.1103/PhysRevLett.102.230403 http://dx.doi.org/10.1103/PhysRevLett.102.230403 http://dx.doi.org/10.1103/PhysRevLett.102.230403 http://dx.doi.org/10.1016/j.cpc.2009.04.015 http://dx.doi.org/10.1016/j.cpc.2009.04.015 http://dx.doi.org/10.1016/j.cpc.2009.04.015 http://dx.doi.org/10.1016/j.cpc.2009.04.015 http://dx.doi.org/10.1016/j.cpc.2012.03.022 http://dx.doi.org/10.1016/j.cpc.2012.03.022 http://dx.doi.org/10.1016/j.cpc.2012.03.022 http://dx.doi.org/10.1016/j.cpc.2012.03.022 http://dx.doi.org/10.1103/PhysRevA.66.023613 http://dx.doi.org/10.1103/PhysRevA.66.023613 http://dx.doi.org/10.1103/PhysRevA.66.023613 http://dx.doi.org/10.1103/PhysRevA.66.023613 http://dx.doi.org/10.1103/PhysRevA.63.053607 http://dx.doi.org/10.1103/PhysRevA.63.053607 http://dx.doi.org/10.1103/PhysRevA.63.053607 http://dx.doi.org/10.1103/PhysRevA.63.053607 http://dx.doi.org/10.1103/PhysRevLett.107.073202 http://dx.doi.org/10.1103/PhysRevLett.107.073202 http://dx.doi.org/10.1103/PhysRevLett.107.073202 http://dx.doi.org/10.1103/PhysRevLett.107.073202 http://dx.doi.org/10.1016/j.physleta.2005.07.044 http://dx.doi.org/10.1016/j.physleta.2005.07.044 http://dx.doi.org/10.1016/j.physleta.2005.07.044 http://dx.doi.org/10.1016/j.physleta.2005.07.044