# Two-dimensional discrete solitons in dipolar Bose-Einstein condensates

###### Abstract

We analyze the formation and dynamics of bright unstaggered solitons in the disk-shaped dipolar Bose-Einstein condensate, which features the interplay of contact (collisional) and long-range dipole-dipole (DD) interactions between atoms. The condensate is assumed to be trapped in a strong optical-lattice potential in the disk’s plane, hence it may be approximated by a two-dimensional (2D) discrete model, which includes the on-site nonlinearity and cubic long-range (DD) interactions between sites of the lattice. We consider two such models, that differ by the form of the on-site nonlinearity, represented by the usual cubic term, or more accurate nonpolynomial one, derived from the underlying 3D Gross-Pitaevskii equation. Similar results are obtained for both models. The analysis is focused on effects of the DD interaction on fundamental localized modes in the lattice (2D discrete solitons). The repulsive isotropic DD nonlinearity extends the existence and stability regions of the fundamental solitons. New families of on-site, inter-site and hybrid solitons, built on top of a finite background, are found as a result of the interplay of the isotropic repulsive DD interaction and attractive contact nonlinearity. By themselves, these solutions are unstable, but they evolve into robust breathers which exist on an oscillating background. In the presence of the repulsive contact interactions, fundamental localized modes exist if the DD interaction (attractive isotropic or anisotropic) is strong enough. They are stable in narrow regions close to the anticontinuum limit, while unstable solitons evolve into breathers. In the latter case, the presence of the background is immaterial.

###### pacs:

03.75.Lm; 05.45.Yv## I Introduction

Disk-shaped two-dimensional (2D) Bose-Einstein condensates (BECs) have been created using a superposition of a tight optical trap, formed by a pair of parallel strongly repelling (blue-detuned) light sheets, and a loose in-plane radial magnetic trap refer1 (). Such configurations may exist with both self-repulsive and self-attractive collisional nonlinearity, which can be controlled by means of the Feshbach-resonance technique. In the low-density limit, the dynamics of the 2D condensate is modeled by the 2D Gross-Pitaevskii equation (GPE) with the cubic nonlinearity. If the density is not very small, a consistent reduction of the underlying three-dimensional cubic GPE to two dimensions leads to the 2D nonpolynomial nonlinear Schrödinger equation (NPSE) luca (); luca2dnovi () (a related but different approach leads to another form of the NPSEs for 1D and 2D configurations, in the case of the self-repulsive nonlinearity Canary ()).

A strong 2D periodic optical-lattice (OL) potential, applied to the disk-shaped BEC, splits it into a planar array of droplets coupled by weak nearest-neighbor linear interactions, due to the tunneling of atoms across potential barriers which separate the droplets ChaosReview (). In this situation, the BEC can be modeled by the 2D discrete GPE, derived from its continuous counterpart by means of an expansion based on localized Wannier functions wannier (); pgk2dfund (). The same discretization procedure, if applied to the continuous equation with the nonpolynomial nonlinearity, leads to the discrete variant of the NPSE, as recently demonstrated in detail in the 1D geometry Luca+we (). A noteworthy feature of the NPSE models, both continuous and discrete ones, is that, in the case of the attractive local nonlinearity, they account for the onset of the collapse even in the framework of the 1D description (recall that the GPE with the self-focusing cubic term does not give rise to the collapse in 1D).

The BEC dynamics is strongly affected by long-range interactions in dipolar condensates, which can be formed by chromium atoms in an external magnetic field, as shown in experiments Cr (). In particular, the dipole-dipole (DD) attraction in the condensate gives rise to a specific (-wave) mode of the collapse in the condensate d-collapse (). Nevertheless, the dipolar condensate can be stabilized against the collapse, adjusting the scattering length of the collisional interactions via the Feshbach-resonance technique experim1 (). It was also proposed to create a condensate dominated by DD interactions between electric dipole moments, that may be induced in atoms by an external electric field dc (). A similar situation may be expected in BEC formed by dipolar molecules hetmol (). In particular, the creation of LiCs dipolar molecules in an ultracold gas mixture has been reported LiCs ().

Once dipolar condensates are available to the experiment, a natural issue is how the nonlocal DD interactions affects localized modes (solitons). Various possibilities to create 2D solitons in dipolar condensates were explored theoretically. It was demonstrated that isotropic solitons Pedri05 () and vortices Ami2 () may exist in the 2D disk-shaped configuration with the moments polarized perpendicular to the disk’s plane, if the sign of the DD interaction is inverted by means of the rapid rotation of dipoles reversal (). The ordinary DD interaction (with the uninverted sign) can give rise to stable anisotropic 2D solitons, if the dipoles are polarized in the plane of the confining disk Ami1 (). It may be relevant to mention that stable 2D solitons supported by nonlocal interactions (different form those corresponding to the DD forces) are known in optical media with the thermal nonlinearity Krolik ().

