Solitary Waves Under the Competition of Linear and Nonlinear Periodic Potentials
In this paper, we study the competition of linear and nonlinear lattices and its effects on the stability and dynamics of bright solitary waves. We consider both lattices in a perturbative framework, whereby the technique of Hamiltonian perturbation theory can be used to obtain information about the existence of solutions, and the same approach, as well as eigenvalue count considerations, can be used to obtained detailed conditions about their linear stability. We find that the analytical results are in very good agreement with our numerical findings and can also be used to predict features of the dynamical evolution of such solutions.
The field of Bose-Einstein condensates (BECs) in atomic physics  has provided a renewal of interest in solitary wave and nonlinear excitations. This is especially so because the leading order, mean-field description of interatomic interactions corresponds to a cubic nonlinearity. This, in turn, is responsible for the emergence and experimental observation of a variety of matter-wave solitons in BECs, including bright [2, 3], dark  and gap  solitons. Furthermore, atom optical devices have been proposed, such as, for instance, the atom chip  that would allow the possibility to controllably manipulate such nonlinear structures.
It is well-known that bright (respectively, dark) matter-wave solitons arise in BECs with attractive (respectively, repulsive) interatomic interactions, i.e., for atomic species with negative (respectively, positive) scattering length .
One of the particularly appealing features of the BEC setting is the existence of a wide variety of experimental “knobs” that can be used to manipulate or control the relevant structures. In particular, interfering laser beams can be used to produce a standing wave pattern, known as the optical lattice, providing a periodic linear potential for the condensate. This type of structure offers a large variety of interesting phenomenology including Bloch oscillations, Landau-Zener tunneling, dynamical instabilities, gap excitations among many others and a considerable amount of review works have already been dedicated to this topic [7, 8, 9].
On the other hand, magnetically-induced Feshbach resonances can be used to modify at will both the magnitude and the sign of the scattering length by tuning the external magnetic field; see e.g.  and also [2, 3] where the Feshbach resonance in Li BECs was used for the formation of bright matter-wave solitons. The ability to modulate the scattering length has led to a large variety of studies where this mechanism has been used. For instance, time-dependent modulations of the scattering length were proposed as a means of preventing collapse in higher-dimensional BECs , or as a way of producing robust matter-wave breathers , among others.
A more recent suggestion has been to add to a constant bias magnetic field a gradient in the vicinity of a Feshbach resonance, allowing for a spatial variation of the scattering length, thereby providing what has come to be termed a “collisionally inhomogeneous environment”. Notice that, given the availability of magnetic and optical (laser-) fields, the external trapping potential and the spatial variation of the scattering length can be adjusted independently (see  for more details on the relevant configuration). In this latter setting, a variety of propositions of interesting dynamical phenomena have been made concerning scenarios for the emission of solitons , delocalizing transition of matter waves , or the dynamics of the waves in random , linear , periodic [18, 19] or localized  spatial modulations. A number of more mathematically minded results on the existence and stability of waves have also appeared in  and a technique for analytically constructing exact solutions in .
The framework of collisionally inhomogeneous environments in combination with external optical lattices provides an ideal environment for competition. This was illustrated in examples of modulational instability of Bloch states  and of the delocalizing transition in one-dimension . In the present work we develop these ideas exploring more general lattice profiles. In particular, external (linear) potentials and collisional (nonlinear) potentials in (as well as out of) phase will be considered. It will be demonstrated that when in phase, these potentials provide a competition leading to a number of interesting effects including stabilization/destabilization thresholds and even the mutual annihilation of the two potentials to provide an effectively quasi-translationally-invariant environment. The effective potential landscape where the solitary waves (of the bright type) live will be obtained following the Lyapunov-Schmidt considerations of . Then, the relevant (translational) eigenvalue of the linearization will be computed based on the curvature of this effective potential landscape and the stability/instability of the waves will be assessed (and the relevant transition points will be obtained). These results will be confirmed by a second independent method based on a direct count of unstable eigenvalues. Finally, we will examine the variation of a number of relevant key parameters (such as the amplitude, the wavenumbers or the relative phase) of the potentials, in order to evaluate the validity of the approach. Furthermore, when the structures are unstable, we will examine what this approach can suggest regarding the actual instability evolution dynamics. The understanding that we will develop will enable us to manipulate the ensuing solitary waves in such a complex territory and to understand their dynamical behavior in the presence of linear and nonlinear lattices.
Our presentation will be structured as follows. In section II, we will present our analytical results; in section III, we will corroborate these results by means of numerical computations. Finally, in section IV, we will summarize our findings and present some interesting directions for future study.
2 Analytical Results: Solitary Wave Statics and Dynamics
The prototypical framework in which we will consider the above discussed competition of linear and nonlinear lattices is that of the perturbed nonlinear Schrödinger equation of the form
In Eq. (2.1), and . While we will keep the presentation of the mathematical results as general as possible, the particular case of interest in the selection of the nonlinear and linear lattice will, respectively, be:
where , , , and are real constants. Notice that the lattices have the same functional form, which will allow us to reveal more lucidly the relevant competition between the corresponding terms.
When , Eq. (2.1) has the well-known stable localized soliton solution given by
where , is the position of the soliton center, is the velocity of the soliton, and .
We presently focus on the stationary modes with . Given the monoparametric nature of the family of the respective soltuions, we can fix in what follows (in fact, we will fix in our numerical computations below). Because of the rotational and translational invariance of the unperturbed equation, this solution is unique only up to rotational and translational symmetry.
On the other hand, when , the translational invariance of the equation is broken, which may naturally lead to the potential destabilization of the localized states, depending on the perturbation parameters. This is the problem that we will examine in what follows under the influence of both linear and nonlinear lattices.
2.1 Hamiltonian Perturbation Approach
The existence and nature of localized solutions to perturbed Hamiltonian systems, of which Eq. (2.1) is a particular case, was studied in  (and subsequently in a broader setting in ). A general perturbative approach was developed in these works based on Lyapunov-Schmidt solvability conditions , and relevant stability calculations were formulated on the basis of the Evans function [26, 27]. Here, we present some of the general features of the theory, adapt our problem to the general framework of [23, 24] and subsequently apply these methods to the problem of interest. In order to apply these criteria, it is convenient to recast Eq. (2.1) as
where . Here,
Then, for fixed , the intuitive condition for the persistence of the wave is given by 
where is the previously-free parameter associated with the invariance (in the case of translation, it is the center of the pulse (2.3)). This condition implies that the wave is going to persist only if centered at the parameter-selected extrema of the energy (which are now going to form, at best, a countably infinite set of solutions, as opposed to the one-parameter infinity of solutions previously allowed by the translational invariance).
Equally importantly, from this expression and from the nature of the wave, one can infer stability information about the solution of interest. In particular, the stability of the perturbed wave is determined by the location of the eigenvalues associated with the translational invariance; previously, the relevant eigenvalue pair was located at the origin of the spectral plane of eigenvalues . On the other hand, we expect the eigenvalues associated with the invariance (i.e., the phase invariance associated with the conservation) to remain at the origin, given the preservation of the latter symmetry under the perturbations considered herein. To compute the relevant eigenvalues we refer to the framework put forth by the works of [23, 24] (adapting the notation of the latter work) in the following form. Using Proposition 6.1 of , we expect that the perturbed system eigenvalues will be given by the matrix equation:
also, is given by
In the formulation of , the relevant eigenvalues are obtained to leading order as , and denotes the solitary wave of Eq. (2.3). Hence, we conclude from Eq. (2.8) that as indicated above, the eigenvalues associated with the rotational invariance will be preserved at , while the translational eigenvalue will be shifted according to:
Hence, the corresponding eigenvalue can be straightforwardly evaluated, provided that we first compute the extrema of the effective energy landscape , which, as a function of , will hereafter be denoted as , or more precisely . This will be the effective energy landscape and the stability or instability of the configuration will be associated with the convexity or concavity of this effective energy landscape.
We now proceed to evaluate the relevant expressions of the general theory for the special case of interest herein, namely for the potentials of Eq. (2.2). The effective energy landscape can be evaluated after performing two straightforward contour integrations, that yield
In the simple case of , the above result reduces to that in , which leads to the well-known conclusion of , according to which a maximum of a linear periodic potential leads to an unstable solitary wave configuration, while the opposite is true for a minimum of a periodic potential. However, in our case, there is an intriguing interplay between the -dependent term stemming from the linear optical lattice and the -dependent term, emerging from the nonlinear optical lattice. This competition leads to the potential for stability-instability transitions for the wave, based on the properties of the trapping (such as ), but also the properties of the wave itself (since the expression of (2.12) is explicitly dependent on ).
2.2 Eigenvalue Count Approach
It is worthwhile to note, however, that the above stability results, based on the formulation of [23, 24] can also be alternatively derived using the approach of [28, 29]. We present this alternative formulation here, since we consider it to be a nice complement to the direct eigenvalue computation of [23, 24] using the Hamiltonian perturbation technique. It is well-known that the stability of the solitary wave  is determined by the number of negative eigenvalues of the operators
In particular, if denotes the count of negative eigenvalues and , then the solitary wave is unstable. The superscript in the operators is to distinguish the and the cases. In the latter, is well known to have a single negative eigenvalue () with an eigenfunction spanned by and a single zero eigenvalue with an eigenfunction spanned by , while has no negative eigenvalues and a single zero eigenvalue with an eigenvector spanned by . In the perturbed case, retains its zero eigenvalue with eigenvector spanned by . Also, from the perturbation theory of Schrödinger operators  it is known that has a negative eigenvalue near , and a second eigenvalue near , which will be denoted by . Both of these eigenvalues are analytic in , at least in a neighborhood of the real axis. Then the stability question (given also that for this branch of solutions ) is completely settled by the count of negative eigenvalues of . In particular:
If has two negative eigenvalues, then the solitary wave will be unstable since ;
If has only one negative eigenvalue, then the coherent structure will be stable.
Hence, the stability issue hinges on the shift of the zero eigenvalue (corresponding to translational invariance, when ) in the presence of the perturbation. Thus, similarly to , we will consider the quantity with being the eigenvector corresponding to the eigenvalue ; here denotes the -inner product of with . can be decomposed as , where is proportional to and . Along the lines of [31, 32] one can show that , and therefore (up to a proportionality factor) as . From perturbation theory  it is also known that the eigenvector of is analytic in . Therefore , and actually , i.e., it will be of the order of the perturbation. Then, we have
However, each of the second and third terms will be of order higher than the first (at least O, while the dominant one will be of O), hence the solitary wave stability will be determined by . But then,
is an appropriate proportionality factor (between and ). By means of a direct computation (differentiating the equation satisfied by the stationary state) one has that
which, in turn, forming the inner product with , and integrating by parts leads to the key result, namely:
The relevant integral term of the right hand side can be seen by direct inspection to be equivalent (up to a negative-definite proportionality factor) to the expression for in Eq. (2.11). Its positivity (indicating a shift of the zero eigenvalue to positive values) will imply stability, while its negativity (indicating a shift of the zero eigenvalue to negative values) will lead to instability. This conclusion is fully equivalent to the ones obtained from Eq. (2.11) [although the latter, in some sense, contains additional information yielding a quantitative measure of the relevant eigenvalue].
2.3 Solitary Wave Dynamics
In order to describe the dynamics of a soliton of Eq. (2.1) at one can employ the perturbation theory for the NLS soliton , or more precisely, the adiabatic approximation. If then and has to be found from the equations of the adiabatic approximations. The straightforward algebra yields
where is the number of particles [it is an integral of motion of (2.1)].
It follows from (2.19) that there exist different types of motion of the soliton. In particular if and the soliton dynamics reproduces the mathematical pendulum. The respective motion of a soliton can be either periodic or translational (i.e., unbounded), depending on the initial conditions. Another special case arises for , in which case the right hand side of (2.19) becomes zero and in the adiabatic approximation , i.e. the motion becomes linear because the periodic nonlinearity effectively compensates exactly the periodic potential.
More sophisticated evolution scenarios can be observed for depending whether they are commensurable or not.
3 Numerical Results
We now proceed to describe our numerical results comparing with the analytical prediction of the previous section. We use (for which we expect the perturbative description to still be meaningful), and vary the relative parameters of the two lattices (linear and nonlinear).
Our first set of numerical results consists of setting and , and varying . In this way, we can test the validity of our predictions for amplitude variations (in this case of the nonlinear lattice). Our results are summarized in Fig. 1. We have varied , finding that there is a stability change within this interval. In particular, the left panel of the figure shows the case of which is unstable, and of the stable . The right panel shows this transition in terms of the real part (and also of the square) of the relevant eigenvalue associated with the translational mode. It is found that this eigenvalue pair starts out as real, for small , and becomes imaginary for . We use Eq. (2.11) to theoretically predict this transition as occurring at . While we see that both qualitatively and fairly quantitatively the dependence of the eigenvalue on the parameter is captured accurately by our theoretical result, it is meaningful to rationalize the error in the critical point estimation. It is, in fact, observed that the soliton does not maintain its amplitude in this continuation process (as a function of ), but rather that its amplitude is reduced from for to for . This clearly shifts the critical point upwards, whose analytical expression in this setting of and can be easily seen to be . This is in agreement (in fact, even quantitatively, if one uses the above amplitude variation) with what is observed in our numerical results.
In the bottom panel of the figure, we show the result of the unstable dynamical evolution for the case of . It can be seen that as a result of the dynamical instability the solitary wave starts moving to the left, eventually executing oscillations between the two maxima of the effective potential of Eq. (2.12). In the same plot, we show the result of the adiabatic soliton perturbation theory in this case (this is a rather “stringent” test of the theory given the unstable dynamical evolution). We observe that the Eq. (2.19) performs well in approximating the soliton trajectory over the first oscillatory cycle. However, for longer times, we observe it to gradually increasingly fail to capture the relevant oscillation. This can be seen to be due to the fact that the solitonic trajectory emits small wakes of radiation as it arrives at the turning points, resulting in a weakly damped oscillation, a feature which is not captured by our present considerations. However, we note in passing that methods similar to those developed by Soffer and Weinstein  can be used to rigorously account for such corrections.
Our second parameter variation involved the role of the wavenumbers. In particular, for the results reported in Fig. 2, we have used , and fixed and , varying . One can see that in this case, the effective potential landscape changes significantly in term of its local structure (in the previous example, it did not change, in that it was simply two cosinusoidal terms with different signs, so it was simply a matter of which had the largest “strength”). The other important feature is that in this case, as well, there is a transition from instability to stability, as is increased. In fact, the theoretical prediction for the critical point is , while the numerical one is . Once again this can be seen for the two different settings of the left panel (the unstable case of and the stable case of ), and is captured extremely accurately by the prediction of Eq. (2.11) about the location of the relevant eigenvalue (associated with translation). Furthermore, once again, the effective potential landscape that can be computed from Eq. (2.12) can provide very useful information not only about the stability of local extrema but also about the instability dynamics. The latter is observed in the bottom panel of the figure. Note, however, that while the effective potential predicts accurately the turning points of the solitary wave dynamics, the situation is more complicated with the dynamical equation of motion of (2.19). While, once again, the latter predicts very accurately the first oscillation cycle, its non-accounting of the radiative corrections of the motion leads to dynamics that overcomes the shallow potential barrier at ; this is not true, however, for the full PDE dynamics. This should serve as a note of caution in regard to using the adiabatic approximation in such (marginal) cases.
We also explored the role of the phase difference between the linear and nonlinear lattice, by varying , for and . One of the particularly interesting features of this example is that while the instability of the original configuration is not modified by this variation, the location of the solitary wave is. This is naturally expected on the basis of Eq. (2.12). In particular, we observe that the bright soliton’s center location features an oscillation around , of period (as expected); for , the wave is shifted to the right, while for , it lies to the left of the origin. In this case, we examine both the prediction for the relevant unstable eigenvalue, as well as the prediction of our theoretical results for the location of the center of the structure. The numerical results once again align extremely well with the theoretical ones, confirming the validity of our theoretical findings.
Finally, we also considered a case in the neighborhood of the complete mutual cancellation of the two contributions of the effective potential. In particular, for the case of , and for , , we examined the dynamics for , whereby in accordance with Fig. 1, the relevant (translational) eigenvalue is numerically found to be returning to , thereby restoring a regime of effective translational invariance. In this setting, we initialized a stationary soliton boosted by a factor of , with , in Fig. 4 (other values of were also used with similar results). We observe that the soliton appears to propagate with a speed near the originally “assigned” one, being submitted only to very weak modulations due to the very weak (in this case) effective potential. These modulations are accompanied by oscillations of the solitary wave amplitude and lead to a speed slightly larger than 0.1 (shown by thick dashed line in Fig. 4). We note that in this way we can induce the robust motion of the waves over the combined linear and nonlinear lattice terrain.
In this paper, we have examined the evolution of bright solitary waves in the presence (and competing effects) of linear and nonlinear lattices. We have computed the effective potential landscape that the wave encounters and have explained how its curvature is associated with the wave stability. The relevant translational eigenvalue has been explicitly computed and the transitions from stability to instability due to the zero crossing of this eigenvalue have been quantified as a function of the system’s parameters. The same threshold condition has been obtained independently from a direct eigenvalue calculation. It has been shown that these theoretical frameworks capture accurately the location of the stationary waves, as well as the pertinent eigenvalues, hence they constitute valuable tools for inferring existence and stability information about the coherent structures of such models. On the other hand, the ensuing potential energy landscape can be used to derive a dynamical equation for the motion of the soliton; however, there is a number of notes of caution that should be made in that regard, as the latter may not capture entirely accurately the dynamical behavior, especially in “marginal” cases, due to the role of radiative corrections. It is also possible to appropriately tune the system parameters so as to nearly mutually cancel the effects of the linear and nonlinear lattice and produce a wave that is propagating at nearly-constant speed.
It would be especially interesting to apply similar considerations to the case of dark solitons in repulsive BECs and examine their impact on the spectrum in the spirit of the recent work of . Similar considerations could subsequently be extended to higher dimensional settings, where an important relevant example would be the influence of the lattice potentials to the existence and stability of structures with vorticity [37, 38].
PGK acknowledges support from NSF-CAREER, NSF-DMS-0505663 and NSF-DMS-0619492, as well as the warm hospitality of MSRI during the final stages of this work. CKRTJ acknowledges the support of NSF-DMS-0410267 and the warm hospitality of MSRI. VVK acknowledges support from Ministerio de Educación y Ciencia (MEC, Spain) under the grant SAB2005-0195 and support of FCT and FEDER under the grant POCI/FIS/56237/2004.
- C.J. Pethick and H. Smith, Bose-Einstein condensation in dilute gases, Cambridge University Press (Cambridge, 2002); L.P. Pitaevskii and S. Stringari, Bose-Einstein Condensation, Oxford University Press (Oxford, 2003).
- K. E. Strecker et al., Nature 417, 150 (2002).
- L. Khaykovich et al., Science 296, 1290 (2002).
- S. Burger et al., Phys. Rev. Lett. 83, 5198 (1999); J. Denschlag et al., Science 287, 97 (2000); B. P. Anderson et al., Phys. Rev. Lett. 86, 2926 (2001); Z. Dutton, et al., Science 293, 663 (2001).
- B. Eiermann et al., Phys. Rev. Lett. 92, 230401 (2004).
- R. Folman et al., Adv. Atom. Mol. Opt. Phys. 48, 263 (2002); J. Reichel, Appl. 469 (2002); J. Fortagh and C. Zimmermann, Science 307 860 (2005).
- P.G. Kevrekidis and D.J. Frantzeskakis, Mod. Phys. Lett. B 18, 173 (2004).
- 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).
- S. Inouye et al., Nature 392, 151 (1998); J. Stenger et al., C. E. Wieman , Phys. Rev. Lett. 81, 5109 (1998); S. L. Cornish et al.,
- F. Kh. Abdullaev et al., Phys. Rev. A 67, 013605 (2003); H. Saito and M. Ueda, Phys. Rev. Lett. 90, 040403 (2003); G. D. Montesinos et al, Physica D 191 193 (2004).
- P. G. Kevrekidis et al., Phys. Rev. Lett. 90, 230401 (2003); D. E. Pelinovsky et al., Phys. Rev. Lett. 91, 240201 (2003); D. E. Pelinovsky et al., Phys. Rev. E 70, 047604 (2004); Z. X. Liang et al., Phys. Rev. Lett. 94, 050402 (2005); M. Matuszewski et al., Phys. Rev. Lett. 95, 050403 (2005).
- G. Theocharis et al., Phys. Rev. A 74, 053614 (2006).
- M.I. Rodas-Verde, H. Michinel, and V.M. Pérez-García, Phys. Rev. Lett. 95, 153903 (2005).
- Yu.V. Bludov, V.A. Brazhhyi, and V.V. Konotop Phys. Rev. A (to be published); cond-mat 0706.0079.
- F.Kh. Abdullaev and J. Garnier, Phys. Rev. A 72, 061605 (2005); J. Garnier and F.Kh. Abdullaev, Phys. Rev. A 74, 013604 (2006).
- G. Theocharis et al., Phys. Rev. A 72, 033614 (2005).
- H. Sakaguchi and B.A. Malomed, Phys. Rev. E 72, 046610 (2005); M.A. Porter et al., Physica D 229, 104 (2007).
- Y. Bludov and V.V. Konotop, Phys. Rev. A 74, 043616 (2006).
- M.T. Primatarowa, K.T. Stoychev and R.S. Kamburova, Phys. Rev. E 72, 036608 (2005).
- G. Fibich, Y. Sivan and M.I. Weinstein, Physica D 217, 31 (2006); Y. Sivan, G. Fibich and M.I. Weinstein, Phys. Rev. Lett. 97, 193902 (2006); P. Torres, Nonlinearity 19, 2103 (2006).
- J. Belmonte-Beitia et al., nlin.PS/0611051.
- T. Kapitula, Physica D 156, 186-200, 2001.
- T. Kapitula, P.G. Kevrekidis and B. Sandstede, Physica D 195, 263-282 (2004).
- M. Golubitsky and D.G. Schaeffer, Singularities and Groups in Bifurcation Theory, vol. 1, (Springer-Verlag, New York, 1985).
- J. Alexander, R. Gardner and C.K.R.T. Jones, J. Reine Agnew. Math. 410, 167-212 (1990).
- T. Kapitula, SIAM J. Math. Anal. 30, 273-297 (1999).
- Y.-G. Oh, Comm. Math. Phys. 121, 11-33 (1989).
- C.K.R.T. Jones, Ergodic Theory and Dynamical Systems 8, 119 (1988).
- T. Kato, Perturbation Theory for Linear Operators. Springer-Verlag, New York (1966).
- A. Floer and A. Weinstein, J. Func. An. 69, 397-408 (1986).
- Y.-G. Oh, Comm. Part. Diff. Eqns. 13, 1499-1519 (1988).
- F. Riesz and B. Sz.-Nagy, Functional Analysis, Dover, New York (1990).
- V. I. Karpman and E. M. Maslov, Zh. Eksp. Teor. Fiz. 73, 537 (1977) [Sov. Phys. JETP 46, 281 (1977)]
- See e.g., A. Soffer and M.I. Weinstein, Inventiones Mathematicae 136, 9 (1999); J. Stat. Phys. 93, 524 (1998).
- D.E. Pelinovsky and P.G. Kevrekidis, cond-mat/0610881.
- P.G. Kevrekidis et al., J. Phys. B: At. Mol. Opt. Phys. 36, 3467 (2003).
- P.G. Kevrekidis et al., Mod. Phys. Lett. B 18, 1481 (2004).