PHYSICAL REVIEW A 87, 063621 (2013) Bright solitons in quasi-one-dimensional dipolar condensates with spatially modulated interactions F. Kh. Abdullaev,1 A. Gammal,2 B. A. Malomed,3 and Lauro Tomio1,4 1Instituto de Fı́sica Teórica, Universidade Estadual Paulista, 01140-070, São Paulo, São Paulo, Brazil 2Instituto de Fı́sica, Universidade de São Paulo, 05508-090, São Paulo, São Paulo, Brazil 3Department of Physical Electronics, School of Electrical Engineering, Faculty of Engineering, Tel Aviv University, Tel Aviv 69978, Israel 4Centro de Ciências Naturais e Humanas, Universidade Federal do ABC, 09210-170, Santo André, Brazil (Received 15 March 2013; revised manuscript received 28 May 2013; published 20 June 2013) We introduce a model for the condensate of dipolar atoms or molecules, in which the dipole-dipole interaction (DDI) is periodically modulated in space due to a periodic change of the local orientation of the permanent dipoles, imposed by the corresponding structure of an external field (the necessary field can be created, in particular, by means of magnetic lattices, which are available to the experiment). The system represents a realization of a nonlocal nonlinear lattice, which has a potential to support various spatial modes. By means of numerical methods and variational approximation (VA), we construct bright one-dimensional solitons in this system and study their stability. In most cases, the VA provides good accuracy and correctly predicts the stability by means of the Vakhitov-Kolokolov criterion. It is found that the periodic modulation may destroy some solitons, which exist in the usual setting with unmodulated DDI and can create stable solitons in other cases, not verified in the absence of modulations. Unstable solitons typically transform into persistent localized breathers. The solitons are often mobile, with inelastic collisions between them leading to oscillating localized modes. DOI: 10.1103/PhysRevA.87.063621 PACS number(s): 67.85.Hj, 03.75.Kk, 03.75.Lm I. INTRODUCTION Studies of Bose-Einstein condensates (BEC) in optical lattices (OLs) are some of the central topics in the current work on cold atoms [1,2]. Two types of the spatially periodic modulation of the parameters of BEC are possible under the action of OLs. First is a periodic (lattice) potential for atoms, which represents the linear OL. The spatial periodicity of another kind can be realized in BEC via a modulation of the atomic scattering length, i.e., the effective local nonlinearity of the condensate, which is known as the nonlinear OL (or a magnetic lattice [3], if the scattering length is affected by the magnetic field). New types of solitons and solitary vortices are possible in nonlinear lattices [4]. Both linear and nonlinear OLs can be realized in conden- sates with atoms interacting via long-range dipolar-dipolar interactions (DDI). Due to the anisotropic and nonlocal character of the DDI, many remarkable new phenomena have been predicted and observed in these ultracold gases, such as new quantum phases, anisotropic collapse and suppression of the collapse, and specific modes of collective excitations [5,6], as well as stabilization of modulated patterns trapped in the OL potential, which are unstable in the BEC with the contact interaction [7]. The DDI may also give rise to localized modes in the form of stable two-dimensional (2D) [8] and one-dimensional (1D) [9,10] solitons and soliton complexes [11]. In fact, the presence of the linear OL potential strongly facilitates the creation and stabilization of solitons in dipolar condensates [9,12–14]. A new type of lattice structure can be created in dipolar condensates, by means of a spatially periodic modulation of the strength of the DDI. In fact, these are nonlocal nonlinear lattices, as the DDI represents long-range cubic interactions. Nonlocal nonlinear lattices of different types are possible too in optical media with thermal nonlinearity [15]. In the dipolar BEC, one possibility for the creation of such a lattice is to use a condensate of atoms or molecules with permanent magnetic or electric dipolar moments, and apply a spatially nonuniform (misaligned) polarizing field, under the action of which the mutual orientation of the moments varies in space periodically. Another approach may make use of a condensate of polarizable atoms or molecules, in which the moments are induced by the external field [16], whose strength varies periodically along the system. In either case, the necessary periodically nonuniform external magnetic fields can be induced by the currently available magnetic lattices [3]. If the condensate is formed by particles with permanent or induced electric moments, the corresponding electric field can be supplied by similarly designed ferroelectric lattices. The electrostriction effects on DDI induced by a far-off-resonant standing laser field can also lead to a periodic structure in the dipolar BEC [17]. Such systems are described by the Gross-Pitaevskii (GP) equation with a spatially varying nonlocal nonlinearity. The objective of the present work is to investigate the existence and stability of effectively 1D [18,19] solitons in the dipolar BEC under the action of nonlocal nonlinear lattices. Here we focus on the model which assumes permanent moments with a periodically varying orientation. The soliton modes will be studied by means of a variational approximation (VA), and with the help of direct numerical simulations of the respective GP equation. For direct simulations, we have employed a relaxation algorithm applied to the GP equation, as described in Ref. [20]. The paper is organized as follows. The model is introduced in Sec. II. Analytical results (the VA, and some simplification of the full model) are reported in Sec. III. Numerical results are summarized and compared to the variational predictions in Sec. IV, and the paper is concluded by Sec. V. 063621-11050-2947/2013/87(6)/063621(12) ©2013 American Physical Society http://dx.doi.org/10.1103/PhysRevA.87.063621 ABDULLAEV, GAMMAL, MALOMED, AND TOMIO PHYSICAL REVIEW A 87, 063621 (2013) II. MODEL Under the assumption that that the external field imposes an orientation of permanent magnetic moments periodically varying along the x axis, the corresponding dipolar BEC obeys the following nonlocal GP equation [9,19,21]: ih̄ ∂� ∂t + h̄2 2m ∂2� ∂x2 − 2a(r) s h̄ω⊥|�|2� − 2d2 l3 ⊥ � ∫ +∞ −∞ dξ f (x,ξ )R(|x − ξ |) |�(ξ )|2 = 0. (1) Here, a(r) s is the renormalized s-wave scattering length, which includes a contribution from the DDI a(r) s = as [ 1 + εDD 2 (1 − 3 cos2 θ0) ] , with εDD = md2 3h̄2as ≡ ad 3as , (2) where as is the proper scattering length, which characterizes collisions between the particles, θ0 is the average value of the dipole’s orientation angle, and ad , defined by the dipole momentum d and the atomic mass m, plays the role of the effective DDI scattering length. In the above expressions, ω⊥ is the transverse trapping frequency, which defines the respective trapping radius l⊥ = √ h̄/mω⊥. The reduction of the full GP equation to its 1D form (2) is valid provided that the peak density of the condensate n0 is small enough to make the non- linear term much weaker than the transverse confinement, and prevent the excitation of transverse models by the DDIs, i.e.,∫ d3r ′Udd (�r − �r ′)|�(�r ′,t)|2 � h̄ω⊥, Udd (�r − �r ′) = d2[1 − 3 cos2(θ0)] 4πε0|�r − �r ′|3 , where ε0 is the vacuum permittivity. In this case, the 3D wave function may be factorized as �(�r,t) ≈ �(x,t)R(ρ), with the transverse radial mode R(ρ) taken as the ground state of the confining harmonic-oscillator potential. After performing the integration over the transverse radius ρ in the GP equation, we arrive at the quasi-1D equation (1). The exact DDI kernel R(x) appearing in Eq. (1) has a rather cumbersome form [19,21], but it may be approximated reasonably well by [9] R(x) = ε3 (x2 + ε2)3/2 , (3) where a cutoff parameter ε, which truncates the formal singularity, is on the order of the transverse radius of the cigar-shaped trap l⊥. In the following, the variables are rescaled as x → l⊥x, ξ → l⊥x ′, t → t/ω⊥, 2ar s / l⊥ → g, �(x,t) → ψ(x,t)/ √ l⊥, (4) with the corresponding dimensionless cutoff parameter fixed to be ε = 1. Further, the function f (x,ξ ) in Eq. (1) describes the modulation of the DDI due to the variation of the dipoles’ orientation along the x axis, f (x,ξ ) = cos [θ (x) − θ (ξ )] − 3 cos [θ (x)] cos [θ (ξ )] , (5) where we assume that the local angle of the orientation of the dipoles with respect to the x axis is periodically modulated as follows: θ (x) = θ0 + θ1 cos (kx) . (6) The periodic variation of θ (x) can be induced by magnetic lattices (MLs), which have been constructed on the basis of ferromagnetic films with permanent periodic structures built into them, the condensate being loaded into the ML [3]. Accordingly, the ground state of the dipolar BEC is formed by an elongated trap with a strong magnetic field B, whose relevant value may be estimated as �200 G. The periodic modulation of the DDI is then induced by the ML with the amplitude of the respective component of the magnetic field B ∼ 40–50 G and a modulation period λ ∼ 0.5 μm, which is within the same order of the magnitude as a typical OL period. The regular ML can be constructed to include ∼1000 sites, which is quite sufficient to observe the localized modes described below [3]. With regard to the rescaling (4), the combined kernel in Eq. (1) is transformed as 2d2/(l3 ⊥h̄ω⊥) f (x,ξ )R(|x − ξ |) → GVDD(x,x ′), where G ≡ 2d2/(l3 ⊥h̄ω⊥) is the DDI effective strength, and VDD(x,x ′) = cos[θ (x) − θ (x ′)] − 3 cos[θ (x)] cos[θ (x ′)] [(x − x ′)2 + ε2]3/2 ≡ − cos[θ (x) − θ (x ′]) + 3 cos[θ (x) + θ (x ′)] 2[(x − x ′)2 + ε2]3/2 . (7) Thus, the underlying GP equation is cast into the form of the following partial-integral differential equation, which combines the contact and dipole-dipole interactions i ∂ψ(x) ∂t = −1 2 ∂2ψ(x) ∂x2 + g|ψ(x)|2ψ(x) +Gψ(x) ∫ +∞ −∞ VDD(x,x ′)|ψ(x ′)|2dx ′. (8) The effective kernel (7) can be further simplified for specific settings. As an example, for the average orientation angles θ0 = πn/2, with n = 0 or n = 1 (the mean orientations parallel or perpendicular to the axis, respectively), we have V (θ0=πn/2) DD (x,x ′) = − cos{θ1[cos(kx) − cos(kx ′)]} + 3(−1)n cos{θ1[cos(kx) + cos(kx ′)]} 2[(x − x ′)2 + ε2]3/2 . (9) Another point of interest is the one when the interaction reduces to zero for θ = θ0. This happens at θ0 ≡ θm ≡ arccos(1/ √ 3) ≈ 0.9553, (10) 063621-2 BRIGHT SOLITONS IN QUASI-ONE-DIMENSIONAL . . . PHYSICAL REVIEW A 87, 063621 (2013) with the respective kernel being V (θ0=θm) DD (x,x ′) ≡ − sin[θ1 cos(kx)] sin[θ1 cos(kx ′)] − √ 2 sin{θ1[cos(kx) + cos(kx ′)]} [(x − x ′)2 + ε2]3/2 . (11) Characteristic shapes of the kernel for four different configura- tions, defined by Eqs. (9) and (11), are displayed in Fig. 1. We choose here some specific configurations of interest. First, in Fig. 1(a) we show the unmodulated case (θ1 = 0), ranging from the most attractive (θ0 = 0) to the most repulsive (θ0 = π/2) conditions. In this case, the maximum of |V (x,y)| is reached when x = y. In view of the symmetry of V (x,y), only positive values for (x − y) are presented. In the other three panels, we fix y = 0. In Fig. 1(b) we show how this interaction varies with θ1, by fixing θ0 = arccos(1/ √ 3), so that the interaction vanishes [as also shown in Fig. 1(a)] at θ1 = 0. In the lower two panels, we plot the effective DDI kernel for the two cases of interest, by varying θ1 with θ0 = 0 [Fig. 1(c)] and θ0 = π/2 [Fig. 1(d)]. Thus, the effective (pseudo)potential of the nonlinear interactions in Eq. (8) is Veff(x; |ψ |2) = g|ψ(x)|2 + G ∫ +∞ −∞ VDD(x,x ′)|ψ(x ′)|2dx ′. (12) 0 0.5 1.0 1.5 2.0 |x − y| −2 −1 0 1 V (x ,y ) π/2 0.9553 π/2−0.9553 0 θ1=0 (a) θ0= 0 0.5 1.0 1.5 2.0 x −2 −1 0 1 2 V (x ,0 ) 0 0.9553 π/2 3π/2 2π−0.9553 θ0=0.9553 θ1= (b) 0 0.5 1.0 1.5 2.0 x −2 −1 0 1 2 V (x ,0 ) 0 1.57 3.14 3.8 4.71 6.28 θ0= 0 θ1= (c) 0 0.5 1.0 1.5 2.0 x −2 −1 0 1 2 V (x ,0 ) 0 1.0 1.57 3.14 4.71 6.28 θ0= π/2 θ1= (d) FIG. 1. (Color online) The effective DDI kernel, as given by Eqs. (7) and (6). The kernel for θ0 = θm = 0.9553 (in the upper panels) corresponds to the more specific expression (11), and the ones for θ0 = 0 and θ = π/2 (in the lower panels) correspond to Eq. (9). In view of the symmetry of V (x,y), we choose to fix y = 0 and take x > 0, except in the case of (a) the unmodulated DDI, where V (x,y) is given as a function of |x − y|. 063621-3 ABDULLAEV, GAMMAL, MALOMED, AND TOMIO PHYSICAL REVIEW A 87, 063621 (2013) Coefficients G, g, ε, and k in Eq. (8) can be easily rescaled. In view of that, in the following we set G ≡ 1, ε ≡ 1, and also choose k = 2π , which fixes the respective period in Eq. (6) as 2π/k = 1. As a starting point, we set g ≡ 0 (no effective local interaction). Then, the remaining parameters are coefficients θ0 and θ1 of periodic modulation (6), and the total norm of the solution N = ∫ +∞ −∞ |ψ(x)|2 dx. (13) The objective is to find a family of bright-soliton solutions in this model, and analyze their stability. In the system with the unmodulated DDI, which also included the linear OL potential, this problem was considered by the authors of Ref. [9]. A related model (the Tonks-Girardeau gas with the attractive DDI) was considered by the authors of Ref. [22]. Bright solitons may exist if the interaction is self-attractive, i.e., the sign of the averaged value of the (pseudo)potential (12) is negative. We start from the case of θ1 = 0 (or k = 0, with θ0 redefined as θ1 + θ0), when the periodic modulation is absent. The respective region of values of θ0, where solitons are likely to occur, can be readily identified. For this case, the interaction kernel VDD(x,y) is shown in Fig. 1(a). The corresponding numerator of the expression (7) is (1 − 3 cos2 θ0), leading to attractive interactions for all values of θ0 with cos2(θ0) � 1/3. Within the period of 0 < θ0 < 2π , this condition holds for 0 < θ0 < θm, π − θm < θ0 < π + θm and 2π − θm < θ0 < 2π [recall θm is given by Eq. (10)]. Next, by switching on θ1, the interaction kernel is made undulating, and the existence of solitons turns out to be possible outside of these intervals, while in some areas inside the intervals the solitons do not exist, as we demonstrate below. A. Approximate analytical approach A simplified version of the model can be introduced by using the well-known Fourier expansions which generate the Bessel functions cos[z cos(η)] = J0(z) + 2 ∞∑ n=1 (−1)nJ2n(z) cos(2nη), (14) sin[z cos(η)] = 2 ∞∑ n=0 (−1)nJ2n+1 cos[(2n + 1)η]. By keeping the lowest-order harmonic in the expansions, we arrive at the following approximation for interaction kernel (7) with the modulation format (6): VDD(x,x ′) ≈ { sin2 θ0 [ J 2 0 (θ1) − 8J 2 1 (θ1) cos(kx) cos(kx ′) ] + 3 sin(2θ0)J0(θ1)J1(θ1)[cos(kx) + cos(kx ′)] + 2 cos2(θ0) [ 2J 2 1 (θ1) cos(kx) cos(kx ′) − J 2 0 (θ1) ]} × [(x − x ′)2 + ε2]−3/2. (15) In the particular cases of the dipoles oriented parallel (θ0 = 0) and perpendicular (θ0 = π/2) to the x axis, the expression (15) reduces, respectively, to V (θ0=0) DD (x,x ′) ≈ 2 [ 2J 2 1 (θ1) cos(kx) cos(kx ′) − J 2 0 (θ1) ] × [(x − x ′)2 + ε2]−3/2, (16) V (θ0=π/2) DD (x,x ′) ≈ [ J 2 0 (θ1) − 8J 2 1 (θ1) cos(kx) cos(kx ′) ] × [(x − x ′)2 + ε2]−3/2. (17) The existence condition can be found in an explicit form for broad solitons, whose width is much larger than that of the kernel (i.e., the GP equation becomes a quasilocal one). In this case, |ψ |2 appearing in the integral of the nonlocal term of the GP equation may be replaced by a constant, hence the necessary condition for the existence of a bright soliton amounts to ∫ +∞ −∞ VDD(0,x ′)dx ′ < 0. (18) For the analysis of this condition, we take into account that∫ +∞ −∞ 1 [(x ′)2 + ε2]3/2 dx ′ = 2 ε2 , (19)∫ +∞ −∞ cos(kx ′) [(x ′)2 + ε2]3/2 dx ′ = 2|kε| ε2 K1(|kε|), where K1(kε) is the modified Bessel functions. Making use of the symmetry properties of Eqs. (6) and (7) and the periodicity, in the following we assume that θ0 and θ1 are positive. To identify existence regions for broad solitons, we consider three illustrative cases. (1) For the dipoles oriented along the x axis, i.e., with θ0 = 0, intervals of θ1 for the existence of the bright soliton follow from Eqs. (16) and (18) in the form of 2|kε| K1(|kε|) J 2 1 (θ1) − J 2 0 (θ1) < 0. (20) Taking, as said above, ε = 1 and k = 2π , we have 2|kε| K1(|kε|) = 0.0124, hence the first three existence in- tervals are 0 < θ1 < 2.29, 2.52 < θ1 < 5.40, 5.63 < θ1 < 8.54. (21) In each interval, the maximum strength of the attractive DDI is attained, severally, at θ1 = 0, 3.83, and 7.0. Virtually the same existence intervals are produced by numerical solutions (see Sec. IV below). (2) For dipoles oriented perpendicular to the x axis, condition (20) is replaced by J 2 0 (θ1) − 8|kε| K1(|kε|) J 2 1 (θ1) < 0, (22) so that, for ε = 1 and k = 2π [8|kε| K1(|kε|) = 0.0496], the first three existence intervals for the broad solitons are 2.17 < θ1 < 2.62, 5.30 < θ1 < 5.73, 8.43 < θ1 < 8.87, (23) cf. Eq. (21), with the strongest DDI attraction attained at θ1 = 2.38,5.51, and 8.65, respectively. (3) In the case of Eq. (10), when the unmodulated DDI vanishes, the expansion converges slowly, so a numerical analysis of the full integral expression in the GP equation (8) 063621-4 BRIGHT SOLITONS IN QUASI-ONE-DIMENSIONAL . . . PHYSICAL REVIEW A 87, 063621 (2013) is necessary. In this case, the solitons exist in much larger intervals of θ1 than in the two previous cases. The first two intervals are 0.61 < θ1 < 3.83, 4.16 < θ1 < 7.01. (24) Further, the consideration of the three above cases reveals overlapping regions where more than one soliton may exist 2.17 < θ1 < 2.29; 2.52 < θ1 < 2.62; 5.30 < θ1 < 5.40; 5.63 < θ1 < 5.73; 8.43 < θ1 < 8.54. (25) A comparison of these predictions with results of the numerical solution of the GP equation is given in Sec. IV. B. Variational approximation To apply the VA, we look for stationary solutions to Eq. (8), with chemical potential μ, as ψ = exp(−iμt)φ(x), where φ(x) satisfies the equation μφ = −1 2 d2φ(x) dx2 + g|φ(x)|2φ(x) +Gφ(x) ∫ +∞ −∞ VDD(x,x ′)|φ(x ′)|2dx ′, (26) which can be derived from the Lagrangian L = ∫ +∞ −∞ Ldx, (27) where density L is given by L = μ|φ|2 − 1 2 ∣∣∣∣dφ dx ∣∣∣∣ 2 − g 2 |φ|4 − G 2 |φ(x)|2 ∫ +∞ −∞ VDD(x,x ′)|φ(x ′)|2dx ′. (28) For a bright-soliton solution, we adopt the simplest Gaussian variational ansatz, with amplitude A and width a φ = A exp ( − x2 2a2 ) , (29) whose norm (13) is N = √ πA2a. The substitution of the ansatz into Lagrangian (27) with density (28) yields the corresponding averaged Lagrangian, which is expressed in terms of N , instead of the amplitude A L̄ = N [ μ − 1 4a2 − gN 2 √ 2πa − GN 2πa2 F (a,k,V0,V1) ] , (30) where we define F (a,k,V0,V1) ≡ ∫ +∞ −∞ e−x2/a2 dx ∫ +∞ −∞ e−x ′2/a2 VDD(x,x ′)dx ′ (31) [recall that VDD(x,x ′) is defined by Eqs. (7) and (6)]. Lastly, the Euler-Lagrange equations dL̄/dN = dL̄/da = 0 are derived from the averaged Lagrangian (30) μ = 1 4a2 + gN√ 2πa + GNF πa2 , (32) N = √ 2π −ga + √ 2/πG (a∂F/∂a − 2F ) . III. NUMERICAL RESULTS In this section, we present the results of the numerical solu- tions of the GP equation (8) for soliton modes, which are com- pared to predictions based on the variational equations (32). A. Dipoles oriented, on average, along the x axis (θ0 = 0) Numerically obtained solitary-mode density profiles for θ1 = 3.8317, which corresponds to a strong attractive DDI in the case of θ0 = 0, are displayed in Fig. 2 for two values of the chemical potential: μ = −1 and −0.3. The solutions were constructed using the full kernel (7), the Bessel approximation (15), as well as the variational ansatz given by Eq. (29) with the amplitude and width found from a numerical solution of algebraic equations (32). A good agreement is observed between the solutions of all the three types. In Fig. 3, two panels are displayed for θ0 = 0, with the chemical potential given as a function of the norm for the orientation of the dipoles parallel to the x axis, for different values of the the modulation amplitude θ1 in Eq. (6). The Vakhitov-Kolokolov (VK) stability criterion dμ/dN < 0 [23] (in particular, the application of the VK criterion to nonlinear OLs was developed by the authors of Ref. [24]) suggests that almost all the localized modes are stable in the region of 0 < θ1 < 3.6, except for a small region of parameter θ1. The unstable region, given by 2.3 < θ1 < 2.5, can be extracted from the results presented in Figs. 3(a) and 3(b). In Fig. 3(a), for θ1 = 2.3, we observe a change of the slope indicating the loss of stability for μ < −0.2. Accordingly, in Fig. 3(b), this change of slope can be observed for θ1 = 2.5, implying the loss of the stability at μ < −0.1. In both limiting cases, we have N ≈ 65. The VA is, in general, consistent in predicting the stability by means of the VK criterion, even considering that the agreement of the predicted shape of the solutions with −8 −6 −4 −2 0 2 4 6 8 x 0.0 0.4 0.8 1.2 1.6 2.0 |ψ| −1.0 −1.0 (A) −1.0 (v) −0.3 −0.3 (A) −0.3 (v) μ = θ0=0 θ1=3.8317 FIG. 2. (Color online) Soliton profiles for dipoles oriented, on the average, parallel to the x axis (θ0 = 0), for θ1 = 3.8317, with μ = −1.0 and −0.3. Full numerical solutions (solid lines, with solid squares pertaining to μ = −0.3) are compared to the approximate solutions (A) obtained from Eq. (16) (dashed lines, with symbols plus indicating μ = −0.3) and with the variational approximation (v). 063621-5 ABDULLAEV, GAMMAL, MALOMED, AND TOMIO PHYSICAL REVIEW A 87, 063621 (2013) 0 10 20 30 40 50 60 70 N −1 −0.8 −0.6 −0.4 −0.2 0 μ 1.6 2.0 2.2 2.3 1.6 (v) 2.0 (v) 2.2 (v) θ0=0 θ1= (a) 0 10 20 30 40 50 60 70 N −1 −0.8 −0.6 −0.4 −0.2 0 μ 2.5 2.6 3.0 3.6 3.6 (v) 3.0 (v) 2.6 (v) θ0=0 θ1= (b) FIG. 3. (Color online) The chemical potential versus the soliton’s norm, for θ0 = 0 with several values of θ1, from 1.6 to 3.6. The solutions are displayed in the intervals where the stability is predicted by Eq. (21). The corresponding variational results, indicated by (v), are also displayed, except in two cases (close to the unstable regions of the parameters), where we cannot obtain consistent results from the VA. their numerical counterparts becomes poor for large values of N and μ. We also note that, close to the region where we have unstable full numerical results (for large μ) we are not able to reach any conclusion by means of the VA, therefore only full numerical results are presented for θ1 = 2.3 and 2.6. In general, for the stable cases, the comparison to the predictions of the VA shows that the agreement is rather good, with the discrepancy amounting to small shifts in the values of N . Following the results presented in Figs. 2 and 3 for θ0 = 0, at larger values of this angle 0 < θ0 < π/4, we have observed a similar stable behavior, produced by numerical solutions and the variational analysis. In general, it is observed that, by increasing the mean angle θ0 from zero, the existence region of stable solitons shrinks. In the next two subsections we consider two other cases of interest. B. Dipoles oriented, on average, perpendicular to the x axis (θ0 = π/2) In the case of θ0 = π/2, when the dipoles are oriented, on average, perpendicular to the x axis, stable solitons are found only in a small interval of parameter θ1. Indeed, by using the VK criterion, one can conclude that almost all values of θ1 correspond to unstable states, except for θ1 smaller than ≈1.5. −4 −2 0 2 4 x 0 1 2 3 4 |ψ| (exact) N=11 (A) (a) (v) θ0=π/2 θ1=2.4045 μ = −1.0 −6 −4 −2 0 2 4 6 x 0 1 2 3 4 5 6 |ψ| (exact) N=16 (A) (v) N=35 (v) N=75 θ0=π/2 θ1=2.4045 μ = −0.2 (b) FIG. 4. (Color online) Soliton profiles for the dipoles oriented, on the average, perpendicular to the x axis (θ0 = π/2 ), with μ = −1 and (a) N ≈ 11 and (b) μ = −0.2. The full-numerical results are shown by solid (blue) lines. The results with the approximate kernel (17), shown by squares, are in perfect agreement with exact numerical results. The VA results (v) are shown by dashed (red) lines and up-triangles. In panel (b), we also show, by small dashed (green) lines and down-triangles, that another VA solution exists, such that we have two possible values for the norm, N ≈ 35 and N ≈ 75 as indicated, both being far from the exact one, which is N ≈ 16. 063621-6 BRIGHT SOLITONS IN QUASI-ONE-DIMENSIONAL . . . PHYSICAL REVIEW A 87, 063621 (2013) −2 −1 0 1 2 x 0 2 4 6 8 |Ψ| (exact) (v) θ0=π/2, θ1=1.0 μ=−0.5 (a) 0 1 2 3 4 x −1 −0.5 0 0.5 1 V ef f(x ) Veff(exact) Veff(v) θ0=π/2 θ1=1.0 μ=−0.5 (b) FIG. 5. (Color online) The soliton profile, for θ0 = π/2, θ1 = 1.0, and μ = −0.5, is shown in the panel (a), with the corresponding effective potential Veff (x), as defined by Eq. (12) with g = 0 and G = 1, shown in panel (b). Numerically exact results are given by solid (black) lines. The corresponding variational (v) calculations are represented by dashed (red) lines with circles. In Fig. 4, we display two panels for soliton profiles, with chemical potentials μ = −1 [Fig. 4(a)] and μ − 0.2 [Fig. 4(b)], for θ1 = 2.408, which corresponds to the strongest attractive DDI. A good agreement is observed between the results produced by the approximate kernel (17) (where the Bessel expansion is used) and the full numerical solutions. As concerns the comparison of the numerically exact results with the VA predictions, the Gaussian ansatz (29) may be too simple in this case to reproduce the exact results. Actually, the Gaussian profile of the VA follows the exact solution approximately in the case of μ = −1.0, shown in Fig. 4(a). But, for μ = −0.2, two solutions with different values of N are predicted by the VA for the same μ, both being very inaccurate, as shown in Fig. 4(b). The effective (pseudo)potential Veff(x), which is defined by Eq. (12), plays an important role in the analysis of the soliton modes. To display the characteristic profiles of the potential, in Fig. 5 we plot the soliton profile [Fig. 5(a)] and the corresponding effective potential [Fig. 5(b)] for θ1 = 1.0 and θ0 = π/2, with μ = −0.5. In this case, the solutions are stable, featuring a very close agreement between the numerical and variational results. In Fig. 6, we present three panels with numerical results for the dependence between the chemical potential and norm at several values of θ1. In terms of the VK criterion, within the interval of μ presented here −1 < μ < 0, the stability of the numerical solutions is quite clear only for θ1 = 1.0, as shown in Fig. 6(a). As we increase θ1, in Fig. 6(b) we see that the system is already at the stability limit for θ1 = 1.5, and becomes unstable for larger values of this parameter. The stability predictions of the VA, which are also displayed in the figure, are in reasonable agreement with the numerical results for larger values of |μ|, where both the VA and full numerical solutions produce similar slopes. This conclusion is also supported by the perfect agreement between the variational and numerical results in Fig. 5. For other solutions, presented in Figs. 6(b) and 6(c), the VK criterion does not predict clear stability regions. The 24 25 26 27 N −1 −0.8 −0.6 −0.4 −0.2 0 μ (exact) (v) (a) θ0=π/2 θ1=1.0 5 10 15 20 N 1.5 2.0 2.5 1.5 (v) 2.0 (v) 2.5 (v) θ0=π/2 θ1= (b) 26 30 34 38 42 46 N 3.0 4.5 4.5 (v) (c) θ0=π/2 θ1= FIG. 6. (Color online) The chemical potential versus the norm for θ0 = π/2 and several values of θ1. The variational results, indicated by label (v), provide a reasonable agreement with the numerical results as |μ| increases, and only for values of θ1 which do not belong to the interval of 2.5 < θ1 < 4. In panel (c), for θ1 = 3.0, there is no VA solution in the displayed interval of μ. 063621-7 ABDULLAEV, GAMMAL, MALOMED, AND TOMIO PHYSICAL REVIEW A 87, 063621 (2013) 0 10 20 30 40 t 0.0 0.5 1.0 1.5 2.0 2.5 |ψ (0 ,t) |2 /N μ=−0.3 (a) −4 −3 −2 −1 0 1 2 3 4 x 0.0 0.5 1.0 1.5 2.0 |ψ (x ,t) |2 /N t=0 t=3.6 μ=−0.3 (b) FIG. 7. (Color online) A typical example of the evolution of (a) the peak density in the case when unstable solitons do not decay, featuring instead oscillations between two different shapes, as shown in panel (b). Here, θ0 = π/2 and θ1 = 2.4048. agreement between the VA and full numerical solutions is completely lost in the interval of 2.5 < θ1 < 4, where no reasonable VA solutions were found for −1 < μ < 0. For the sake of the comparison with the VA, we keep the plot with θ1 = 2.5 in Fig. 6(b), which shows how the agreement between the VA and numerical data starts to deteriorate. Therefore, for θ1 = 3.0, we can only show the exact numerical results. Another point to be observed in these plots is that a conspicuous region of the reasonable agreement of the VA with numerical results (although without stabilization of the 0 10 20 30 40 N −1 −0.8 −0.6 −0.4 −0.2 0 μ 0 −0.1 −0.15 −0.25 −0.15(v) −0.25(v) g = θ1=3.0 θ0=π/2 FIG. 8. (Color online) The role of the local self-attractive nonlin- earity, accounted for by coefficient g < 0, in stabilizing the solitons at θ0 = π/2 and θ1 = 3.0 (which are shown in Fig. 6 to be unstable at g = 0). As suggested by the VK criterion, and is corroborated by numerical simulations (not shown here in detail), the solutions are stabilized at g < −0.15. The corresponding VA-predicted μ(N ) branches are shown for g = −0.15 and −0.25. For smaller |g|, the simple Gaussian ansatz (29) cannot provide reliable results. solitons) reappears at large values of θ1, as seen in Fig. 6(c) for θ1 = 4.5. Concerning the evolution of unstable solitons, there are generic situations in which the solitons are not exactly stable, but they do not decay either. This happens, for example, at θ0 = π/2 and θ1 = 2.4048, as shown in Fig. 7, where the profiles are presented in Fig. 6(b) for two instants of time, with μ = −0.3. In Fig. 6(a) we show the peak of the density at the origin (x = 0). In this case, the unstable soliton is transformed into a persistent breather oscillating between two different 4 8 12 16 20 24 28 32 N −1.0 −0.8 −0.6 −0.4 −0.2 μ 1.0 2.0 2.5 3.0 3.5 4.0 4.5 θ0=0.9553 θ1= FIG. 9. (Color online) The chemical potential versus the norm, for θ0 = θm [see Eq. (10)], at several values of θ1. In agreement with analytical predictions for the region of the soliton existence, see Eq. (24), in most cases we have stable results, except for the small region of 2 < θ1 < 3, where the results were represented by dashed (red) lines. The respective variational results are chiefly consistent, in terms of the stability prediction, with the above, but not shown, as they yield much larger values of N . 063621-8 BRIGHT SOLITONS IN QUASI-ONE-DIMENSIONAL . . . PHYSICAL REVIEW A 87, 063621 (2013) −4 −2 0 2 4 6 8 x, t 0.0 0.1 0.2 0.3 0.4 |ψ (x ,t) |2 /N (a) 0 1 2 3 4 x −0.4 −0.2 0.0 0.2 0.4 V ef f(x ) Veff(exact) Veff(v) θ0=0 θ1=3.8 μ=−0.5 (b) FIG. 10. (Color online) Numerical results for a soliton moving at velocity vs = 1 are shown in panel (a), by means of the juxtaposition of density profiles, with the horizontal axis simultaneously showing values of x and t . The parameters are θ0 = 0, θ1 = 3.8, and μ = −0.5. In panel (b) we have the corresponding effective (pseudo)potential. configurations. While such breathers are characteristic modes of the present system with the spatially modulated DDI, they have not been observed in models of the dipolar BEC with the constant DDI. Next, we proceed to the consideration of the GP equation (8) which includes the contact interactions, accounted for by g = 0. The results are presented in Fig. 8 in terms of the μ(N ) curves, which show, by means of the VK criterion (whose validity has been tested in this case by means of direct simulations), the role of this contact interaction in stabilizing solitons that were found to be unstable for g = 0, with θ0 = π/2. To this end, we here consider the case of θ1 = 3.0. The corresponding results for g = 0 are presented above in Fig. 6(c), and are included here too for the sake of the comparison. It is observed that the solutions get stabilized by the self-attractive contact nonlinearity with g < −0.15. The results are reported here for full numerical solutions, except in the cases where stable solutions were found, where we also consider the corresponding VA (see the cases of g = −0.15 and −0.25). At smaller |g|, the VA predictions are far from the numerical results, at least for large values of μ. A large discrepancy between the variational and numerical curves’ behavior is actually observed in Fig. 8 for g = −0.15. C. Average dipoles’ orientation corresponding to θ0 = θm ≡ arccos(1/ √ 3) It is also interesting to consider the configuration with the mean orientation of the dipoles corresponding to the zero interaction (if θ1 = 0), as per Eq. (10). The corresponding dependence of μ on N is plotted in Fig. 9, for several values of θ1, from 1.0 to 4.5. In agreement with the predictions of the VA, in this case most soliton solutions are stable, considering condition (24). However, the VA results yield much larger values of N , and it also predicts stable solitons in the region where the numerical solutions are unstable, e.g., for θ1 = 2.5. -20 -15 -10 -5 0 5 10 15 20 x 0 -10 -20 t |ψ(x,t)| 2 /N (a) -30 -20 -10 0 10 20 30 x 0 10 20 30 t |ψ(x,t)| 2 /N (b) FIG. 11. (Color online) (Color on-line) The two-soliton collision at the same parameters as in Fig. 10. Originally, the solitons’ centers are placed at points x = ±20, as shown in the panel (a), with the collision point being x = 0. The collision gives rise to internal excitations of the solitons, which manifest themselves as periodic oscillations of their amplitudes, observed in the panel (b). 063621-9 ABDULLAEV, GAMMAL, MALOMED, AND TOMIO PHYSICAL REVIEW A 87, 063621 (2013) Summarizing our numerical analysis, exemplified by the above subsections for θ0 = 0, θ0 = π/2, and θ0 = 0.9553, we conclude that the case of θ0 = 0 is the one where the stability is more likely to be reached for a wide range of values of the modulation parameter θ1 (see Fig. 3). On the other hand, the case of θ0 = π/2 is the one where the stability may occur only in a very specific region of values of θ1. D. Soliton motion and collisions By means of numerical simulations, we have also in- vestigated the mobility of solitons in the framework of the present model. As usual, the soliton was set in motion, multiplying its wave function by the kicking factor exp (ivsx). For a configuration with θ0 = 0, θ1 = 3.8, and μ = −0.5, and by taking vs = 1, our simulations demonstrate the soliton propagation with a practically undistorted shape, as shown in Fig. 10(a). Radiation losses are practically absent in this case, as well as no excitation of internal modes is observed in the soliton. In Fig. 10(b), the effective potential, as given by Eq. (12), is plotted for the same parameters. It gives rise to an effective Peierls-Nabarro potential experienced by the moving soliton [11,25]. For a broad soliton, whose width is large in comparison to the internal scale of the potential (this is the case in Fig. 10), the effective obstacle created by the potential is exponentially small, hence one may expect practically free motion of the soliton, which is confirmed by the simulations, even for an initial speed reduced to vs = 0.05. Additionally, at such a low speed the soliton propagates along the lattice without changing its shape. Concluding the numerical analysis of the moving solitons, in Fig. 11 the collision of two solitons is displayed for the same set of parameters as given in Fig. 10. Figures 11(a) and 11(b) display the evolution of the density profiles, before and after the collision, respectively. As one observes in Fig. 11(b), the collision does not destroy the solitons, but rather induces fluctuations in the shapes of the colliding solitons (probably −3 −2 −1 0 1 2 3 x t 5.8 5.6 5.4 5.2 5.0 4.8 4.6 4.4 4.2 4.0 |Ψ(x,t)|2/N 4.0 < t < 5.8 −1 0 1 2 3 x+(t−13) 13.0 13.2 13.4 13.6 13.8 14.0 14.2 14.4 14.6 14.8 t 13.0 < t < 14.8|Ψ(x,t)| 2 /N 4 3 2 1 0 -1 -2 -3 -4 0 2 4 6 8 10 t x 12 14 4.5 3.5 2.5 1.5 0.5 0 1 2 3 4 FIG. 12. (Color online) Illustrated in the upper panel, we present the evolution of density profiles |ψ(x,t)|2/N , for the collision of two solitons, for θ0 = 0, θ1 = 3.8, with chemical potential μ = −2.5. The initial velocities and positions of the solitons are vs ± 0.3 and x = ∓3, respectively. The juxtapositions of density profiles for their time evolution are shown in the bottom rows, considering two time intervals, as shown inside the frames. The formation of a bound state trapped around x = 0, along with the strong emission of radiation, is clearly observed as the result of this inelastic collision. 063621-10 BRIGHT SOLITONS IN QUASI-ONE-DIMENSIONAL . . . PHYSICAL REVIEW A 87, 063621 (2013) −15 −10 −5 0 5 10 15 x 0 1 2 3 |ψ (x ,t )| 2 /N μ=−5 θ0=0 θ1=3.8 −> −> −> −> <− <− <− <− t = 0 2 4 6 8 8 6 4 2 0 0 < t < 10 (a) −5 0 5 x 0 1 2 3 4 5 |Ψ (x ,t )| 2 /N 10 < t < 14 (b) FIG. 13. The inelastic collision of two solitons for θ0 = 0, θ1 = 3.8, with chemical potentials μ = −5 and initial positions x = ±10. Panel (a) shows a sequence of soliton profiles at 0 � t � 10, through time intervals �t = 2. In panel (b), a sequence of soliton profiles is shown at 10 � t � 14. The merger of the colliding solitons into a bound state and strong emission of radiation are clearly identified by the plots. related to the excitation of internal modes), eventually showing an almost elastic collision. It is obvious that the present system is far from any integrable limit, hence inelastic collisions leading to a merger of the colliding solitons are expected as well. Two examples of such strongly inelastic collisions are shown in Figs. 12 and 13, for the parameters θ0 = 0 and θ1 = 3.8. In Fig. 12, the results are plotted for two solitons with μ = −2.5 and initial velocity vs = 0.3. The merger of the solitons into an oscillating localized wave packet is observed here as the result of the collision, accompanied by the emission of radiation. In Fig. 13, we present results for the collision of two narrow solitons with μ = −5. Figures 13(a) and 13(b) display the picture before and after the collision, for time intervals 0 � t � 10 and 10 � t � 14, respectively. The merger into an oscillating localized mode, surrounded by a cloud of emitted radiation, is clearly seen in this case also. IV. CONCLUSION In this work, we propose the realization of a nonlinear nonlocal lattice in quasi-one-dimensional dipolar BECs. The lattice is introduced by a spatially periodic modulation of the DDI, that, in turn, may be induced by an external field, which periodically changes the local orientation of permanent dipoles. The necessary setting may be created in experiments using the currently available MLs. We have focused on the consideration of bright solitons, which were constructed in the numerical form, and by means of the VA. The stability of the solitons is studied by direct simulations of their perturbed evolution and also by means of the VK criterion [23]. The VA provides, overall, a good accuracy in comparison with numerical results. Stable solitons exist below a critical value of the angle (θ0) between the average orientation of the dipoles and the system’s axis. The stability area shrinks with the increase of the modulation amplitude (θ1). Unstable solitons are, typically, spontaneously transformed into persistent localized breathers (which are not observed in the system with the uniform DDI). In addition, the mobility of the solitons and collisions between them have been studied. Strongly inelastic collisions lead to the merger of the solitons into oscillatory localized modes, surrounded by conspicuous radiation fields. The next objective, which will be considered elsewhere, is to study the model of the condensate composed of polarizable atoms or molecules, in which the magnitude of the external polarizing field, which induces the local electric or magnetic moments (in the absence of a permanent moment), varies periodically along the system’s axis. A challenging possibility is to extend this class of models to two-dimensional settings. ACKNOWLEDGMENTS B.A.M. appreciates support through a visitor’s grant from the South American Institute for Fundamental Research, and the hospitality of the Instituto de Fı́sica Teórica at the Universidade Estadual Paulista (São Paulo, Brazil). We also acknowledge support from Brazilian agencies Fundação de Amparo à Pesquisa do Estado de São Paulo and Conselho Nacional de Desenvolvimento Cientı́fico e Tecnológico. [1] M. Greiner, O. Mandel, T. Esslinger, T. W. Hansch, and I. Bloch, Nature (London) 415, 39 (2002); D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. 81, 3108 (1998); V. A. Brazhnyi and V. V. Konotop, Mod. Phys. Lett. B 18, 627 (2004); O. Morsch and M. Oberthaler, Rev. Mod. Phys. 78, 179 (2006); M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A. Sen(De), and U. Sen, Adv. Phys. 56, 243 (2007); I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 929 (2008). [2] M. Lewenstein, A. Sanpera, and V. Ahufinger, Ultracold Atoms in Optical Lattices: Simulating Quantum Many-Body Systems (Oxford University Press, Oxford, England, 2012). 063621-11 http://dx.doi.org/10.1038/415039a http://dx.doi.org/10.1103/PhysRevLett.81.3108 http://dx.doi.org/10.1142/S0217984904007190 http://dx.doi.org/10.1142/S0217984904007190 http://dx.doi.org/10.1103/RevModPhys.78.179 http://dx.doi.org/10.1103/RevModPhys.78.179 http://dx.doi.org/10.1080/00018730701223200 http://dx.doi.org/10.1103/RevModPhys.80.885 http://dx.doi.org/10.1103/RevModPhys.80.885 ABDULLAEV, GAMMAL, MALOMED, AND TOMIO PHYSICAL REVIEW A 87, 063621 (2013) [3] S. Ghanbari, T. D. Kieu, A. Sidorov, and P. Hannaford, J. Phys. B 39, 847 (2006); A. Abdelrahman, P. Hannaford, and K. Alameh, Opt. Express 17, 24358 (2009). [4] Y. V. Kartashov, B. A. Malomed, and L. Torner, Rev. Mod. Phys. 83, 247 (2011). [5] M. Baranov, Phys. Rep. 464, 71 (2008). [6] T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, and T. Pfau, Rep. Prog. Phys. 72, 126401 (2005). [7] A. Maluckov, G. Gligorić, Lj. Hadžievski, B. A. Malomed, and T. Pfau, Phys. Rev. Lett. 108, 140402 (2012); Phys. Rev. A 87, 023623 (2013). [8] P. Pedri and L. Santos, Phys. Rev. Lett. 95, 200404 (2005); R. Nath, P. Pedri, and L. Santos, Phys. Rev. A 76, 013606 (2007); I. Tikhonenkov, B. A. Malomed, and A. Vardi, Phys. Rev. Lett. 100, 090406 (2008); Phys. Rev. A 78, 043614 (2008). [9] J. Cuevas, B. A. Malomed, P. G. Kevrekidis, and D. J. Frantzeskakis, Phys. Rev. A 79, 053608 (2009). [10] L. E. Young-S., P. Muruganandam, and S. K. Adhikari, J. Phys. B 44, 101001 (2011). [11] K. Lakomy, R. Nath, and L. Santos, Phys. Rev. A 86, 013610 (2012). [12] G. Gligorić, A. Maluckov, L. Hadžievski, and B. A. Malomed, Phys. Rev. A 78, 063615 (2008); 79, 053609 (2009); J. Phys. B: At. Mol. Opt. Phys. 42, 145302 (2009). [13] P. Muruganandam, R. K. Kumar, and S. K. Adhikari, J. Phys. B: At. Mol. Opt. Phys. 43, 205305 (2010); S. K. Adhikari and P. Muruganandam, ibid. 44, 121001 (2011); 45, 045301 (2012); Phys. Lett. A 376, 2200 (2012). [14] S. Rojas-Rojas, R. A. Vicencio, M. I. Molina, and F. Kh. Abdullaev, Phys. Rev. A 84, 033621 (2011). [15] Y. V. Kartashov, V. A. Vysloukh, and L. Torner, Opt. Lett. 33, 1774 (2008). [16] S. Yi and L. You, Phys. Rev. A 61, 041604 (2000); B. Deb and L. You, ibid. 64, 022717 (2001); T. J. McCarthy, M. T. Timko, and D. R. Herschbach, J. Chem. Phys. 125, 133501 (2006); Z. D. Li, Q. Y. Li, P. B. He, J. Q. Liang, W. M. Liu, and G. S. Fu, Phys. Rev. A 81, 015602 (2010); A. E. Golomedov, G. E. Astrakharchik, and Yu. E. Lozovik, ibid. 84, 033615 (2011). [17] S. Giovanazzi, D. O’Dell, and G. Kurizki, Phys. Rev. Lett. 88, 130402 (2002). [18] M. Fattori, G. Roati, B. Deissler, C. D’Errico, M. Zaccanti, M. Jona-Lasinio, L. Santos, M. Inguscio, and G. Modugno, Phys. Rev. Lett. 101, 190405 (2008). [19] S. Sinha and L. Santos, Phys. Rev. Lett. 99, 140406 (2007). [20] M. Brtka, A. Gammal, and L. Tomio, Phys. Lett. A 359, 339 (2006). [21] Y. Cai, M. Rosenkranz, Z. Lei, and W. Bao, Phys. Rev. A 82, 043623 (2010). [22] B. B. Baizakov, F. Kh. Abdullaev, B. A. Malomed, and M. Salerno, J. Phys. B: At. Mol. Opt. Phys. 42, 175302 (2009). [23] M. Vakhitov and A. Kolokolov, Radiophys. Quantum Electron. 16, 783 (1973); L. Bergé, Phys. Rep. 303, 259 (1998). [24] F. Kh. Abdullaev, A. Gammal, M. Salerno, and L. Tomio, Phys. Rev. A 77, 023615 (2008). [25] F. Kh. Abdullaev, A. Gammal, and L. Tomio, J. Phys. B 37, 635 (2004); F. Kh. Abdullaev, R. M. Galimzyanov, M. Brtka, and L. Tomio, Phys. Rev. E 79, 056220 (2009). 063621-12 http://dx.doi.org/10.1088/0953-4075/39/4/009 http://dx.doi.org/10.1088/0953-4075/39/4/009 http://dx.doi.org/10.1364/OE.17.024358 http://dx.doi.org/10.1103/RevModPhys.83.247 http://dx.doi.org/10.1103/RevModPhys.83.247 http://dx.doi.org/10.1016/j.physrep.2008.04.007 http://dx.doi.org/10.1088/0034-4885/72/12/126401 http://dx.doi.org/10.1103/PhysRevLett.108.140402 http://dx.doi.org/10.1103/PhysRevA.87.023623 http://dx.doi.org/10.1103/PhysRevA.87.023623 http://dx.doi.org/10.1103/PhysRevLett.95.200404 http://dx.doi.org/10.1103/PhysRevA.76.013606 http://dx.doi.org/10.1103/PhysRevA.76.013606 http://dx.doi.org/10.1103/PhysRevLett.100.090406 http://dx.doi.org/10.1103/PhysRevLett.100.090406 http://dx.doi.org/10.1103/PhysRevA.78.043614 http://dx.doi.org/10.1103/PhysRevA.78.043614 http://dx.doi.org/10.1103/PhysRevA.79.053608 http://dx.doi.org/10.1088/0953-4075/44/10/101001 http://dx.doi.org/10.1088/0953-4075/44/10/101001 http://dx.doi.org/10.1103/PhysRevA.86.013610 http://dx.doi.org/10.1103/PhysRevA.86.013610 http://dx.doi.org/10.1103/PhysRevA.78.063615 http://dx.doi.org/10.1103/PhysRevA.79.053609 http://dx.doi.org/10.1088/0953-4075/42/14/145302 http://dx.doi.org/10.1088/0953-4075/42/14/145302 http://dx.doi.org/10.1088/0953-4075/43/20/205305 http://dx.doi.org/10.1088/0953-4075/43/20/205305 http://dx.doi.org/10.1088/0953-4075/44/7/075301 http://dx.doi.org/10.1088/0953-4075/45/4/045301 http://dx.doi.org/10.1016/j.physleta.2012.05.030 http://dx.doi.org/10.1103/PhysRevA.84.033621 http://dx.doi.org/10.1364/OL.33.001774 http://dx.doi.org/10.1364/OL.33.001774 http://dx.doi.org/10.1103/PhysRevA.61.041604 http://dx.doi.org/10.1103/PhysRevA.64.022717 http://dx.doi.org/10.1063/1.2202829 http://dx.doi.org/10.1103/PhysRevA.81.015602 http://dx.doi.org/10.1103/PhysRevA.84.033615 http://dx.doi.org/10.1103/PhysRevLett.88.130402 http://dx.doi.org/10.1103/PhysRevLett.88.130402 http://dx.doi.org/10.1103/PhysRevLett.101.190405 http://dx.doi.org/10.1103/PhysRevLett.101.190405 http://dx.doi.org/10.1103/PhysRevLett.99.140406 http://dx.doi.org/10.1016/j.physleta.2006.05.067 http://dx.doi.org/10.1016/j.physleta.2006.05.067 http://dx.doi.org/10.1103/PhysRevA.82.043623 http://dx.doi.org/10.1103/PhysRevA.82.043623 http://dx.doi.org/10.1088/0953-4075/42/17/175302 http://dx.doi.org/10.1007/BF01031343 http://dx.doi.org/10.1007/BF01031343 http://dx.doi.org/10.1016/S0370-1573(97)00092-6 http://dx.doi.org/10.1103/PhysRevA.77.023615 http://dx.doi.org/10.1103/PhysRevA.77.023615 http://dx.doi.org/10.1088/0953-4075/37/3/009 http://dx.doi.org/10.1088/0953-4075/37/3/009 http://dx.doi.org/10.1103/PhysRevE.79.056220