Solitons in 1D models of the dipolar BEC have been recently predicted too, using the respective GPE Cuevas () or NPSE we-latest (). In the latter case, the model enables to study the influence of the long-range DD nonlinearity on the onset of the collapse in the condensate with the self-attractive nonlinearity. Stable 1D solitons were also predicted in the model of the Tonks-Girardeau (TG) gas of bosons carrying permanent moments, using the known model of the TG gas with the self-repulsive quintic local term, to which the long-range DD attraction was added BBB ().

A natural extension of the theoretical studies of localized modes in the dipolar BEC is to include a strong OL potential, which leads to the discrete model with the DD interactions between remote lattice sites. Recently, this analysis was performed in the 1D setting gpe (); Santos (); NPSE (). In particular, unstaggered solitons in such models, with both the cubic and nonpolynomial on-site nonlinearities, were studied, respectively, in Refs. gpe () and NPSE (), with a conclusion that the DD interactions might enhance the stability of the discrete solitons.

The objective of the present work is to introduce discrete models for the dipolar condensate trapped in a deep 2D lattice, and analyze the existence and stability of localized modes in these settings (discrete bright solitons). We will focus on the most fundamental case of unstaggered solitons, performing the analysis in terms of two different lattice models. The first model includes the ordinary cubic on-site nonlinearity, i.e., it corresponds to the 2D version of the cubic discrete GPE (to be abbreviated as DGPE), while the second model deals with the nonpolynomial local nonlinearity, which corresponds to the 2D discrete NPSE (alias DNPSE). Both models incorporate the cubic nonlocal term accounting for the DD interactions. The results obtained for discrete solitons in the frameworks of both models will be compared.

The paper is structured as follows. The two discrete models are formulated in Section II. Families of on-site, inter-site and hybrid 2D solitons in the case of the attractive contact interaction, combined with isotropic and anisotropic DD interactions, are presented in Section III. Their stability is analyzed in the same section. Properties of fundamental localized structures that may be formed as a result of the competition between the repulsive local nonlinearity and attractive DD isotropic interactions (i.e., with the inverted sign of the DD force, as mentioned above), as well as between the repulsive local and anisotropic DD interactions with the natural sign, are considered in Section IV. The paper is concluded by Section V.

## Ii The models

The scaled form of the continuous GPE, which includes both contact and DD interactions in the 3D geometry, is well known Cr (); Pedri05 ():

(1) |

where is the external potential, is the normalized strength of the local interactions, the number of bosonic atoms, the scattering length of inter-atomic collisions, and is the scale of the trapping provided by potential , which acts perpendicular to plane of the disk-shaped condensate. The total potential, , includes both the latter term and the in-plane 2D OL potential, . The notation used in Eq. (1) assumes that wave function is subject to the usual normalization,

(2) |

Coefficient in Eq. (1) defines the strength and sign of the DD interactions between atomic dipoles, and is the angle between vector and the orientation of dipoles, which is fixed by a strong external magnetic field.

The reduction of Eq. (1) to the 2D form is performed by assuming the factorization of the 3D wave function luca (); luca2dnovi ():

(3) |

where is the width of the distribution along axis ( may depend on coordinates and time ), is an arbitrary in-plane (2D) wave function, and is its normalized axial counterpart,

(4) |

With regard to the underlying 3D normalization condition (2), the factorized ansatz based on Eqs. (3) and (4) yields the following normalization condition for the 2D wave function,

(5) |

Substituting the factorized expression (3) in Eq. (1), assuming the presence of the strong 2D optical lattice potential in the plane of , and averaging the result along coordinate , the following 2D NPSE with the nonlocal term accounting for the DD interactions is derived:

(6) | |||||

where is the sign of the contact interaction ( and correspond for the self-attraction and repulsion, respectively), and is the relative strength of the nonlocal (DD) interaction versus its local counterpart which can take positive and negative values. In particular, the nonlocal term in Eq. (6) is obtained by the integration in the direction of , assuming the long-distance limit Fisher (); Parker (), i.e., that characteristic length scales in the plane are essentially larger than the above-mentioned transverse-confinement size, .

In the presence of the strong in-plane OL potential, that will make it possible to reduce the model to the discrete form (see below), the latter condition amounts to the assumption that is essentially smaller that the lattice period. For instance, the trapping frequency kHz yields m for chromium atoms, while the OL period can be made equal to m, hence the latter condition can be readily satisfied. Actually, this means that the DD coupling may be considered as interaction between dipolar droplets with collective magnetic moments, trapped in local potential wells of the OL gpe (); Santos (), while inside the strongly confined droplets the DD interaction effectively reduces to a contact form Sanlew ()

Continuing the derivation of the 2D model, axial width , introduced above in Eq. (4), is determined, in the lowest approximation by a local relation which reduces to an algebraic equation,

(7) |

In the low-density limit, , Eq. (6) is tantamount to the GPE with the ordinary cubic local nonlinearity:

