PHYSICAL REVIEW A 87, 013618 (2013) Dipolar droplet bound in a trapped Bose-Einstein condensate Luis E. Young-S.* and S. K. Adhikari† Instituto de Fı́sica Teórica, UNESP–Universidade Estadual Paulista, 01.140-070 São Paulo, São Paulo, Brazil (Received 9 November 2012; published 17 January 2013) We study the statics and dynamics of a dipolar Bose-Einstein condensate (BEC) droplet bound by interspecies contact interaction in a trapped nondipolar BEC. Our findings are demonstrated in terms of stability plots of a dipolar 164Dy droplet bound in a trapped nondipolar 87Rb BEC with a variable number of 164Dy atoms and interspecies scattering length. A trapped nondipolar BEC of a fixed number of atoms can bind only a dipolar droplet containing fewer atoms than a critical number for the interspecies scattering length between two critical values. The shape and size (statics) as well as the small breathing oscillation (dynamics) of the dipolar BEC droplet are studied using numerical and variational solutions of a mean-field model. We also suggest an experimental procedure for achieving such a 164Dy droplet by relaxing the trap on the 164Dy BEC in a trapped binary 87Rb-164Dy mixture. DOI: 10.1103/PhysRevA.87.013618 PACS number(s): 03.75.Hh, 03.75.Mn, 03.75.Kk I. INTRODUCTION The experimental observation of dipolar Bose-Einstein condensates (BECs) of 52Cr [1–5], 164Dy [6,7], and 168Er [8] atoms with large magnetic dipole moments has opened new directions of research in cold atoms in the quest of novel and in- teresting features related to the anisotropic long-range dipolar interaction. Polar molecules, with much larger (electric) dipole moment, are also being considered [9] for BEC experiments. The atomic interaction in a dilute BEC of alkali metal and other types of atoms (with negligible dipole moment) is represented by an S-wave contact (δ-function) potential. However, the nonlocal anisotropic long-range dipolar interaction acts in all partial waves and is attractive in certain directions and repulsive in others. An untrapped three-dimensional (3D) BEC with attractive interaction does not exist in nature due to collapse instability [10]. However, the collapse of an untrapped BEC can be avoided due to the interspecies contact interaction in a binary mixture with a trapped BEC. The shape of such a stable droplet bound in a trapped BEC is controlled by the interspecies interactions, whereas that of a trapped BEC is determined by the underlying trap. Here we consider such a dipolar droplet bound in a trapped nondipolar BEC. The effect of the underlying atomic interaction, especially that of the anisotropic dipolar interaction, will easily manifest itself in such a dipolar droplet. These droplets will be termed quasifree as they can easily move around inside the larger trapped BEC responsible for their binding. By taking the trapped BEC to be nondipolar, one can easily study the effect of intraspecies dipolar interaction on the quasifree droplet in the absence of any interspecies dipolar interaction. Dipolar BECs are immediately distinguishable from those with purely contact interactions by their strong shape and stability dependence [11] on the trapping geometry. Here, we are proposing a different way to trap a dipolar BEC in the form of a quasifree dipolar droplet, using an attractive interspecies mean-field potential, which introduces a unique *lyoung@ift.unesp.br †adhikari@ift.unesp.br; URL: http://www.ift.unesp.br/users/adhikari trapping geometry that results in further interesting shape and stability characteristics of dipolar BECs. We study the statics and dynamics of a quasifree dipolar droplet using numerical solution and variational approximation of a mean- field model [12,13]. For this purpose we consider a binary mixture of nondipolar 87Rb and dipolar 164Dy where the trapped nondipolar 87Rb BEC could be in a cigar or disk shape. Among the available dipolar BECs [4,7,8], 164Dy atoms have the largest (magnetic) dipolar interaction strength [1,6]. The existence of stable 3D 164Dy droplets is illustrated by stability plots involving the number of atoms of the two species and the interspecies scattering length. For a fixed number N (Rb) of 87Rb atoms the dipolar droplet could be bound below a critical number N (Dy) of 164Dy atoms between two limiting values of interspecies attraction. If the interspecies attraction is too small, the 164Dy atoms in the droplet cannot be bound and, for a very large interspecies attraction, the droplet is destroyed by collapse instability. Usually, the dipolar droplet has the same shape as the trap acting on the nondipolar BEC. However, for a large number of 164Dy atoms, as the interspecies attraction approaches the collapse instability, the dipolar droplet always is of cigar shape even if the external trap is disk shaped, and eventually the dipolar droplet collapses on the axial z axis from the cylindrical configuration. For a small number of atoms, the dipolar droplet collapses to the center, maintaining its shape, rather than on the axial z axis. The variational approximation to the sizes and chemical potentials of the stationary droplets is compared with the numerical solution of the mean-field model. The numerical study of breathing oscillation of the stable dipolar droplet is found to be in reasonable agreement with a time-dependent variational model calculation. We also demonstrate a viable experimental way of creating a 164Dy droplet bound in a trapped 87Rb BEC, e.g., by slowly removing the trap on a trapped binary 164Dy-87Rb mixture, while the trapped 164Dy BEC evolves into a quasifree droplet. In Sec. II the mean-field model for a dipolar BEC droplet bound in a trapped nondipolar BEC is developed. A time- dependent, analytic, Euler-Lagrange Gaussian variational ap- proximation of the model is also presented. The results of numerical calculation are shown in Sec. III. Finally, in Sec. IV we present a brief summary of our findings. 013618-11050-2947/2013/87(1)/013618(7) ©2013 American Physical Society http://dx.doi.org/10.1103/PhysRevA.87.013618 LUIS E. YOUNG-S. AND S. K. ADHIKARI PHYSICAL REVIEW A 87, 013618 (2013) II. MEAN-FIELD MODEL FOR A QUASIFREE DIPOLAR DROPLET We consider a binary BEC, where one of the species is dipolar and the other nondipolar, interacting via inter- and intraspecies interactions with the mass, number, magnetic dipole moment, and scattering length for the two species i = 1,2, denoted by mi, Ni, μi, and ai, respectively. The first species (87Rb) is taken to be nondipolar (μ1 = 0) and trapped while the second species (164Dy) is dipolar (μ2 �= 0) and polarized along the axial z direction. The angular frequencies for the axially symmetric trap on 87Rb along the x, y, and z directions are taken as ωx = ωy = ω1 and ωz = λ1ω1. The inter- (V12) and intra- (Vi) species interactions for two atoms at positions r and r′ are taken as V12(R) = 2πh̄2a12 mR δ(R), V1(R) = 4πh̄2a1 m1 δ(R), (1) V2(R) = μ0μ 2 2 4π 1 − 3 cos2 θ |R|3 + 4πh̄2a2 m2 δ(R), (2) where R = r − r′, μ0 is the permeability of free space, θ is the angle made by the vector R with the polarization z direction, and mR = m1m2/(m1 + m2) is the reduced mass of the two species of atoms. With these interactions, the coupled Gross-Pitaevskii (GP) equations for the binary dipolar BEC can be written as [13] ih̄ ∂φ1(r,t) ∂t = [ − h̄2 2m1 ∇2 + 1 2 m1ω 2 1 ( ρ2 + λ2 1z 2 ) + 4πh̄2 m1 a1N1|φ1(r,t)|2 + 2πh̄2 mR a12N2|φ2(r,t)|2 ] φ1(r,t), (3) ih̄ ∂φ2(r,t) ∂t = [ − h̄2 2m2 ∇2 + 4πh̄2 m2 a2N2|φ2(r,t)|2 + 2πh̄2 mR a12N1|φ1(r,t)|2 + N2 μ0 μ2 2 4π × ∫ Vdd (R)|φ2(r′,t)|2dr′ ] φ2(r,t), (4) Vdd (R) = 1 − 3 cos2 θ R3 , ρ2 = x2 + y2, i = √−1. (5) To compare the dipolar and contact interactions, the intraspecies dipolar interaction is expressed in terms of the length scale add , defined by add = μ0μ 2 2m2/(12πh̄2). We express the strength of the dipolar interaction in Eq. (4) by this length scale and transform Eqs. (3) and (4) into the following dimensionless form [13]: i ∂φ1(r,t) ∂t = [ − ∇2 2 + 1 2 ( ρ2 + λ2 1z 2 ) + g1|φ1|2 + g12|φ2|2 ] φ1(r,t), (6) i ∂φ2(r,t) ∂t = [ − m12 ∇2 2 + g2|φ2|2 + g21|φ1|2 + gdd ∫ Vdd (R)|φ2(r′,t)|2dr′ ] φ2(r,t), (7) where m12 = m1/m2,g1 = 4πa1N1,g2 = 4πa2N2m12, g12 = 2πm1a12N2/mR,g21 = 2πm1a12N1/mR, and gdd = 3N2addm12. In Eqs. (6) and (7), length is expressed in units of the oscillator length of the first species l0 = √ h̄/m1ω1, energy in units of the oscillator energy h̄ω1, density |φi |2 in units of l−3 0 , and time in units of t0 = ω−1 1 . A convenient analytic variational approximation to Eqs. (6) and (7) can be obtained with the following ansatz for the wave functions [14–16]: φi(r,t) = π−3/4 wρi √ wzi exp [ − ρ2 2w2 ρi − z2 2w2 zi + iαiρ 2 + iβiz 2 ] , (8) where r = { �ρ,z}, �ρ = {x,y}, wρi and wzi are the widths, and αi and βi are additional variational parameters. The effective Lagrangian for the binary system is L = ∫ dr 1 2 [ ∑ i {iNi(φiφ̇ ∗ i − φ∗ i φ̇i) + Nigi |φi(r)|4} +N1|∇φ1(r)|2 + m12N2|∇φ2(r)|2 + N1 ( ρ2 + λ2 1z 2 ) × |φ1(r)|2 ] + N1g12 ∫ dr|φ1(r)|2|φ2(r)|2 + N2 2 gdd ∫∫ Vdd (R)|φ2(r′)|2|φ2(r)|2dr′dr, (9) = 2∑ i=1 Ni 2 ( 2w2 ρi α̇i + w2 zi β̇i ) + N1 2 [ 1 w2 ρ1 + 1 2w2 z1 + 4w2 ρ1α 2 1 + 2w2 z1β 2 1 ] + N2m12 2 [ 1 w2 ρ2 + 1 2w2 z2 + 4w2 ρ2α 2 2 + 2w2 z2β 2 2 ] + N1 2 [ w2 ρ1 + λ2 1 w2 z1 2 ] + N2 1 a1√ 2πw2 ρ1wz1 + N2 2 m12√ 2πw2 ρ2wz2 [a2 − addf (κ)] + CN1N2 2AB1/2 , κ = wρ2 wz2 , (10) f (κ) = 1 + 2κ2 − 3κ2d(κ) 1 − κ2 , d(κ) = arctan( √ κ2 − 1)√ κ2 − 1 , where A = w2 ρ1 + w2 ρ2, B = w2 z1 + w2 z2, and C = 4a12m1/ ( √ πmR). In these equations the overdot denotes the time derivative. The Euler-Lagrange variational equations for the widths for the effective Lagrangian (10), obtained in the usual fashion [16], are ẅρ1 + wρ1 = 1 w3 ρ1 + 1√ 2π 2N1a1 w3 ρ1wz1 + CN2wρ1 A2B1/2 , (11) ẅz1 + λ2 1wz1 = 1 w3 z1 + 1√ 2π 2N1a1 w2 ρ1w 2 z1 + CN2wz1 AB3/2 , (12) ẅρ2 = m12 w3 ρ2 + N2m12[2a2 − addg(κ)]√ 2πw3 ρ2wz2 + CN1wρ2 A2B1/2 , (13) ẅz2 = m12 w3 z2 + 2N2m12[a2 − addh(κ)]√ 2πw2 ρ2w 2 z2 + CN1wz2 AB3/2 , (14) 013618-2 DIPOLAR DROPLET BOUND IN A TRAPPED BOSE- . . . PHYSICAL REVIEW A 87, 013618 (2013) g(κ) = 2 − 7κ2 − 4κ4 + 9κ4d(κ) (1 − κ2)2 , (15) h(κ) = 1 + 10κ2 − 2κ4 − 9κ2d(κ) (1 − κ2)2 . (16) The solution of the time-dependent Eqs. (11)–(14) gives the dynamics of the variational approximation. If μi is the chemical potential with which the stationary wave function φi(r,t) propagates in time, e.g., φi(r,t) ∼ exp(−iμit)φi(r), then the variational estimate for μi is μ1 = 1 2 [ 1 w2 ρ1 + 1 2w2 z1 ] + 1 2 [ w2 ρ1 + λ2 1 w2 z1 2 ] + 2N1a1√ 2πw2 ρ1wz1 + CN2 2AB1/2 , (17) μ2 = m12 2 [ 1 w2 ρ2 + 1 2w2 z2 ] + 2N2m12√ 2πw2 ρ2wz2 [a2 − addf (κ)] + CN1 2AB1/2 . (18) The total energy of the system is given by E = N1 2 [ 1 w2 ρ1 + 1 2w2 z1 ] + N1 2 [ w2 ρ1 + λ2 1 w2 z1 2 ] + N2m12 2 [ 1 w2 ρ2 + 1 2w2 z2 ] + N2 1 a1√ 2πw2 ρ1wz1 + N2 2 m12√ 2πw2 ρ2wz2 [a2 − addf (κ)] + CN1N2 2AB1/2 . (19) The widths of the stationary state can be obtained from the solution of Eqs. (11)–(14) by setting the time derivatives of the widths equal to zero. This procedure is equivalent to a minimization of the energy (19), provided the stationary state is stable and corresponds to an energy minimum. III. NUMERICAL RESULTS We solve Eqs. (6) and (7) by a split-step Crank-Nicolson discretization scheme using a space step of 0.1 and the time step 0.001 [15,17]. The contribution of the dipolar interaction is calculated in momentum space by Fourier transformation [15]. For both species of atoms 87Rb and 164Dy we take the intraspecies scattering length as ai = 110a0, and the strength of dipolar interaction as add = 131a0 [7,10]. The yet unknown interspecies scattering length a12 is taken as a variable. The variation of a12 can be achieved experimentally by the Feshbach resonance technique [18]. We consider the trap frequency ω1 = 2π × 115 Hz, so that the length scale l0 ≡ √ h̄/m1ω1 = 1 μm and time scale t0 ≡ ω−1 1 = 1.38 ms. We find that a quasifree 164Dy droplet is achievable for a moderately attractive interspecies attraction (negative a12) and for appropriate values of the number of atoms of the two species N1 and N2. We illustrate the domain of existence of a stable 164Dy droplet in terms of stability plots in Figs. 1(a), 1(b), and 1(c) for N (Rb) = 2000, 10 000, and 50 000, and for λ = 4 (disk-shaped trap) and 0.25 (cigar-shaped trap). We consider cigar- and disk-shaped 164Dy droplets in this study, where the effect of dipolar interaction is expected to be 0 40 80 120 160 102 103 104 |a 12 | ( un its o f a 0) N(Dy) unstable stable N(Rb) = 2000 (a) add = 131a0 λ = 4 0 40 80 120 160 102 103 104 |a 12 | ( un its o f a 0) N(Dy) unstable stable N(Rb) = 2000 (a) add = 131a0 λ = 4 λ = 0.25 0 40 80 120 160 102 103 104 |a 12 | ( un its o f a 0) N(Dy) unstable stable N(Rb) = 2000 (a) add = 131a0 λ = 4 λ = 0.25 0 40 80 120 160 102 103 104 |a 12 | ( un its o f a 0) N(Dy) unstable stable N(Rb) = 2000 (a) add = 131a0 λ = 4 λ = 0.25 0 40 80 120 160 102 103 104 |a 12 | ( un its o f a 0) N(Dy) unstable stable N(Rb) = 10000 (b) add = 131a0 λ = 4 0 40 80 120 160 102 103 104 |a 12 | ( un its o f a 0) N(Dy) unstable stable N(Rb) = 10000 (b) add = 131a0 λ = 4 λ = 0.25 0 40 80 120 160 102 103 104 |a 12 | ( un its o f a 0) N(Dy) unstable stable N(Rb) = 10000 (b) add = 131a0 λ = 4 λ = 0.25 0 40 80 120 160 102 103 104 |a 12 | ( un its o f a 0) N(Dy) unstable stable N(Rb) = 10000 (b) add = 131a0 λ = 4 λ = 0.25 0 40 80 120 160 102 103 104 |a 12 | ( un its o f a 0) N(Dy) unstable stable N(Rb) = 50000 (c) add = 131a0 λ = 4 0 40 80 120 160 102 103 104 |a 12 | ( un its o f a 0) N(Dy) unstable stable N(Rb) = 50000 (c) add = 131a0 λ = 4 λ = 0.25 0 40 80 120 160 102 103 104 |a 12 | ( un its o f a 0) N(Dy) unstable stable N(Rb) = 50000 (c) add = 131a0 λ = 4 λ = 0.25 0 40 80 120 160 102 103 104 |a 12 | ( un its o f a 0) N(Dy) unstable stable N(Rb) = 50000 (c) add = 131a0 λ = 4 λ = 0.25 FIG. 1. (Color online) |a12|-N (Dy) stability plot showing the domain of stable 164Dy droplets in a binary 87Rb-164Dy mixture for (a) N (Rb) = 2000, (b) N (Rb) = 10 000, and (c) N (Rb) = 50 000, and for λ = 0.25 and 4. The intraspecies scattering lengths are taken as ai = 110a0. more prominent, and not spherically symmetric ones where the effect of dipolar interaction is expected to be a minimum [1,5]. In all plots of Fig. 1, for a fixed N (Rb), the 164Dy droplet can be bound for a maximum of N (Dy). This maximum of N (Dy) increases with N (Rb), all other parameters remaining fixed. For a fixed N (Rb), the disk-shaped trap can accommodate a larger maximum number of 164Dy atoms than the cigar-shaped trap. For N (Dy) smaller than this maximum number, the 164Dy droplet can be bound for |a12| between two limiting values. For |a12| above the upper limit, there is too much interspecies attraction on the the 164Dy droplet, leading to its collapse. For |a12| below the lower limit, there is not enough attraction and the 164Dy droplet cannot be bound. The lower limit of |a12| is small and tends to zero as N (Rb) tends to infinity. For fixed values of N (Rb) and N (Dy) and for small |a12| near the lower limit, the disk-shaped configuration is favored, and it can bind the 164Dy droplet for smaller values of |a12|. However, for a fixed N (Rb) and for a large |a12| near the upper limit, the disk-shaped configuration is favored only for a large N (Dy) allowing a 164Dy droplet for larger |a12|, whereas the cigar-shaped configuration is favored for a small N (Dy) allowing a 164Dy droplet for larger |a12|. In Fig. 2 we illustrate the numerical and variational results for chemical potentials and root-mean-square (rms) sizes 〈x〉 013618-3 LUIS E. YOUNG-S. AND S. K. ADHIKARI PHYSICAL REVIEW A 87, 013618 (2013) -2 -1 0 1 2 3 4 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 4 N(Rb) = 10000 (a) μ1/4(v) -2 -1 0 1 2 3 4 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 4 N(Rb) = 10000 (a) μ1/4(v) (v) -2 -1 0 1 2 3 4 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 4 N(Rb) = 10000 (a) μ1/4(v) (v) (v) -2 -1 0 1 2 3 4 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 4 N(Rb) = 10000 (a) μ1/4(v) (v) (v) μ2(v) -2 -1 0 1 2 3 4 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 4 N(Rb) = 10000 (a) μ1/4(v) (v) (v) μ2(v) (v) -2 -1 0 1 2 3 4 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 4 N(Rb) = 10000 (a) μ1/4(v) (v) (v) μ2(v) (v) (v) -2 -1 0 1 2 3 4 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 4 N(Rb) = 10000 (a) μ1/4(n) -2 -1 0 1 2 3 4 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 4 N(Rb) = 10000 (a) μ1/4(n) (n) -2 -1 0 1 2 3 4 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 4 N(Rb) = 10000 (a) μ1/4(n) (n) (n) -2 -1 0 1 2 3 4 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 4 N(Rb) = 10000 (a) μ1/4(n) (n) (n) μ2(n) -2 -1 0 1 2 3 4 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 4 N(Rb) = 10000 (a) μ1/4(n) (n) (n) μ2(n) (n) -2 -1 0 1 2 3 4 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 4 N(Rb) = 10000 (a) μ1/4(n) (n) (n) μ2(n) (n) (n) -1 0 1 2 3 4 5 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 0.25 N(Rb) = 10000 (b) μ1(v) -1 0 1 2 3 4 5 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 0.25 N(Rb) = 10000 (b) μ1(v) (v) -1 0 1 2 3 4 5 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 0.25 N(Rb) = 10000 (b) μ1(v) (v) (v) -1 0 1 2 3 4 5 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 0.25 N(Rb) = 10000 (b) μ1(v) (v) (v) μ2(v) -1 0 1 2 3 4 5 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 0.25 N(Rb) = 10000 (b) μ1(v) (v) (v) μ2(v) (v) -1 0 1 2 3 4 5 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 0.25 N(Rb) = 10000 (b) μ1(v) (v) (v) μ2(v) (v) (v) -1 0 1 2 3 4 5 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 0.25 N(Rb) = 10000 (b) μ1(n) -1 0 1 2 3 4 5 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 0.25 N(Rb) = 10000 (b) μ1(n) (n) -1 0 1 2 3 4 5 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 0.25 N(Rb) = 10000 (b) μ1(n) (n) (n) -1 0 1 2 3 4 5 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 0.25 N(Rb) = 10000 (b) μ1(n) (n) (n) μ2(n) -1 0 1 2 3 4 5 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 0.25 N(Rb) = 10000 (b) μ1(n) (n) (n) μ2(n) (n) -1 0 1 2 3 4 5 102 103 104 μ/ (- h ω 1) , < x> ( l 0 ), < z> ( l 0 ) N(Dy) a12 = -30a0 λ = 0.25 N(Rb) = 10000 (b) μ1(n) (n) (n) μ2(n) (n) (n) FIG. 2. (Color online) Variational (v) and numerical (n) results for chemical potential μ and rms sizes 〈x〉 and 〈z〉 for the trapped 87Rb BEC of 10 000 atoms and the bound 164Dy droplet versus the number of 164Dy atoms for (a) λ = 4 and (b) λ = 0.25. The interspecies scattering length a12 = −30a0 and the subscript 1 refer to 87Rb atoms and 2 to 164Dy atoms. and 〈z〉 of the trapped 87Rb BEC of 10 000 atoms and of the bound 164Dy droplet versus N (Dy) for (a) λ = 4 and (b) 0.25. Considering that the 164Dy droplet may not have a Gaussian shape as assumed in the variational approximation, the agreement between the numerical and variational results is quite satisfactory. The existence of the quasifree 164Dy droplet can be studied qualitatively by a minimization of the energy (19). However, the widths of the 87Rb BEC for a fixed N (Rb) do not vary much as N (Dy) or a12 is varied. Hence, for a qualitative understanding of the existence of a stable 164Dy droplet for a fixed N (Rb), we can make a further approximation in the energy (19) and take the widths of the trapped 87Rb BEC to be the same as the widths of the 87Rb BEC in the absence of 164Dy atoms under otherwise identical conditions. The widths of the 164Dy droplet so obtained for N (Rb) = 50 000 are compared in Table I for several values of a12, N (Dy), and λ with the widths obtained from exact energy minimization. We find from Table TABLE I. Root-mean-square sizes of the binary 87Rb-164Dy system of 50 000 87Rb atoms from the variational approximation (11)–(14) (Var), and an approximate energy minimization (Approx). a12 N (Dy) λ 〈x1〉 〈z1〉 〈x2〉 〈z2〉 Approx −40a0 1000 4 2.7597 0.7106 1.0755 0.4847 Var −40a0 1000 4 2.7475 0.7091 1.0703 0.4834 Approx −80a0 1000 0.25 2.2780 8.9525 0.6468 2.5799 Var −80a0 1000 0.25 2.2549 8.8558 0.6378 2.5525 Approx −60a0 500 4 2.7597 0.7106 0.8478 0.3768 Var −60a0 500 4 2.7492 0.7090 0.8443 0.3757 Approx −100a0 100 0.25 2.2780 8.9525 0.6384 1.5803 Var −100a0 100 0.25 2.2752 8.9392 0.6373 1.5779 (b) (c) FIG. 3. (Color online) The isodensity contour of (a) a disk-shaped (λ = 4) 87Rb BEC of 10 000 atoms and that of the bound 164Dy droplet of (b) 1000 atoms (a12 = −30a0), (c) 1000 atoms (a12 = −80a0), and (d) 4000 atoms (a12 = −30a0). The isodensity contour of the 87Rb BEC is practically the same in all three cases. The density on the contour in all cases is 0.001l−3 0 and length is measured in units of l0 (=1 μm). I that the size of the 87Rb for a fixed N (Rb) remains reasonably constant with respect to the variation of N (Dy) and a12 and that the approximate energy minimization with fixed widths for 87Rb BEC leads to a good approximation to the widths of the 164Dy droplet. In the numerical solution of the binary GP equations (6) and (7), we also verified that the shape and size of the trapped 87Rb BEC are fairly independent of the presence or absence of the 164Dy droplet. Hence in the study of the shape, size, and dynamics of the dipolar 164Dy droplet bound in the trapped 87Rb BEC, the trapped BEC will have only a passive role and we shall highlight only the shape, size, and dynamics of the 164Dy droplet in the following. We consider the shape of a stable 164Dy droplet in the stability plots of Fig. 1 and study its change as the collapse instability is approached. In Fig. 3(a) we show the isodensity contour plot [density |φ1(r)|2] of a disk-shaped (λ = 4) 87Rb BEC of 10 000 atoms and that [density |φ2(r)|2] of the bound 164Dy droplet in Fig. 3(b) for N (Dy) = 1000 (a12 = −30a0), in Fig. 3(c) for N (Dy) = 1000 (a12 = −80a0), and in Fig. 3(d) for N (Dy) = 4000 (a12 = −30a0). The isodensity of the 87Rb BEC is practically the same in all three cases. Of the isodensities of the 164Dy droplet, the one in Fig. 3(b) is deep inside the stability region far away from collapse instability. The parameters of Figs. 3(c) and 3(d) are close to the region of collapse instability. Of these two, Fig. 3(c) corresponds to a medium number of 164Dy atoms and Fig. 3(d) to a large number. Independent of the initial shape, a 164Dy droplet with a medium number of 164Dy atoms always collapses towards the center with shrinking size. However, a 164Dy droplet containing a large number of atoms, independent of the associated trap symmetry, always first takes a cigar shape and then collapses on the axial z direction. A strong dipolar interaction prohibits a collapse to center of a large 164Dy droplet and favors a cigar shape. In Fig. 3 the trap is disk shaped. The medium-sized 164Dy droplet of Fig. 3(c) collapses to center, maintaining the disk shape, as the net dipolar interaction is smaller in this case. The large 164Dy 013618-4 DIPOLAR DROPLET BOUND IN A TRAPPED BOSE- . . . PHYSICAL REVIEW A 87, 013618 (2013) (b) (c) FIG. 4. (Color online) The isodensity contour of (a) a cigar- shaped (λ = 0.25) 87Rb BEC of 10 000 atoms and that of the bound 164Dy droplet of (b) 1000 atoms (a12 = −30a0), (c) 1000 atoms (a12 = −65a0), and (d) 3000 atoms (a12 = −30a0). The isodensity contour of the 87Rb BEC is practically the same in all three cases. The density on the contour in all cases is 0.001l−3 0 and length is measured in units of l0 (=1 μm). droplet in Fig. 3(d) has changed its shape from the supporting disk-shaped trap to a cigar shape. If the attraction is further increased by increasing |a12|, the 164Dy droplet of Fig. 3(d) would collapse from the cigar shape on the axial z axis. In Fig. 4(a) we illustrate the isodensity contour of a cigar- shaped (λ = 0.25) 87Rb BEC of 10 000 atoms, and the same for the 164Dy droplet are shown in Fig. 4(b) for N (Dy) = 1000 (a12 = −30a0), in Fig. 4(c) for N (Dy) = 1000 (a12 = −65a0), and in Fig. 4(d) for N (Dy) = 3000 (a12 = −30a0). The density of the 87Rb BEC is virtually the same in all three cases. Of the three 164Dy droplets, the one in Fig. 4(b) is deep inside the stability region far away from collapse instability. The parameters for Figs. 4(c) and 4(d) are close to the region of collapse instability. In all cases the 164Dy droplet maintains the cigar shape of the accompanying trap acting on the 87Rb BEC. The size of the 164Dy droplet in Fig. 4(b) is larger than that of Fig. 4(c) with the same number of atoms due to a strong interspecies attraction acting on the latter. The 164Dy droplet of Fig. 4(d) is larger than those in Figs. 4(b) and 4(c) due to the larger number of 164Dy atoms in it. We investigate the dynamics of a bound 164Dy droplet in a trapped 87Rb BEC. First, to test the present scheme we consider a small 164Dy droplet of 100 atoms bound in a disk-shaped (λ = 4) 87Rb BEC of 2000 atoms with a12 = −50a0. This corresponds to the stable region of Fig. 1(a). In real-time evolution of the system, the scattering length a12 is suddenly changed to 1.01a12, thus starting the breathing oscillation. In Fig. 5 we show the resultant oscillation of the rms sizes (a) 〈ρ〉 and (b) 〈r〉 during time evolution as obtained from numerical solution of Eqs. (6) and (7) and the variational approximation (11)–(14). After obtaining a satisfactory result for the dynamics with a small system, we venture with a larger system of experimental interest, e.g., a stable 164Dy droplet of 1000 atoms bound in a cigar-shaped (λ = 0.25) 87Rb BEC of 10 000 atoms with a12 = −50a0 corresponding to the stable region of Fig. 1(b). 0 10 20 30 40 50 < r> ( ar bi tr ar y un its ) time (units of t0) λ = 4 a12 = -50a0 N(Rb) = 2000 N(Dy) = 100 (b) Dy (num) 0 10 20 30 40 50 < r> ( ar bi tr ar y un its ) time (units of t0) λ = 4 a12 = -50a0 N(Rb) = 2000 N(Dy) = 100 (b) Dy (num) Dy (var) 0 10 20 30 40 50 < ρ> ( ar bi tr ar y un its ) time (units of t0) λ = 4 a12 = -50a0 N(Rb) = 2000 N(Dy) = 100 (a) Dy (num) 0 10 20 30 40 50 < ρ> ( ar bi tr ar y un its ) time (units of t0) λ = 4 a12 = -50a0 N(Rb) = 2000 N(Dy) = 100 (a) Dy (num) Dy (var) FIG. 5. (Color online) The rms sizes (a) 〈ρ〉 and (b) 〈r〉 from numerical simulation (num) and variational approximation (var) during breathing oscillation initiated by the sudden change a12 → 1.01a12 versus time in units of t0 (=1.38 ms) of a 164Dy droplet of 100 atoms bound in a trapped 87Rb BEC of 2000 atoms for a12 = −50 and λ = 4. In real-time evolution of the system, the scattering length a12 is again suddenly changed to 1.01a12, thus starting the breathing oscillation. In Fig. 6 we show the resultant oscillation of the rms sizes (a) 〈z〉 and (b) 〈x〉 during time evolution as obtained from a numerical solution of Eqs. (6) and (7) and the variational approximation (11)–(14). The agreement between the numerical simulation and variational approximation in 0 10 20 30 40 50 < x> ( ar bi tr ar y un its ) time (units of t0)(b) λ = 0.25 a12 = -50a0 N(Rb) = 10000 N(Dy) = 1000 Dy (num) 0 10 20 30 40 50 < x> ( ar bi tr ar y un its ) time (units of t0)(b) λ = 0.25 a12 = -50a0 N(Rb) = 10000 N(Dy) = 1000 Dy (num) Dy (var) 0 10 20 30 40 50 < z> ( ar bi tr ar y un its ) time (units of t0)(a) λ = 0.25 a12 = -50a0 N(Rb) = 10000 N(Dy) = 1000 Dy (num) 0 10 20 30 40 50 < z> ( ar bi tr ar y un its ) time (units of t0)(a) λ = 0.25 a12 = -50a0 N(Rb) = 10000 N(Dy) = 1000 Dy (num) Dy (var) FIG. 6. (Color online) The rms sizes (a) 〈z〉 and (b) 〈x〉 from numerical simulation (num) and variational approximation (var) during breathing oscillation initiated by the sudden change a12 → 1.01a12 versus time in units of t0 (=1.38 ms) of a 164Dy droplet of 1000 atoms bound in a trapped 87Rb BEC of 10 000 atoms for a12 = −50 and λ = 0.25. 013618-5 LUIS E. YOUNG-S. AND S. K. ADHIKARI PHYSICAL REVIEW A 87, 013618 (2013) both cases of breathing oscillation illustrated in Figs. 5 and 6 for disk and cigar shapes is quite satisfactory considering the facts that (i) during oscillation the bound 164Dy droplet might have a density distribution different from the Gaussian distribution assumed in the variational approximation and (ii) the large nonlinearity in the 87Rb BEC will also make its density distribution non-Gaussian. The stable oscillation under small perturbation as shown in Figs. 5 and 6 confirms the dynamical stability of the 164Dy droplet for both disk- and cigar-type environments. The present quasifree dipolar droplet is not just of theoretical interest, but can be realized experimentally by initially preparing a binary 87Rb-164Dy mixture where both the components are harmonically trapped. The trap on 164Dy can then be ramped to zero exponentially during a few milliseconds. The 164Dy cloud will initially expand and finally form the quasifree 164Dy droplet. To illustrate numerically the viability of this procedure, we consider an initial 87Rb-164Dy binary mixture with the following parameters: N (Rb) = 10 000, N (Dy) = 1000, and a12 = −50a0. We take the same axial trap with λ = 4 acting on both components. Then we perform real-time propagation of the binary GP equation with this initial state and ramp down the 164Dy trap by Vtrap(Dy) → exp[−(t − 10)]Vtrap(Dy) for time t > 10. The trap on 164Dy is practically zero for times t � 20. The trapped 164Dy BEC first expands for t > 10 and eventually it emerges as the quasifree 164Dy droplet. The time evolution of the rms sizes of the 164Dy droplet is shown in Fig. 7(a). The widths of the 164Dy droplet first increase and eventually execute small 0 1 2 3 4 0 10 20 30 40 50 60 70 80 90 100 < x> , < z> ( un its o f l 0) time (units of t0) λ = 4 a12 = -50a0 N(Rb) = 10000 N(Dy) = 1000 (a) (Dy) 0 1 2 3 4 0 10 20 30 40 50 60 70 80 90 100 < x> , < z> ( un its o f l 0) time (units of t0) λ = 4 a12 = -50a0 N(Rb) = 10000 N(Dy) = 1000 (a) (Dy) (Dy) 0 1 2 3 4 0 10 20 30 40 50 60 70 80 90 100 < x> , < z> ( un its o f l 0) time (units of t0) λ = 4 a12 = -50a0 N(Rb) = 10000 N(Dy) = 1000 (a) (Rb) 0 1 2 3 4 0 10 20 30 40 50 60 70 80 90 100 < x> , < z> ( un its o f l 0) time (units of t0) λ = 4 a12 = -50a0 N(Rb) = 10000 N(Dy) = 1000 (a) (Rb) (Rb) 0 0.2 0.4 0.6 -4 -2 0 2 4 (units of l0) (c) |ϕ(x)|2 t = 20 0 0.2 0.4 0.6 -4 -2 0 2 4 (units of l0) (c) |ϕ(x)|2 t = 20 40 0 0.2 0.4 0.6 -4 -2 0 2 4 (units of l0) (c) |ϕ(x)|2 t = 20 40 60 0 0.2 0.4 0.6 -4 -2 0 2 4 (units of l0) (c) |ϕ(x)|2 t = 20 40 60 80 0 0.2 0.4 0.6 -4 -2 0 2 4 (units of l0) (c) |ϕ(x)|2 t = 20 40 60 80 100 0 0.2 0.4 0.6 -4 -2 0 2 4 (units of l0) (c) |ϕ(x)|2 t = 20 40 60 80 100 sta 0 0.2 0.4 0.6 -4 -2 0 2 4 |ϕ |2 ( un its o f l 0-1 ) (units of l0) (b) |ϕ(z)|2 t = 20 0 0.2 0.4 0.6 -4 -2 0 2 4 |ϕ |2 ( un its o f l 0-1 ) (units of l0) (b) |ϕ(z)|2 t = 20 40 0 0.2 0.4 0.6 -4 -2 0 2 4 |ϕ |2 ( un its o f l 0-1 ) (units of l0) (b) |ϕ(z)|2 t = 20 40 60 0 0.2 0.4 0.6 -4 -2 0 2 4 |ϕ |2 ( un its o f l 0-1 ) (units of l0) (b) |ϕ(z)|2 t = 20 40 60 80 0 0.2 0.4 0.6 -4 -2 0 2 4 |ϕ |2 ( un its o f l 0-1 ) (units of l0) (b) |ϕ(z)|2 t = 20 40 60 80 100 0 0.2 0.4 0.6 -4 -2 0 2 4 |ϕ |2 ( un its o f l 0-1 ) (units of l0) (b) |ϕ(z)|2 t = 20 40 60 80 100 sta FIG. 7. (Color online) (a) The rms sizes 〈x〉 and 〈z〉 of the 164Dy and 87Rb BECs versus time during the passage of the trapped 164Dy BEC in the binary 87Rb-164Dy mixture with N (Rb) = 10 000, N (Dy) = 1000, a12 = −50a0, and λ = 4 as the confining trap on the 164Dy BEC is relaxed for t > 10 exponentially as Vtrap → exp[−(t − 10)]Vtrap, so that this trap is practically zero for t > 20. The linear 1D densities (b) |ϕ(z)|2 and (c) |ϕ(x)|2 at times t = 20,40,60,80,100 (in units of t0) during the expansion together with the numerically calculated density of the stationary (sta) 164Dy droplet. oscillations illustrating the stable 164Dy droplet. To investigate the shape of the 164Dy droplet, we plot in Figs. 7(b) and 7(c) the one-dimensional densities |ϕ(z)|2 ≡ ∫∫ dxdy|φ(x,y,z)|2 and |ϕ(x)|2 ≡ ∫∫ dydz|φ(x,y,z)|2 at times t = 20, 40, 60, 80, and 100 along with the corresponding numerical densities of the stationary droplet calculated from the solution of Eqs. (6) and (7). The 1D numerical densities of the stationary 164Dy droplet are in agreement with the densities obtained from dynamical simulation of the passage of the trapped 164Dy BEC to a quasifree droplet. In Figs. 7(b) and 7(c) one can easily identify the small oscillation of the evolving dynamical droplet around its stationary shape. IV. SUMMARY AND DISCUSSION Using variational approximation and numerical solution of a set of coupled mean-field GP equations, we demonstrate the existence of a stable dipolar 164Dy droplet bound by interspecies attraction in a trapped nondipolar BEC of 87Rb atoms. The domain of stability of the 164Dy droplet is highlighted in stability plots of the number of 164Dy atoms and interspecies scattering length a12 for both cigar- and disk- shaped traps acting on the 87Rb BEC. Results of the variational approximation and numerical solution for the statics (sizes and chemical potentials) and dynamics (breathing oscillation) of the 164Dy droplet are found to be in satisfactory agreement with each other. We also demonstrate numerically that such droplets can be obtained experimentally by considering a trapped binary 87Rb-164Dy BEC and then removing the trap on 164Dy. The 164Dy BEC then expands and transforms into a bound dipolar droplet. Dipolar interaction among atoms is quite different from the normal short-range atomic interaction and manifests itself in different ways in a trapped dipolar BEC. However, in a trapped dipolar BEC, the confining constraints could be too strong and could make the effect of dipolar interaction difficult to observe. A special experimental set up and theoretical formulation might be necessary to study the effect of dipolar interaction. However, it will be much easier to see the effect of dipolar interaction in the present binary mixture of nondipolar 87Rb and dipolar 164Dy. Strong dipolar interaction in the axial polarization z direction should elongate the dipolar BEC along this direction, thus transforming it to a cigar shape. However, it will be difficult to observe this change in the presence of a harmonic trap along z. The effect of dipolar interaction is clearly seen in the present quasifree 164Dy droplet of Fig. 3(d), where due to a strong dipolar interaction the shape of the 164Dy droplet has changed from the disk shape–the shape of the 87Rb BEC of Fig. 3(a) responsible for its binding–to the cigar shape. The thin disk-shaped 87Rb BEC stays only near the central z = 0 plane of the 164Dy droplet. Most of the 164Dy atoms lying outside the 87Rb BEC are bound due to the intraspecies dipolar interaction. In the dipolar 164Dy droplet of Fig. 3(d) the dipolar interaction is playing a more important role in its binding and shape. The present quasifree dipolar droplet may also find other interesting applications in the BEC phenomenology. Among the interesting features in a trapped dipolar BEC, one can mention its peculiar shape and stability properties [11], D- wave collapse [19], anisotropic solitons, vortex solitons [20] 013618-6 DIPOLAR DROPLET BOUND IN A TRAPPED BOSE- . . . PHYSICAL REVIEW A 87, 013618 (2013) and vortex lattices [21], anisotropic shock- and sound-wave propagation [22], anisotropic Landau critical velocity [23], stable checkerboard, stripe, and star configurations in a two- dimensional optical lattice such as a stable Mott insulator [24], as well as superfluid soliton [25] states. It would be of great interest to find out how these features and properties of a trapped dipolar BEC would manifest in a 3D quasifree dipolar droplet. For example, a quasifree dipolar droplet could be used in the experimental study of anisotropic sound- and shock-wave propagation [22], collapse dynamics [19], anisotropic Landau critical velocity [23], formation of vortex dipoles and vortex lattices [21], etc., in a different setting of confinement which will facilitate the observation of the effect of the anisotropic dipolar interaction. ACKNOWLEDGMENTS We thank FAPESP and CNPq (Brazil) for partial support. [1] T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, and T. Pfau, Rep. Prog. Phys. 72, 126401 (2009). [2] T. Lahaye, T. Koch, B. Fröhlich, M. Fattori, J. Metz, A. Griesmaier, S. Giovanazzi, and T. Pfau, Nature (London) 448, 672 (2007); A. Griesmaier, J. Stuhler, T. Koch, M. Fattori, T. Pfau, and S. Giovanazzi, Phys. Rev. Lett. 97, 250402 (2006). [3] J. Stuhler, A. Griesmaier, T. Koch, M. Fattori, T. Pfau, S. Giovanazzi, P. Pedri, and L. Santos, Phys. Rev. Lett. 95, 150406 (2005). [4] K. Goral, K. Rzazewski, and T. Pfau, Phys. Rev. A 61, 051601 (2000). [5] T. Koch et al., Nat. Phys. 4, 218 (2008). [6] 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). [7] M. Lu, N. Q. Burdick, Seo Ho Youn, and B. L. Lev, Phys. Rev. Lett. 107, 190401 (2011). [8] K. Aikawa, A. Frisch, M. Mark, S. Baier, A. Rietzler, R. Grimm, and F. Ferlaino, Phys. Rev. Lett. 108, 210401 (2012). [9] J. Deiglmayr, A. Grochola, M. Repp, K. Mortlbauer, C. Gluck, J. Lange, O. Dulieu, R. Wester, and M. Weidemuller, Phys. Rev. Lett. 101, 133004 (2008); M. H. G. de Miranda, A. Chotia, B. Neyenhuis, D. Wang, G. Quéméner, S. Ospelkaus, J. L. Bohn, J. Ye, and D. S. Jin, Nat. Phys. 7, 502 (2011); J. Doyle, B. Friedrich, R. V. Krems, and F. Masnou-Seeuws, Eur. Phys. J. D 31, 149 (2004). [10] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 71, 463 (1999). [11] N. G. Parker, C. Ticknor, A. M. Martin, and D. H. J. O’Dell, Phys. Rev. A 79, 013617 (2009); R. M. Wilson, S. Ronen, and J. L. Bohn, ibid. 80, 023614 (2009); N. G. Parker and D. H. J. O’Dell, ibid. 78, 041601 (2008); L. Santos, G. V. Shlyapnikov, P. Zoller, and M. Lewenstein, Phys. Rev. Lett. 85, 1791 (2000); C. Ticknor, N. G. Parker, A. Melatos, S. L. Cornish, D. H. J. O’Dell, and A. M. Martin, Phys. Rev. A 78, 061607 (2008); R. M. W. van Bijnen, A. J. Dow, D. H. J. O’Dell, N. G. Parker, and A. M. Martin, ibid. 80, 033617 (2009); A. Junginger, J. Main, G. Wunner, and T. Bartsch, ibid. 86, 023632 (2012); M. Asad-uz-Zaman and D. Blume, ibid. 80, 053622 (2009); S. Ronen, D. C. E. Bortolotti, and J. L. Bohn, Phys. Rev. Lett. 98, 030406 (2007). [12] 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). [13] Luis E. Young-S. and S. K. Adhikari, Phys. Rev. A 86, 063611 (2012). [14] S. Yi and L. You, Phys. Rev. A 63, 053607 (2001). [15] K. Goral and L. Santos, Phys. Rev. A. 66, 023613 (2002). [16] V. M. Perez-Garcia, H. Michinel, J. I. Cirac, M. Lewenstein, and P. Zoller, Phys. Rev. A 56, 1424 (1997). [17] 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). [18] S. Inouye et al., Nature (London) 392, 151 (1998). [19] T. Lahaye, J. Metz, B. Frohlich, T. Koch, M. Meister, A. Griesmaier, T. Pfau, H. Saito, Y. Kawaguchi, and M. Ueda, Phys. Rev. Lett. 101, 080401 (2008). [20] I. Tikhonenkov, B. A. Malomed, and A. Vardi, Phys. Rev. Lett. 100, 090406 (2008); S. K. Adhikari and P. Muruganandam, J. Phys. B 45, 045301 (2012); Luis E. Young-S., P. Muruganandam, and S. K. Adhikari, ibid. 44, 101001 (2011); P. Muruganandam and S. K. Adhikari, ibid. 44, 121001 (2011); P. Köberle, D. Zajec, G. Wunner, and B. A. Malomed, Phys. Rev. A 85, 023630 (2012). [21] R. M. W. van Bijnen, D. H. J. O’Dell, N. G. Parker, and A. M. Martin, Phys. Rev. Lett. 98, 150401 (2007); R. Kishor Kumar and P. Muruganandam, J. Phys. B 45, 215301 (2012); M. Abad, M. Guilleumas, R. Mayol, M. Pi, and D. M. Jezek, Phys. Rev. A 79, 063622 (2009). [22] P. Muruganandam and S. K. Adhikari, Phys. Lett. A 376, 480 (2012); C. Krumnow and A. Pelster, Phys. Rev. A 84, 021608 (2011). [23] R. M. Wilson, S. Ronen, and J. L. Bohn, Phys. Rev. Lett. 104, 094501 (2010). [24] B. Capogrosso-Sansone, C. Trefzger, M. Lewenstein, P. Zoller, and G. Pupillo, Phys. Rev. Lett. 104, 125301 (2010). [25] K. Lakomy, R. Nath, and L. Santos, Phys. Rev. A 85, 033618 (2012); S. K. Adhikari and P. Muruganandam, Phys. Lett. A 376, 2200 (2012). 013618-7 http://dx.doi.org/10.1088/0034-4885/72/12/126401 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.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.1038/nphys887 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/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.108.210401 http://dx.doi.org/10.1103/PhysRevLett.101.133004 http://dx.doi.org/10.1103/PhysRevLett.101.133004 http://dx.doi.org/10.1038/nphys1939 http://dx.doi.org/10.1140/epjd/e2004-00151-x http://dx.doi.org/10.1140/epjd/e2004-00151-x http://dx.doi.org/10.1103/RevModPhys.71.463 http://dx.doi.org/10.1103/RevModPhys.71.463 http://dx.doi.org/10.1103/PhysRevA.79.013617 http://dx.doi.org/10.1103/PhysRevA.80.023614 http://dx.doi.org/10.1103/PhysRevA.78.041601 http://dx.doi.org/10.1103/PhysRevLett.85.1791 http://dx.doi.org/10.1103/PhysRevA.78.061607 http://dx.doi.org/10.1103/PhysRevA.80.033617 http://dx.doi.org/10.1103/PhysRevA.86.023632 http://dx.doi.org/10.1103/PhysRevA.80.053622 http://dx.doi.org/10.1103/PhysRevLett.98.030406 http://dx.doi.org/10.1103/PhysRevLett.98.030406 http://dx.doi.org/10.1103/PhysRevA.86.033606 http://dx.doi.org/10.1103/PhysRevLett.102.230403 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.63.053607 http://dx.doi.org/10.1103/PhysRevA.66.023613 http://dx.doi.org/10.1103/PhysRevA.56.1424 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.1038/32354 http://dx.doi.org/10.1103/PhysRevLett.101.080401 http://dx.doi.org/10.1103/PhysRevLett.100.090406 http://dx.doi.org/10.1103/PhysRevLett.100.090406 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/44/10/101001 http://dx.doi.org/10.1088/0953-4075/44/12/121001 http://dx.doi.org/10.1103/PhysRevA.85.023630 http://dx.doi.org/10.1103/PhysRevA.85.023630 http://dx.doi.org/10.1103/PhysRevLett.98.150401 http://dx.doi.org/10.1088/0953-4075/45/21/215301 http://dx.doi.org/10.1103/PhysRevA.79.063622 http://dx.doi.org/10.1103/PhysRevA.79.063622 http://dx.doi.org/10.1016/j.physleta.2011.11.054 http://dx.doi.org/10.1016/j.physleta.2011.11.054 http://dx.doi.org/10.1103/PhysRevA.84.021608 http://dx.doi.org/10.1103/PhysRevA.84.021608 http://dx.doi.org/10.1103/PhysRevLett.104.094501 http://dx.doi.org/10.1103/PhysRevLett.104.094501 http://dx.doi.org/10.1103/PhysRevLett.104.125301 http://dx.doi.org/10.1103/PhysRevA.85.033618 http://dx.doi.org/10.1103/PhysRevA.85.033618 http://dx.doi.org/10.1016/j.physleta.2012.05.030 http://dx.doi.org/10.1016/j.physleta.2012.05.030