Demixing and symmetry breaking in binary dipolar Bose-Einstein condensate solitons S. K. Adhikari∗1 1 Instituto de F́ısica Teórica, UNESP - Universidade Estadual Paulista, 01.140-070 São Paulo, São Paulo, Brazil We demonstrate 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 repul- sive 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 the species occupies the central region, whereas the other species separates into two equal parts and stay pre- dominantly 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. PACS numbers: 03.75.Hh, 03.75.Mn, 03.75.Kk, 03.75.Lm I. INTRODUCTION Quasi-one-dimensional (quasi-1D) bright soliton and soliton train were created and investigated experimen- tally in Bose-Einstein condensate (BEC) of 7Li [1] and 85Rb atoms [2]. These quasi-1D solitons appear for at- tractive contact interaction in an axially-free BEC under radial harmonic trap [3]. However, due to collapse insta- bility, in three dimensions (3D), such solitons are fragile and can accommodate only a small number of atoms. Recent observation of BECs of 164Dy [4, 5], 168Er [6] and 52Cr [7, 8] atoms with large magnetic dipole mo- ments has opened new directions of research in BEC soli- tons. In dipolar BECs, in addition to the quasi-1D soli- tons [9, 10], one could also have quasi-two-dimensional (quasi-2D) solitons [11] and vortex solitons [12]. More- over, 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] config- urations. More interestingly, one can have dipolar BEC solitons for fully repulsive contact interaction [9]. Hence these solitons stabilized solely by long-range dipolar in- teraction are less vulnerable to collapse in 3D for a large number of atoms due to the repulsive contact interac- tion. This creates a new scenario for robust solitons of large number of atoms stabilized by short-range repul- sion and long-range inter- and intraspecies dipolar attrac- tion. One can study the interplay between the long-range anisotropic dipolar interactions and variable contact in- teractions 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 re- pulsive interspecies contact interaction, we demonstrate fully demixed (separated) bright dipolar binary solitons with a minimum of spatial overlap between the two ∗adhikari@ift.unesp.br; URL: http://www.ift.unesp.br/users/adhikari species stabilized by long-range dipolar interactions for repulsive intraspecies contact interactions. The demixed binary soliton could be spatially symmetric while one of the species occupies the central region (|z| < z0) and the other species stays symmetrically in both peripheral regions (|z| > z0). For the same set of parameters (num- ber 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 at- tractive 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 re- alistic values of inter- and intraspecies contact and dipo- lar 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 critical number of atoms and the interspecies scattering length. The profiles of the binary soliton are displayed in 3D iso- density 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 con- tact interaction is perturbed with a small time-dependent oscillating part. A steady propagation of the binary soli- ton 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 studied in laboratory and we suggest a way of achieving this goal. First, a highly cigar-shaped dipo- lar BEC with appropriate number of atoms has to be ar X iv :1 40 1. 31 84 v1 [ co nd -m at .q ua nt -g as ] 1 4 Ja n 20 14 http://www.ift.unesp.br/users/adhikari 2 created. Then the axial trap is to be removed slowly and linearly is a small interval of time while 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 a very slow turn-off of the axial trap will always lead to a stable soliton. However, a rea- sonably quick turn-off 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 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 viabil- ity of preparing these solitons experimentally by relaxing the axial trap on a quasi-1D binary mixture is demon- strated by real-time simulation of the mean-field model. A similar demonstration for a single-component dipolar quasi-1D soliton was also made. Finally, in Sec. IV a brief summary of our findings is presented. II. MEAN-FIELD MODEL FOR BINARY MIXTURE The extension of the mean-field Gross-Pitaevskii (GP) equation to a 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, num- ber of atoms, magnetic dipole moment, and scattering length for the two species i = 1, 2, given by mi, Ni, µi, ai, respectively. The intra- (Vi) and interspecies (V12) inter- actions 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 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 imi 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) ddm12, g (1) dd = 3N1a (1) dd , g (12) dd = 3N2a (12) dd m1/(2mR), g (21) dd = 3N1a (12) dd m1/(2mR). 3 In Eqs. (7) and (8), length is expressed in units of oscilla- tor length l0 = √ ~/(m1ω1), energy in units of oscillator energy ~ω1, probability density |φi|2 in units of l−30 , and time in units of t0 = 1/ω1. 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, add the dipolar length. III. NUMERICAL RESULTS We demonstrate stable demixed spatially-symmetric and asymmetric binary dipolar solitons for realistic val- ues 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 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, with a0 the Bohr radius. For the 164Dy-162Dy mixture, the 164Dy atoms are la- beled i = 1 and the 162Dy atoms are labeled i = 2. The magmetic moment of a 162Dy atom is µ2 = 10µB . Consequently the dipolar lengths are: 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: µB = 9.27401 × 10−24 Am2, a0 = 5.292 × 10−11 m, µ0 = 4π × 10−7 N/A2, ~ = 1.05457 × 10−34 m2kg/s, 1 amu = 1.66054× 10−27 kg. The dipolar interaction in 164Dy atoms is roughly double of 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 fre- quencies ω1 = ω2 = 2π×61.6 Hz, so that the length scale l0 ≡ √ ~/m1ω1 ≈ 1 µm, and time scale t0 ≡ ω−11 ≈ 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 ther- malization 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 coordi- nates independent of the underlying trap symmetry using 10 2 10 3 10 4 10 5 0 20 40 60 80 100 120 140 N c ri t a/a0 Collapse Stable Stable Collapse Stable U n b o u n d U n b o u n d a dd /a0 Dy (n) 10 2 10 3 10 4 10 5 0 20 40 60 80 100 120 140 N c ri t a/a0 Collapse Stable Stable Collapse Stable U n b o u n d U n b o u n d a dd /a0 Dy (n) Dy (v) 10 2 10 3 10 4 10 5 0 20 40 60 80 100 120 140 N c ri t a/a0 Collapse Stable Stable Collapse Stable U n b o u n d U n b o u n d a dd /a0 Dy (n) Dy (v) Er (n) 10 2 10 3 10 4 10 5 0 20 40 60 80 100 120 140 N c ri t a/a0 Collapse Stable Stable Collapse Stable U n b o u n d U n b o u n d a dd /a0 Dy (n) Dy (v) Er (n) Er (v) 10 2 10 3 10 4 10 5 0 20 40 60 80 100 120 140 N c ri t a/a0 Collapse Stable Stable Collapse Stable U n b o u n d U n b o u n d a dd /a0 10 2 10 3 10 4 10 5 0 20 40 60 80 100 120 140 N c ri t a/a0 Collapse Stable Stable Collapse Stable U n b o u n d U n b o u n d a dd /a0 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 sys- tem 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 l0 = 1 µm. 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 usual fashion [20]. We study demixing in binary solitons supported by long-range dipolar interactions for repulsive interspecies and intraspecies contact interactions. The repulsive con- tact interactions make the dipolar solitons less vulnerable to collapse for a larger number of atoms compared to the nondipolar solitons. 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 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 max- imum number of atoms, beyond which the system col- lapses. 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 devel- oped 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 of contact 4 |φ(z,t)|2 0 100 200 300 t/t0 -150 -100 -50 0 50 100 150 z/l0 0.01 0.02 0.03 0.04 0.05 FIG. 2: (Color online) Integrated 1D density |φ(z, t)|2 =∫ dx ∫ dy|φ(r, t)|2 of 164Dy atoms during real-time propaga- tion when the axial trap of angular frequency ω = 2π × 2 Hz on a quasi-1D dipolar BEC of 8000 164Dy atoms is re- moved linearly for 20 > t/t0 > 0 and the resultant axially- free soliton is propagated. Parameters used in simulation: add = 132.7a0, a = 120a0, ωρ = 2π × 61 Hz, l0 = 1 µm. repulsion as a → add, which is counterintuitive. The solitons are bound by long-range dipolar interaction and increase of 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) con- tact repulsion dominates over dipolar attraction and the soliton cannot be bound. In the collapse region, the op- posite 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 shall take a . add. In this fashion, a binary soliton of large number of atoms could be created and this would be of greater experimen- tal interest. In this study we take for the 164Dy atoms a1 = 120a0, and for the 168Er atoms 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 rel- atively 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 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 is gradually (linearly) 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 large interval of time. Long sustained propa- gation of the central peak establishes the soliton nature of the axially-free dipolar BEC. The result of this sim- ulation is presented in Fig. 2, where we show the inte- grated 1D density |φ(z, t)|2 = ∫ dx ∫ dy|φ(r, t)|2 during real-time propagation. Upon the 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 contact attraction only cannot be easily prepared in this fashion. B. Binary 164Dy-168Er soliton Scattering lengths play important roles in the prepa- ration of binary solitons and can be experimentally con- trolled independently 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 numerical simulation we realized that a reasonably large interspecies scattering length a12 (& 70a0) and the num- ber of atoms of the species below a critical value is 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 val- ues of the number of atoms, the system collapses due to a strong net dipolar interaction even for repulsive in- traspecies interaction. The demixing takes place due to interspecies contact repulsion. For smaller values of in- terspecies scattering length, the repulsion is not sufficient for demixing and a mixed or overlapping configuration of binary dipolar BEC soliton is realized [16]. A spatially-symmetric binary soliton is obtained by solving Eqs. (7) and (8) with imaginary-time propa- gation with the following spatially-symmetric Gaussian input for both the 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 soli- ton are spatially-symmetric around z = 0. To obtain a symmetry-broken binary soliton, a symmetry-broken ini- tial input is needed. Now we display the typical profile of a spatially- symmetric demixed binary dipolar BEC soliton in Fig. 5 FIG. 3: (Color online) 3D 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, 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. 3 for 2000 164Dy atoms and 5000 168Er atoms for the interspecies scattering length a12 = a(Dy-Er) = 90a0 and for large intraspecies scattering lengths: a1 = a(Dy) = 120a0, a2 = a(Er) = 60a0. The species containing the larger number of atoms (168Er) stay at the center and the species containing the smaller number of atoms (164Dy) breaks into two equal parts, leave the central region oc- cupied by 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 pa- per 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 ap- pearance of the demixed binary dipolar soliton in the 164Dy-168Er mixture for different number of atoms and variable interspecies scattering length a12. In Figs. 4 (a) and (b) we show the critical number of 164Dy atoms N1 = N(Dy) versus a12 for the numbers of 168Er atoms N2 = N(Er) = 1000, 5000, respectively. The binary sys- tem 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 in- terspecies contact repulsion. A Gaussian variational ap- proximation 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 Figs. 4. As the variational approximation overbinds the binary soli- ton, 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 the contact repulsion between the two species of atoms increases and the system gets demixed. There are two ways this demix- ing takes place. The demixed binary soliton can either 10 2 10 3 10 4 10 5 60 70 80 90 100 110 120 N (D y ) a(Dy-Er)/a0 N(Er) = 1000 collapse 10 2 10 3 10 4 10 5 60 70 80 90 100 110 120 N (D y ) a(Dy-Er)/a0 N(Er) = 1000 collapse 10 2 10 3 10 4 10 5 60 70 80 90 100 110 120 N (D y ) a(Dy-Er)/a0 N(Er) = 1000 collapse 10 2 10 3 10 4 10 5 60 70 80 90 100 110 120 N (D y ) a(Dy-Er)/a0 N(Er) = 1000 collapse demixed demixed mixed (a) variational 10 2 10 3 10 4 10 5 60 70 80 90 100 110 120 N (D y ) a(Dy-Er)/a0 N(Er) = 5000 collapse 10 2 10 3 10 4 10 5 60 70 80 90 100 110 120 N (D y ) a(Dy-Er)/a0 N(Er) = 5000 collapse 10 2 10 3 10 4 10 5 60 70 80 90 100 110 120 N (D y ) a(Dy-Er)/a0 N(Er) = 5000 collapse 10 2 10 3 10 4 10 5 60 70 80 90 100 110 120 N (D y ) a(Dy-Er)/a0 N(Er) = 5000 collapse demixed demixed mixed(b) variational FIG. 4: (Color online) Stability phase diagram for a binary 164Dy-168Er soliton for N(Er) = (a) 1000 and (b) 5000 atoms. The variational mixed-collapse border is shown by solid cir- cles. The oscillator length l0 = 1 µm, a(Dy)= 120a0, a(Er)= 60a0. Stable demixed binary soliton is observed in the darker (pink and green) region to the right while stable mixed bi- nary soliton appears in the lighter (yellow) region to the left. For a spatially-symmetric demixed binary soliton, in the up- per 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). be spatially-symmetric around z = 0, or it can be spa- tially asymmetric around z = 0. The two configurations are possible for identical parameters of the binary sys- tem, e.g., number of atoms of the two components, and all dipolar and contact interaction strengths. However, there is a region aroundN1 ≈ 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 compo- nents (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 configura- tion 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 configura- tion. In Figs. 5 (a) and (b) we show the critical number of 168Er atoms N2 = N(Er) versus a12 for the numbers of 164Dy atoms N1 = N(Dy) = 1000, 5000, respectively. A similar scenario emerges in Figs. 5 as in Figs. 4. 6 10 2 10 3 10 4 10 5 60 70 80 90 100 110 120 N (D y ) a(Dy-Er)/a0 N(Dy) = 1000 collapse 10 2 10 3 10 4 10 5 60 70 80 90 100 110 120 N (D y ) a(Dy-Er)/a0 N(Dy) = 1000 collapse 10 2 10 3 10 4 10 5 60 70 80 90 100 110 120 N (D y ) a(Dy-Er)/a0 N(Dy) = 1000 collapse 10 2 10 3 10 4 10 5 60 70 80 90 100 110 120 N (D y ) a(Dy-Er)/a0 N(Dy) = 1000 collapse demixed demixed mixed (a) variational 10 2 10 3 10 4 10 5 60 70 80 90 100 110 120 N (D y ) a(Dy-Er)/a0 N(Dy) = 5000 collapse 10 2 10 3 10 4 10 5 60 70 80 90 100 110 120 N (D y ) a(Dy-Er)/a0 N(Dy) = 5000 collapse 10 2 10 3 10 4 10 5 60 70 80 90 100 110 120 N (D y ) a(Dy-Er)/a0 N(Dy) = 5000 collapse 10 2 10 3 10 4 10 5 60 70 80 90 100 110 120 N (D y ) a(Dy-Er)/a0 N(Dy) = 5000 collapse 10 2 10 3 10 4 10 5 60 70 80 90 100 110 120 N (D y ) a(Dy-Er)/a0 N(Dy) = 5000 collapse 10 2 10 3 10 4 10 5 60 70 80 90 100 110 120 N (D y ) a(Dy-Er)/a0 N(Dy) = 5000 collapse demixed mixed demixed (b) variational FIG. 5: (Color online) Same as in Fig. 4 for (a) N(Dy) = 1000 and (b) 5000 atoms. In Fig. 3 we presented a spatially-symmetric demixed binary dipolar soliton with 2000 164Dy and 5000 168Er atoms. In this case, 168Er atoms lie in the central re- gion 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 iso- density contour of the binary soliton the roles of 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-symmetric configuration the 164Dy atoms stay at the center and 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 small variation of the parameters of the initial state. To achieve this, we con- sider the final state of the spatially-symmetric binary soli- ton 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 scatter- ing length a12 = 90a0 + 10a0 cos(t/t0), while 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 no- FIG. 6: (Color online) 3D 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, 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. FIG. 7: (Color online) Contour plot of integrated 1D densi- ties |φi(z, t)|2 = ∫ dx ∫ dy|φ(r, t)|2 of 164Dy and 168Er atoms when the binary soliton of Fig. 3 is propagated 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. ticeable 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 = ∫ ∫ dxdy|φ(x, y, z, t)|2 during real- time propagation. The visual profile of the solitons in Fig. 7 maintain a constant width and do 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 dipo- lar 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 7 0 0.01 0.02 0.03 0.04 0.05 -60 -40 -20 0 20 40 60 1 D d e n s it y | φ i( z )| 2 z/l0 Dy (i) 0 0.01 0.02 0.03 0.04 0.05 -60 -40 -20 0 20 40 60 1 D d e n s it y | φ i( z )| 2 z/l0 Dy (i) Er (i) 0 0.01 0.02 0.03 0.04 0.05 -60 -40 -20 0 20 40 60 1 D d e n s it y | φ i( z )| 2 z/l0 Dy (f) 0 0.01 0.02 0.03 0.04 0.05 -60 -40 -20 0 20 40 60 1 D d e n s it y | φ i( z )| 2 z/l0 Dy (f) Er (f) FIG. 8: (Color online) Initial (i) at t = 0 and final (f) 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. FIG. 9: (Color online) Same as in Fig. 3 for the symmetry- broken binary soliton of 2000 164Dy and 5000 168Er atoms. by imaginary-time propagation one should consider spatially-asymmetric initial states. In particular we con- sider the following spatial-symmetry broken (asymmetric around z = 0) initial states with a spatial displacement 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 con- verges to the lowest-energy state maintaining the sym- metry of the initial state. For example, in the 1D har- monic oscillator problem the imaginary-time propaga- tion with a spatially-symmetric initial state will con- verge to the ground state of the problem, whereas with a spatially-antisymmetric initial state it will converge to FIG. 10: (Color online) Same as in Fig. 3 for the symmetry- broken binary soliton of 5000 164Dy and 1000 168Er atoms. 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 bi- nary solitons. We tried some other forms of symmetry- broken initial state and found that the solution converges always to the same final state. After some experimenta- tion, we could locate for the same set of parameters of the problem 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 ini- tial 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 resul- tant 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 sec- ond 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 pa- rameters 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. In case of demixed solitons, it is interesting to ask which of the two configurations − symmetric or asym- metric − has the lower energy and hence will be exper- imentally favorable. To this end we calculated numeri- cally the respective energies in several cases and found that consistently the spatially asymmetric configuration has slightly lower energy. However, this will not be of phenomenological interest as the energy difference is very small, being of 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 8 FIG. 11: (Color online) Same as in Fig. 7 for the symmetry- broken binary soliton of Fig. 10. and use these as the initial states in 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 conden- sates 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, e 10 are ideal examples of soliton molecules of dif- ferent species where the participating solitons maintain their identities with a minimum of overlap. The binding between the two solitons come from the long-range dipo- lar interaction. In nondipolar systems such binding can only come from short-range interspecies attraction and a composite soliton molecule can appear possibly only with complete overlap between participating solitons [16, 22]. The demixed binary dipolar soliton or the binary soli- ton molecule can be prepared and observed experimen- tally. We illustrate the possibility of experimental real- ization 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 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 pro- gram. 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 contin- 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 angu- lar 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 is removed linearly for 20 > t/t0 > 0 and the resultant soliton is propagated. Parameters used in simulation: a(Dy) = 120a0, a(Er) = 60a0, a(Dy-Er) = 90a0, ωρ = 2π × 61 Hz, l0 = 1 µm. ued 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 propaga- tion. A similar simulation for the spatially-asymmetric binary soliton of Fig. 9 is shown in Fig. 12 (b). Af- ter 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 (b). 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 laboratory [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) 9 10 2 10 3 10 4 10 5 60 70 80 90 100 110 120 130 140 N (1 6 4 D y ) a( 162 Dy- 164 Dy)/a0 N( 162 Dy) = 1000 collapse 10 2 10 3 10 4 10 5 60 70 80 90 100 110 120 130 140 N (1 6 4 D y ) a( 162 Dy- 164 Dy)/a0 N( 162 Dy) = 1000 collapse 10 2 10 3 10 4 10 5 60 70 80 90 100 110 120 130 140 N (1 6 4 D y ) a( 162 Dy- 164 Dy)/a0 N( 162 Dy) = 1000 collapse 10 2 10 3 10 4 10 5 60 70 80 90 100 110 120 130 140 N (1 6 4 D y ) a( 162 Dy- 164 Dy)/a0 N( 162 Dy) = 1000 collapse demixed demixed mixed variational FIG. 13: (Color online) 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. Parameters used: l0 = 1 µm, a(164Dy)= a(164Dy)= 120a0. =120a0, viz. Fig. 1. The interspecies scattering length is considered as a variable. First we consider the sta- bility phase plot for this binary soliton. In Fig. 13 we show the number of 164Dy atoms in the binary soliton for 1000 162Dy atoms. Again there appears mixed and demixed binary solitons. As the mass and dipolar lengths are almost the same for the two isotopes the binary plot is quasi symmetric. The plot for 1000 164Dy atoms in the binary soliton will lead to a practiaclly identical plot as that in Fig. 13 with the role of the two isotopes in- terchanged. 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 appears for a larger value of interspecies scattering length to compen- sate 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 also one has mixed, spatially-symmetric demixed, and spatial symmetry-broken demixed binary solitons. Without repeating a detailed discussion of dif- ferent types of solitons, in Fig. 14 we show the iso- density 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. FIG. 14: (Color online) 3D 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, 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 inter- actions 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 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 bi- nary solitons. In the spatially-symmetric case the species with smaller number of atoms breaks up into two equal parts and the parts leave the central region occupied by the species with larger number of atoms and consolidate on both sides of the BEC component with 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 in- teraction 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 10 amplitude. A stable binary soliton should then execute small oscillations. The converged solution of the binary soliton in imaginary-time propagation is used as the ini- tial state of real-time propagation. The solitons are found to exhibit small oscillation over a long time thus estab- lishing their stability. The solitons considered in this paper are stabilized by long-range dipolar attraction and short-range contact re- pulsion. The dipolar interaction is attractive in the axial direction and the dipolar repulsion in the transverse di- rection is compensated by a harmonic trap. Hence unlike normal BEC solitons stabilized by short-range contact attraction, the present dipolar BECs will be more im- mune to collapse due to short-range repulsion and can easily accommodate 10000 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 these 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 col- lapse instability. In a stable nondipolar soliton, the max- imum 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 experi- mental 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, while the demixed binary BEC turns into a demixed bi- nary soliton. With the present experimental techniques, such binary dipolar BECs can be observed and the con- clusions of this study verified. Acknowledgments We thank FAPESP and CNPq (Brazil) for partial sup- port. [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, Phys. Rev. Lett. 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, Seo Ho Youn, and B. L. Lev, Phys. Rev. Lett. 107, 190401 (2011). [6] K. Aikawa et al., Phys. Rev. Lett. 108, 210401 (2012). [7] T. Lahaye et al., Nature (London) 448, 672 (2007); A. Griesmaier et al., Phys. Rev. Lett. 97, 250402 (2006). [8] J. Stuhler et al., Phys. Rev. Lett. 95, 150406 (2005); K. Goral, K. Rzazewski, and T. Pfau, Phys. Rev. A 61, 051601 (2000); T. Koch et al., Nature Phys. 4, 218 (2008); T. Lahaye et al., Rep. Prog. Phys. 72, 126401 (2009). [9] Luis 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, Phys. Rev. Lett. 95, 200404 (2005); I. I. Tikhonenkov, B. A. Malomed, and A. Vardi, Phys. Rev. Lett. 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, Phys. Rev. A 86, 053611 (2012); R. Eichler, J. Main, and G. Wunner, Phys. Rev. A 83, 053604 (2011); Ai-Xia Zhang and Ju-Kui Xue, Phys. Rev. A 82, 013606 (2010); V. M. Lashkin, Phys. Rev. A 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); Phys. Rev. A 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. Timmer- mans, 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, Com- put. Phys. Commun. 183, 2021 (2012). [20] K. Goral and L. Santos, Phys. Rev. A. 66, 023613 (2002); S. Yi and L. You, Phys. Rev. A 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, (2013). I Introduction II Mean-field model for binary mixture III Numerical Results A Single component 164Dy and 168Er solitons B Binary 164Dy-168Er soliton C Binary 164Dy-162Dy soliton IV Summary and Discussion Acknowledgments References