(8) | |||||

In this limit, Eq. (7) amounts to .

The choice of potential in the form corresponding to a deep OL suggests to approximate the wave function by a superposition of localized Wannier modes, wannier (), . Projecting, as usual Luca+we (), equation (6) onto the set of this modes, Eq. (6) is reduced to the set of equations

(9) | |||||

where play the role of 2D discrete coordinates replacing , the summation in the DD part excludes the term with , and is an effective lattice coupling constant, determined by the overlap integral of functions centered at two adjacent sites of the lattice. In addition to Eq. (9), discrete width function obeys the respective counterpart of algebraic equation (7),

(10) |

On the other hand, the 2D DGPE, which represents the discretization of Eq. (8), amounts to a single equation pgk2dfund (),

(11) | |||||

Stationary solutions to Eqs. (9) and (11) with chemical potential are sought for as , with discrete function satisfying, respectively, the following stationary equations:

(12) | |||||

(13) | |||||

In the former case, are related to by Eq. (10), with replaced by .

Generally, the expression for angles in the term which describes the long-range DD interaction in the above discrete equations is cumbersome. However, following Refs. Pedri05 (); Ami2 (); Ami1 (), it makes sense to focus on two most fundamental cases, with magnetic dipoles oriented either parallel or perpendicular to the plane. In the former case, we choose the orientation of the dipole moments along the . For the parallel orientation, the DD term in Eqs. (12) and (13) is anisotropic, being attractive in one in-plane direction and repulsive in the other Ami1 ():

(14) |

For the perpendicular orientation of the dipoles, the DD term is isotropic Pedri05 (); Ami2 (), taking a simple form, cf. the above expression:

(15) |

The character of the DD interaction in both cases is defined by the sign of , which is related to the sign of coefficient in the underlying GPE (1). Normally, is positive. However (as mentioned above), its sign may be reversed by means of a rapidly rotating magnetic field reversal (), allowing to take negative values Pedri05 (); Ami2 (). Thus, implies the repulsive isotropic DD interaction, while the anisotropic DD interaction is repulsive along discrete coordinate , and attractive, being twice as strong in the absolute value, along . For , the isotropic DD interaction is attractive, while in the anisotropic configuration the DD interaction is attractive along and repulsive in the direction of . The integral quantity that characterizes localized discrete modes is their norm (or power, in terms of optical models), , which is the dynamical invariant of Eqs. (9) and (11).

For both in-plane and out-of-plane (perpendicular) dipole orientations, stationary equations (12) and (13) were solved by means of a numerical algorithm based on a modified Powell minimization method, which uses a finite-difference expression for the underlying Jacobian. In the analysis reported below, we focus on the three most fundamental unstaggered localized modes (“staggered” are configurations with opposite signs of the field at adjacent sites, which is relevant to local models with the repulsive nonlinearity): on-site modes, i.e., those centered on a lattice site, inter-site ones, which are centered between lattice sites in both directions and , and hybrid modes, centered on-site along one lattice direction, and between sites along the other. This nomenclature was adopted in previous studies of local discrete models pgk2dfund (). Results reported below were obtained using a square lattice composed of x or x elements for on-site modes, and of x or x elements for inter-site and hybrid ones. It was checked that the results did not alter if an essentially larger lattices were used.

## Iii Fundamental unstaggered solitons in the case of attractive contact interactions

The nonlinear localization of matter waves occurs in gaps of the linear spectrum of the underlying system. In the case of the attractive local nonlinearity, fundamental solitons populate the semi-infinite gap wannier (); kivshargap (), which corresponds to regions and , in the GPE and NPSE models, respectively. Discrete models, derived in the tight-binding approximation, cannot correctly describe the entire linear spectrum of the underlying continuum tightb (). However, they can capture localized modes existing inside the semi-infinite gap. We will consider fundamental localized modes in this case (their topologically modified counterparts, such as discrete vortices discrete-vortices (), will be considered elsewhere).

To outline the entire existence region of the fundamental solitons, we followed the known approach, identifying regions where continuous-wave (CW) solutions with a given chemical potential, , are subject to the modulational instability, MI (it is assumed that the MI splits unstable CW states into arrays of solitons) gpe ()-NPSE (). Further details of the respective numerical procedure are reported elsewhere acta (). It is relevant to notice that the instability of uniform states in the dipolar gas is determined by two different factors: the local self-focusing, which is the usual driving mechanism for the MI at small wavenumbers of perturbations, and, in addition, the roton instability at finite wavenumbers, induced by the long-range DD interactions roton ().

### iii.1 The effect of the DD interaction on fundamental solitons

As mentioned above, in the local model with the attractive contact interactions (), CW solutions are modulationally unstable in regions and for the DNPSE and DGPE models, respectively. Note that these values correspond to the upper boundaries of the semi-infinite gaps in the corresponding continuous NPSE and GPE models. The presence of the nonlocal DD interaction extends the MI domain of the CW states, thus enabling the existence of localized modes beyond the above-mentioned limit values, and .

In the case of the attractive contact interactions, families of fundamental unstaggered bright solitons of the three above-mentioned types – on-site, inter-site, and hybrid (see Fig. 1) – have been found, as expected, precisely in the domain of the parameter space featuring the MI of the CW states. To this end, stationary equations (12) and (13) were solved for the DNPSE and DGPE models, respectively.

In the case of the attractive contact interaction, the attractive isotropic DD nonlinearity (that with the inverted sign) acts in the same direction and cannot change qualitative properties of the solitons. Therefore, only the case of the attractive on-site nonlinearity competing with its repulsive DD isotropic counterpart is really interesting, similar to the situation 1D counterpart of the present model Cuevas (), and in the Salerno model (a combined Ablowitz-Ladik/discrete nonlinear-Schödinger system) with competition between the attractive on-site and repulsive inter-site nonlinearities Zaragoza (). For this reason, in the case of the attractive contact interaction, we fix (the natural sign of the DD forces). In parallel, the interplay of the attractive on-site and anisotropic DD interactions will also be considered for .

Basic characteristics of soliton families are presented by the dependence, which shows the norm versus the chemical potential. First, in Fig. 2 we display curves generated by the DNPSE (a) and DGPE (b) models for the on-site, hybrid and inter-site discrete solitons in the presence of the attractive contact interaction, while the DD interaction is still absent. These plots may be used as the reference point for the comparison with results obtained in the presence of the DD interactions.

The dependencies for the soliton families affected by the isotropic (repulsive) and anisotropic DD interactions are displayed in Figs. 3 and 4, respectively. The comparison of the dependencies displayed in Figs. 3 and 4 with those in Fig. 2 shows that the repulsive isotropic DD interaction slightly extends the region in which the fundamental localized structures can be found, which is consistent with the results of the MI analysis for the CW solutions acta (). Simultaneously, the obtained dependencies show that there is no qualitative difference between the DGPE and DNPSE models in this case.

Usually, the Vakhitov-Kolokolov (VK) criterion, alias the slope condition, , is considered as a necessary condition for the stability of soliton VK (); isrl (). Although it is strictly applicable only to systems with a local power-law nonlinearity, it may be valid in more general situations too. Here we consider the slope criterion in both cases combining the cubic or nonpolynomial local nonlinearity with the nonlocal DD-interaction term. The slope criterion is violated at values of close to the right edge of the existence region of the fundamental solitons, see Figs. 2, 3, and 4.

The general condition which must be fulfilled for the stability of solitons is the absence of eigenvalues with positive real parts in the spectrum of small perturbations around the solitons (the spectral criterion); in fact, this means the eigenvalues must have zero real parts isrl (). This criterion was implemented numerically by computing the corresponding eigenvalues, using linearized equations for small perturbations, similar to how it was done in Ref. Luca+we (). In the case of the on-site fundamental solitons, the real part of the eigenvalues vanishes in certain regions of parameter plane (, ), as shown in Figs. 5 (a) and (c) for both discrete models. This indicates the presence of a wide stability region for the on-site modes, as reported earlier in papers dealing with 2D solitons in the discrete nonlinear Schrödinger equation, without DD interactions pgk2dfund (), aniso () and rotating (). The VK criterion is also satisfied in these cases.

In the presence of the isotropic repulsive DD interaction the region of the existence of fundamental solitons expands. These extended regions are inside the domains where CW solutions are subject to the MI, in accordance with the above statement. Two significant effects are generated by the nonlocal DD interaction in this case. First, the isotropic repulsive DD interaction makes the stability regions wider for on-site fundamental solitons with respect to the model without the nonlocal interaction, as seen in Figs 5 (b) and (d) for the DGPE and DNPSE models, respectively. Second, the nonlocal interaction generates new localized structures rising above a finite background. Properties of these new stationary structures, which are on-site centered, are considered in the following subsection. Here we only mention that the linear stability analysis shows that these solitons-on-the-background are unstable. On the other hand, in the presence of the anisotropic DD interaction, the stability domain of the on-site solitons shrinks (not shown here in detail).

In the absence of the nonlocal interactions, a known result is that the fundamental solitons of the inter-site and hybrid types, which are characterized by greater values of the norms in comparison with the on-site modes, are unstable in the whole existence region pgk2dfund (). We have checked that DD interaction of either kind (isotropic repulsive or anisotropic) do not stabilize these modes.

### iii.2 Solitons on a finite background

As briefly mentioned above, a new feature revealed by the analysis in the presence of the repulsive isotropic DD interactions (i.e., those with the natural sign) is the role of the background in the formation of localized modes. This effect is most significant near the edge of the soliton’s existence region (at close to or in the DGPE and DNPSE models, respectively), where, in the presence of the isotropic DD repulsion, narrow localized modes are found, as stationary solutions, on top of a finite uniform background, see Fig. 6. We stress that such structures cannot be found in the absence of the DD interactions. When these interactions are present, the background may play a significant role in the formation of the localized structures, because the background as a whole couples to the central peak via the long-range DD forces. Actually, the long-range coupling gives rise to these newly formed solitons on the finite background (SFBs) when its strength exceeds a certain threshold level, see Fig. 7. The threshold becomes lower with the increase of coupling constant , as seen in Fig. 7. Based on the waveform observed above the background, one can identify SFBs of the on-site, inter-site, and hybrid types.

Localized SFB structures are generated by the MI of a finite-amplitude CW states in the nonlinear lattice with the DD interaction. In the course of the development of the MI, the interplay of the local (contact) and nonlocal (DD) nonlinearities arrests the exponential growth of the instability and imposes global correlations, which enable the creation of the SFB as an outcome of the nonlinear development of the MI. However, the SFBs themselves turn out to be unstable, eventually evolving into robust breathers, which consist of a localized vibrating peak and finite oscillating background, as shown below.

### iii.3 Dynamical properties of the fundamental localized modes

In addition to the modes described above, strongly pinned on-site solitons are found in certain parametric areas in the model with the attractive contact and isotropic repulsive DD interactions. Being very narrow, these modes are not significantly affected by the DD interaction. In general, these solutions correspond to early flat parts of the plots in Figs. (2)-(4) (at and in the DGPE and DNPSE model, respectively), and they have their counterparts ion the local model. These localized structures, such as the one displayed in Fig. (8) for the DNPSE, remain ”frozen” under the action of arbitrary perturbations.

The evolution of all the unstable fundamental solitons considered in the previous subsections follows approximately the same scenario, in the presence of the DD interactions. Direct simulations of Eqs. (9) and (11) demonstrate that the outcome of the evolution under the action of perturbations is the emergence of persistent strongly pinned on-site breathers, surrounded by a finite oscillating background, see an example for the DNPSE in Fig. 9. The evolution may feature different transient stages, depending on the initial structure. For example, unstable on-site solitons directly evolve towards single-peaked tightly localized breathers, while inter-site solitons temporarily form localized structures with two peaks, see Fig. 10. The eventually formed breathers are very robust and survive perturbations of any kind.

A general conclusion is that the repulsive nonlocal DD interaction in the 2D discrete models does not prevent the formation of the extremely narrow localized structures, which are strongly pinned to the underlying lattice. These very narrow, dynamically stable modes have been associated with the quasi-collapse in the 2D lattice 2ddiscoll (), pitaevcollapse (). Nevertheless, the DD interactions produce an essential effect, as they give rise to an oscillating background which surrounds the narrow peak, being coupled to it by the long-range DD interactions, see Figs. 9 and 10; in the usual local model, with , the narrow peak would not be coupled to any background.

With the increase the strength of the DD interaction, the evolution of the corresponding solitons of the SFB type becomes sensitive to the form of perturbations. Under small asymmetric perturbations (which actually correspond to small shifts of the initial soliton’s location), or random perturbations, the SFBs decay. On the other hand, under small symmetric perturbations they evolve towards tightly bound single-peaked localized structures, surrounded by an oscillating background, similar to those displayed in Figs. 9 and 10.

The simulations demonstrate that the newly formed narrow breathers, built on top of the oscillating background, are robust modes which survive any type of perturbations. The intrinsic structure of these of the breathers depends on the initial perturbations which initiated the formation of the breather from an unstable soliton. Namely, for asymmetric or random perturbations, the repulsion exerted by the DD interactions is more significant than in the case of symmetric perturbations, making the emerging breather broader, with a richer internal structure, see Fig. 11.

Finally, we have found that the fundamental solitons obtained in the model combining attractive contact and anisotropic DD interactions with a finite strength (recall they correspond to the in-plane orientation of the dipole moments) are unstable in their entire existence region, except the on-site solitons in the region of very small (the anti-continuum limit, as concerns the local coupling). We could conclude that the anisotropy of the DD interaction makes the solitons more sensitive to small perturbations (i.e., more unstable), in comparison to the case of the isotropic DD interactions.

## Iv Fundamental solitons generated by the DD interactions in the case of the repulsive local nonlinearity

In the local model with the repulsive contact interactions (), unstaggered localized modes cannot exist. However, the presence of the DD interactions, either attractive isotropic or anisotropic, opens a possibility to create such modes, as recently shown in discrete 1D models gpe (); NPSE (), and was earlier predicted for quasi-2D continuous models in Refs. Pedri05 () (the competition of the repulsive contact and isotropic attractive DD interactions) and Ami1 () (the interplay of the repulsive contact and anisotropic DD interactions).

We have found that the long-range DD interactions between lattice sites may indeed create discrete 2D solitons in the case when the on-site nonlinearity is repulsive. Numerical calculations show that, as in the situations considered above, these discrete 2D solitons appear exactly in a part of the parameter space where the CW states are modulationally unstable (this is shown in detail in Ref. acta ()). This region belongs to the semi-infinite gap of the linear spectrum, which corresponds to and , in the GPE and NPSE models, respectively. Again, three families of fundamental unstaggered discrete solitons – of the on-site, inter-site, and hybrid types – have been found in this situation as numerical solutions of stationary equations (12) and (13) for the DNPSE and DGPE models, respectively.

It is relevant to emphasize that the solitons are also obtained in the limit of the vanishing local interaction, i.e., the discrete solitons may be supported solely by the DD interaction, which is relevant to the experiment experim1 (). In other words, the attractive nonlocal DD nonlinearity is by itself sufficient for the formation of trapped localized modes in the 2D lattices.

### iv.1 Solitons generated by the attractive isotropic DD interaction

In the case of the competition between the repulsive contact and attractive isotropic DD interactions, all unstaggered fundamental solitons have a bell-like shape, which is a consequence of the nonlocality of the interaction which creates them. This feature is also similar to what was found in the 1D counterpart of the present models gpe ().

The curves for the on-site, hybrid and inter-site unstaggered solitons in the present case are much closer to each other than it was in the case of the attractive contact interaction, see Fig. 12. On the other hand, numerical computations demonstrate a violation of the spectral stability criterion for all types of the unstaggered solitons in an almost entire existence region, except for the case of a small coupling constant, , and certain values of and . With the increase of the strength of the isotropic DD attraction, the stability window shrinks with respect to the values of , as shown in Fig. 13 for .

Direct numerical simulations show that the unstable bell-shaped localized structures evolve into the corresponding breathers, but, on contrary to the situations considered above, without significant generation of an oscillating background, as shown in Fig. 14 (almost all the initial soliton’s norm remains localized, rather than being transferred to the background, unlike the case of the attractive local interactions). The width of the localized breather remains nearly the same as that of the original unstable soliton.

### iv.2 Solitons generated by the anisotropic DD interaction

The anisotropic DD interaction can also provide for the formation of unstaggered fundamental discrete solitons of the on-site, hybrid and inter-site types in the case of the repulsive contact interactions. Examples of the characteristics for such discrete solitons at fixed and two different values of the coupling constant, (close to the anti-continuum limit) and (close to the continuum limit), are shown in Fig. 15. The respective curves for solitons of the on-site and hybrid types are very close to each other.

The linear-stability analysis indicates similar stability properties of the on-site and hybrid solitons. They are unstable in almost the whole existence region. As above, an exception is found close to the anti-continuum limit (), where stable solitons exist. The respective stability diagram is displayed in Fig. 13 (b). In this figure, the stable solitons of the on-site and hybrid types are found in the white window. The stability window for the inter-site solitons existing at is smaller than for their on-site and hybrid counterparts.

We have also considered the limit of large values of coupling constant , i.e., the near-continuum limit, for which a narrow stability region for solitons had been reported in Ref. Ami1 (), produced by the variational analysis and direct simulations of the 3D continuous GPE. Taking into account the geometry, stable quasi-2D solitary modes, characterized by long-lived, slowly decaying oscillations were observed in certain parameter regions. However, in the discrete 2D DGPE and DNPSE models considered here (without the third dimension, that was explicitly present in the analysis presented in Ref. Ami1 ()), stable anisotropic solitons were not found in this part of the existence region.

To summarize, in the model with the repulsive local interactions the long-range attractive isotropic or anisotropic DD interactions give rise to discrete solitons, but (marginally) stable ones can be found only close to the anti-continuum limit. This conclusion may be considered as a consequence of the long-range character of the DD interaction, which becomes more significant as the system is getting more discrete Parker (); Sanlew ().

## V Conclusions

In this paper we have introduced the 2D discrete model for the dipolar BEC trapped in the deep optical lattice. We have analyzed the structure and dynamics of fundamental discrete solitons in two versions of the model, which are based, respectively, on the ordinary discrete form of the GPE (Gross-Pitaevskii equation) with the cubic nonlinearity, or the 2D discrete NPSE (nonpolynomial Schrödinger equation), both including nonlocal DD (dipole-dipole) interactions. Fundamental localized modes studied in this work are counterparts of solitons existing in the semi-infinite gap, in terms of the continuum limit of the discrete models. Only the situations with competing local and nonlocal interactions (attractive/repulsive, or vice versa) appear to be interesting. We have shown close similarity in the behavior of the fundamental unstaggered solitons in the discrete models of both the GPE and NPSE types in these cases.

The presence of the repulsive isotropic DD interaction competing with the local attraction extends the existence and stability regions for fundamental on-site solitons. Due to the action of the nonlocal DD interactions, the background plays an essential role in the formation of a new type of unstaggered localized structures: near the edge of the semi-infinite gap of the corresponding continuous system, the generation of SFBs (solitons on a finite background) of the three types (on-site, inter-site and hybrid) was revealed by the systematic numerical analysis. These solitons are unstable, but they spontaneously transform themselves into robust breathers, surrounded by a finite oscillating background. Actually, the presence of the background is a feature specific to the dipolar model.

In the case of repulsive contact interactions, unstaggered fundamental solitons may exist due to the DD attraction. In the case of the isotropic attractive DD interaction, stable on-site, hybrid and inter-site solitons, with close values of the norm, can be found in the strongly discrete limit. In general, the fundamental solitons generated by the isotropic DD interaction are broad bell-shaped modes. Being unstable, they evolve into localized breathers with almost no loss of the norm.

###### Acknowledgements.

G.G., A.M., M. S. and Lj.H. acknowledge support from the Ministry of Science, Serbia (Project 141034). The work of B.A.M. was supported, in a part, by grant No. 149/2006 from the German-Israel Foundation. This author appreciates hospitality of the Vinča Institute of Nuclear Sciences (Belgrade, Serbia).## References

- (1) A. Görlitz, J. M. Vogels, A. E. Leanhardt, C. Raman, T. L. Gustavson, J. R. Abo-Shaeer, A. P. Chikkatur, S. Gupta, S. Inouye, T. Rosenband, and W. Ketterle, Phys. Rev. Lett. 87, 130402 (2001); J. Chen, J. G. Story, J. J. Tollet, and R. G. Hulet, Phys. Rev. Lett. 69, 1344 (1992).
- (2) L. Salasnich, A. Parola, and L. Reatto, Phys. Rev. A 65, 043614, (2002).
- (3) L. Salasnich and B. A. Malomed, Phys. Rev. A 79, 053620 (2009).
- (4) A. Muñoz Mateo and V. Delgado, Phys. Rev. A 77, 013617 (2008).
- (5) M. A. Porter, R. Carretero-González, P. G. Kevrekidis, and B. A. Malomed, Chaos 15, 015115 (2005).
- (6) A. Trombettoni and A. Smerzi, Phys. Rev. Lett. 86, 2353 (2001); F. Kh. Abdullaev, B. B. Baizakov, S. A. Darmanyan, V. V. Konotop, and M. Salerno, Phys. Rev. A 64, 043606 (2001); G. L. Alfimov, P. G. Kevrekidis, V. V. Konotop, and M. Salerno, Phys. Rev. E 66, 046608 (2002).
- (7) P. G. Kevrekidis, K. Ø. Rasmussen, and A. R. Bishop, Phys. Rev. E 61, 2006 (2000).
- (8) A. Maluckov, Lj. Hadžievski, B. A. Malomed, and L. Salasnich, Phys. Rev. A 78, 013616 (2008).
- (9) K. Goral, K. Rzazewski, and T. Pfau, Phys. Rev. A 61, 051601 R (2000); M. Greiner, C. Regal, and D. S. Jin, Nature 426, 537 (2003); A. Griesmaier, J. Werner, S. Hensler, J. Stuhler, and T. Pfau, Phys. Rev. Lett. 94, 160401 (2005); J. Stuhler, A. Griesmaier, T. Koch, M. Fattori, T. Pfau, S. Giovanazzi, P. Pedri, and L. Santos, ibid. 95, 150406 (2005); J. Werner, A. Griesmaier, S. Hensler, J. Stuhler, and T. Pfau, ibid. 94, 183201 (2005); A. Griesmaier, J. Stuhler, T. Koch, M. Fattori, T. Pfau, and S. Giovanazzi, ibid. 97, 250402 (2006); T. Lahaye, T. Koch, B. Fröhlich, M. Fattori, J. Metz, A. Griesmaier, S. Giovanazzi, and T. Pfau, Nature (London) 448, 672 (2007).
- (10) T. Lahaye, J. Metz, B. Froehlich, T. Koch, M. Meister, A. Griesmaier, T. Pfau, H. Saito, Y. Kawaguchi, and M. Ueda, Phys. Rev. Lett. 101, 080401 (2008).
- (11) T. Koch, T. Lahaye, J. Metz, B. Fröhlich, A. Griesmaier, T. Pfau, Nature Physics 4, 218 (2008).
- (12) M. Marinescu and L. You, Phys. Rev. Lett. 81, 4596 (1998); S. Giovanazzi, D. O’Dell, and G. Kurizki, Phys. Rev. Lett. 88, 130402 (2002); I. E. Mazets, D. H. J. O’Dell, G. Kurizki, N. Davidson, and W. P. Schleich, J. Phys. B 37, S155 (2004).
- (13) J. Sage, S. Sainis, T. Bergeman, and D. DeMille, Phys. Rev. Lett. 94, 203001 (2005); C. Ospelkaus, L. Humbert, P. Ernst, K. Sengstock, and K. Bongs, ibid. 97, 120402 (2006); T. Köhler, K. Góral, and P. S. Julienne, Rev. Mod. Phys. 78, 1311 (2006).
- (14) J. Deiglmayr, A. Grochola, M. Repp, K. Mörtlbauer, C. Glück, J. Lange, O. Dulieu, R. Wester, and M. Weidemüller, Phys. Rev. Lett. 101, 133004 (2008).
- (15) P. Pedri and L. Santos, Phys. Rev. Lett. 95, 200404 (2005); R. Nath, P. Pedri, and L. Santos, Phys. Rev. A 76, 013606 (2007).
- (16) I. Tikhonenkov, B.A. Malomed, and A. Vardi, Phys. Rev. A 78, 043614 (2008).
- (17) S. Giovanazzi, A. Görlitz, and T. Pfau, Phys. Rev. Lett. 89, 130401 (2002); A. Micheli, G. Pupillo, H. P. Büchler, and P. Zoller, Phys. Rev. A 76, 043604 (2007).
- (18) I. Tikhonenkov, B. Malomed, and A. Vardi, Phys. Rev. Lett. 100, 090406 (2008).
- (19) C. Rotschild, O. Cohen, O. Manela, and M. Segev, Phys. Rev. Lett. 95, 213904 (2005); D. Briedis, D. E. Petersen, D. Edmundson, W. Królikowski, and O. Bang, Opt. Exp. 13, 435 (2005).
- (20) J. Cuevas, B. A. Malomed, P. G. Kevrekidis, and D. J. Frantzeskakis, Phys. Rev. A 79, 053608 (2009).
- (21) G. Gligorić, A. Maluckov, Lj. Hadžievski, and B. A. Malomed, J. Phys. B 42, 145302 (2009).
- (22) B. B. Baizakov, F. Kh. Abdullaev, B. A. Malomed, and M. Salerno, J. Phys. B: At. Mol. Opt. Phys. 42, 175302 (2009).
- (23) G. Gligorić, A. Maluckov, Lj. Hadžievski, and B. A. Malomed, Phys. Rev. A 78, 063615 (2008).
- (24) M. Klawunn and L. Santos, Phys. Rev. A 80, 013611 (2009).
- (25) G. Gligorić, A. Maluckov, Lj. Hadžievski, and B. A. Malomed, Phys. Rev. A 79, 053609 (2009).
- (26) U. R. Fisher, Phys. Rev. A 73, 031602R (2006).
- (27) N. G. Parker and D. H. J. O’Dell, Phys. Rev. A 78, 041601R (2008).
- (28) L. Santos, G. V. Shlypanikov, and M. Lewenstein, Phys. Rev. Lett. 90, 250403 (2003).
- (29) B. B. Baizakov,V. V. Konotop, and M. Salerno, J. Phys. B: At. Mol. Opt. Phys. 35, 5105 (2002); E. A. Ostrovskaya, C. M. Savege, and Yu. S. Kivshar, Phys. Rev. A 67, 013602 (2003).
- (30) A. A. Sukhorukov, Yu. S. Kivshar, H. S. Eisenberg, and Y. Silberberg, IEEE J. Quant. Elect. 39, 31 (2003).
- (31) B. A. Malomed and P. G. Kevrekidis, Phys. Rev. E 64, 026601 (2001).
- (32) A. Wöllert, G. Gligorić, M. Škorić, A. Maluckov, and Lj. Hadžievski, Acta Physica Polonica A 116, 519 (2009).
- (33) L. Santos, G. V. Shlyapnikov, and M. Lewenstein, Phys. Rev. Lett. 90, 250403 (2003); S. Komineas and N. R. Cooper, Phys. Rev. A 75, 023623 (2007); R. M. Wilson, S. Ronen, J. L. Bohn, and H. Pu, Phys. Rev. Lett. 100, 245302 (2008).
- (34) J. Gómez-Gardeñes, B. A. Malomed, L. M. Floría, and A. R. Bishop, Phys. Rev. E 73, 036608 (2006); ibid. 74, 036607 (2006).
- (35) N. G. Vakhitov and A. A. Kolokolov, Izv. Vyssh. Uchebn. Zaved., Radiofiz. 16, 10120 (1973) [in Russian; English translation: Radiophys. Quantum Electron. 16, 783 (1973)].
- (36) Y. Sivan, B. Ilan, and G. Fibich, Phys. Rev. E 78, 046602 (2008).
- (37) P. G. Kevrekidis, D. J. Frantzeskakis, R. Carretero-González, B. A. Malomed, and A. R. Bishop, Phys. Rev. E 72, 046613 (2005).
- (38) J. Cuevas, B. A. Malomed, and P. G. Kevrekidis, Phys. Rev. E 76, 046608 (2007).
- (39) V. K. Mezentsev, S. L. Musher, I. V. Ryzhenkova, and S. K. Turitsyn, JETP Lett. 60, 829 (1994).
- (40) L. P. Pitaevskii, Phys. Lett. A 221 14 (1996).