Modulational instability and localized breather in discrete Schrödinger equation with helicoidal hopping and a power-law nonlinearity

JID:PLA AID:25076 /DIS Doctopic: Nonlinear science

[m5G; v1.235; Prn:25/04/2018; 14:56] P.1 (1-9)

Physics Letters A ••• (••••) •••–•••

Contents lists available at ScienceDirect

Physics Letters A www.elsevier.com/locate/pla

Discussion

Modulational instability and localized breather in discrete Schrödinger equation with helicoidal hopping and a power-law nonlinearity Wei Qi a,∗ , Zi-Hao Li a , Jiang-Tao Cai a , Guan-Qiang Li a , Hai-Feng Li b a b

School of Arts and Sciences, Shaanxi University of Science and Technology, Xi’an 710021, China School of Science, Xi’an Technological University, Xi’an 710032, China

a r t i c l e

i n f o

Article history: Received 23 November 2017 Received in revised form 15 March 2018 Accepted 17 April 2018 Available online xxxx Communicated by A.P. Fordy Keywords: Modulational instability Localized breather DNLS equation Power-law nonlinearity Long-range hopping

a b s t r a c t We investigate the discrete nonlinear Schrödinger model with helicoidal hopping and a power-law nonlinearity, motivated by the tunable nonlinearity in the model of DNA chain and ultra-cold atoms trapped in a helix-shaped optical trap. In the study of modulational instability, we find a successive destabilization along with increasing nonlinear-power. In particular, the critical amplitudes of secondstage instability decrease as nonlinear-power increases. Furthermore, it is shown that information on the stability properties of weakly localized solutions can be inferred from the plane-wave modulational instability results. This link enable us to analytically estimate the critical parameters at which the breather solutions turn unstable, and find these parameters are dramatically influenced by the nonlinearpower. The stability properties of localized breathers perform an obvious change when the nonlinear power crosses a critical value γcr . It is reflected that at weak nonlinearity the breathers exhibit monostability, while exceeding γcr the bistability and instability will set in. The interplay between nonlinear-power and long-range hopping on the stability properties of breathers is also discussed in detail. © 2018 Elsevier B.V. All rights reserved.

1. Introduction As one of the most fundamental nonlinear lattice dynamical models, the discrete nonlinear Schrödinger (DNLS) [1] equation originally appeared in the connection of polarons in molecular crystals [2], but finding immediate modern applications in fiber optics and quantum condensates trapped in the periodic optical lattices. In atomic physics, the DNLS equation served as a model for Bose–Einstein condensates (BECs) [3] or superfluid Fermi gases [4] trapped in a deep optical lattice. The researches on this field have experienced a huge growth in the past decades, including Bloch oscillations [5] and discrete breather [6]. In biology physics field, the DNLS equation can be used to study the energy and charge transport or storage in DNA chain [7], and also used to describe protein in the collapsed phase [8]. The researches on DNLS equation are not only restricted to nearest neighbors (NN) dispersive hopping, but also long-range hopping coming from geometric deformations of lattice are of interest. A specific proposal for this deformation is by introducing a helicoidal structure, by which lattice sites which are far apart in the one-dimensional model can be close enough in the

*

Corresponding author. E-mail address: [email protected] (W. Qi).

https://doi.org/10.1016/j.physleta.2018.04.038 0375-9601/© 2018 Elsevier B.V. All rights reserved.

three-dimensional case. This idea can be extended to model ultracold neutral atoms in nanofiber-based double-helix dipole trap [9]. In molecular biology, this hopping can also be found in helicoidal Peyrard–Bishop (PB) model, or called Peyrard–Bishop– Dauxois (PBD) model which takes the double-helix structure of DNA into consideration [10]. Many researches have been made in this field, including Bloch dynamics in lattices with long-range hopping [11–13] and self-trapping [14], which indicate the significant impact of helicoidal hopping terms. Besides that, the nonlinearity also plays an important role in the DNLS equation, many new physics have induced by the nonlinear term. Besides the commonly studied cubic nonlinearity, sometimes we encounter a more complex on-site power-law nonlinearity of the form ||2γ , where γ is the nonlinearity power. On the one hand, in ultracold atomic fields, interatomic interactions lead to many meaningful nonlinear effects: such as γ = 1, the model describes the s-wave scattering between bosons, while γ = 2, the model describes Tonk–Girardeau gas in optical lattices [15]. When 2/3 < γ < 1, the model describes superfluid Fermi gas along the BEC to BCS crossover [16,17]. On the other hand, in the study of DNA chain, the power-law nonlinearity is used to describe the strength of hydrogen bond between base-pairs [7]. In the DNLS equation, the fundamental states are discrete breathers. In the nearest neighbor DNLS model with cubic [1] or

JID:PLA

AID:25076 /DIS Doctopic: Nonlinear science

[m5G; v1.235; Prn:25/04/2018; 14:56] P.2 (1-9)

W. Qi et al. / Physics Letters A ••• (••••) •••–•••

2

power-law nonlinearity [16,18], the breathers have been discussed in detail. Then, the DNLS equation with only cubic nonlinearity is extended by an N-th neighbor hopping term to describe the dynamics of the helicoidal shaped lattices, this model was also successfully used to describe the collective dynamics of ultracold atoms trapped in a helix-shaped optical traps [19,20]. Up to now, to the best of our knowledge, combing the helicoidal hopping and power-law nonlinearity in DNLS equation to study the modulational instability (MI) and breather is still absent. Therefore, in this work, we study the influence of nonlinear power γ on the DNLS model with the presence of long-range hopping in one dimensional case, from the aspects of both the MI and the localized breather solutions. In our mainly concerned situation, we find that the interplay of long-range hopping and nonlinearity dramatically change the condition of MI and the stability properties of breather. We find that there is a critical value of the2 nonlinearity power, for γ < γcr , each value of norm P = j |ψ j | possesses exactly one stable soliton, but when γ > γcr , bistability and instability arise in the branch of breathers, these results are conformed with analytical analysis, and showing that the stability properties of breathers are dramatically modified by the nonlinear power. We also find that the second-stage MI is linked to the stability properties of localized breathers [22] in this DNLS equation with general power-law nonlinearity. This link enables us to obtain the approximation analytical expressions of breather critical parameter, which reveals how the power-law nonlinearity and long-range hopping term change the breather stability significantly. The paper is structured as follows. Section 2 provides a brief introduction of the discrete nonlinear Schrödinger model with helicoidal hopping and power-law nonlinearity and discusses its application in molecular biology and ultracold atomic physics. Section 3 demonstrates the dispersion relation and study the MI features through linear stability analysis. In Section 4, we discuss the breather stability properties and their relationship with the MI, then investigate the conditions for the occurrence of bistability. The conclusions are presented in Section 5. The details derivations of variational and continuum approximations are given in Appendix A. 2. Discrete nonlinear Schrödinger model with helicoidal hopping and power-law nonlinearity We consider a generalized DNLS model with long-range hopping and arbitrary onsite power-law nonlinearity is

i

d j dτ

= −t 1 ( j +1 +  j −1 ) − t N ( j + N +  j − N ) − | j |2γ  j , (1)

where  j is the complex wave function at site j ∈ Z, τ denotes time, t 1 and t N are positive hopping terms for the nearest and fixed N-th neighbor, respectively. γ > 0 is the parameter to tune the nonlinearity. In the specific case γ = 1, the power-law nonlinearity reduces to cubic nonlinearity. When N-th neighbor hopping t N = 0, the equation represents the model with only NN hopping. By rescaling time the right side of DNLS equation can be multiplied by a constant, thus one can set t N = 1 without loss of generality. Although our analytical results are for most general cases, that is, for arbitrary N and t N . In order to investigate the role of power-law nonlinearity γ combing with long-range hopping term t N playing in this helicoidal lattices, in our numerical simulation we mainly focus on the cases of N = 2 as example, it is the so-called the zigzag lattices as refereed in Ref. [21]. One wants to get the detailed analysis of how the helicoidal hopping term t N and general N to influence the lattice dynamics, one can refer to Ref. [22], these situations are beyond our concerned scope.

Fig. 1. (a) Helicoidal lattice including hopping coupling, black solid lines denote NN hopping, red dashed lines denote the N-th neighbor hopping (N = 4 here); (b) PBD model of DNA chain, black and brown balls denote two strands of nucleotides, black and brown dashed lines denote coupling of NN nucleotides, red dashed lines denote N-th neighbor hopping of nucleotides at different strands, and the double-blue lines represent the hydrogen bonds between base-pairs. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

As shown in Fig. 1(a), where the N-th neighbor hopping term comes from helicoidal structure of lattice, with lattice sites represented by black ball placed equidistantly along the helicoidal curve. The black solid lines denote the interactions between the adjacent lattice sites. One of the consequences of introducing helicoidal structure is that lattice sites which are far apart in the one dimensional model can be close enough in three dimensional structure, which can form the interaction between adjacent windings of the helix, indicated by the red dashed lines. Here, we first discuss its application in molecular physics. The DNLS equation may serve as a PBD model for DNA chain which presents as a double helix structure. In this model, the DNA chain is spirally coiled by two nucleotide strands as shown in the Fig. 1(b). The NN coupling between adjacent nucleotides represented by the black or brown dashed lines and the nonlinear interaction of hydrogen bond represented by the blue dashed line, where the helicoidal hopping terms are not generated from the interactions between the nucleotides on the same strand, but from the interactions between the j-th nucleotide of one strand and ( j ± N )th nucleotides of the other strand, represented by the red dashed line. It is assumed that this additional term is harmonic, just like the NN interactions. Usually the nonlinear potential that model the hydrogen bonds between base-pair is either a Morse potential or a Lennard–Jones potential [23], but as these potentials are all very complicated to obtain any analytical results, it is suggested to use the power-law nonlinear potential as a substitute. The nonlinear power is included to have the possibility to tune the coupling strength of hydrogen bonds [7,24]. Therefore, the PBD model is formulated as Eq. (2), where t 1 is the NN hopping term to describe the coupling of adjacent nucleotides, t N is the N-th neighbor hopping term to describe the helicoidal hopping introduced by helicoidal structure, and γ is the power to describe the nonlinearity of hydrogen bond. Meanwhile, in the ultra-cold atomic fields, for deeper understanding of the nonlinear quantum dynamics of ultracold atomic gases trapped in helix-shaped optical lattices, within mean-field theory, the interaction between ultracold atoms can be expressed as a power-law nonlinearity as discussed in Section 1, hence the above mentioned DNLS is also suitable to describe the nonlinear dynamic properties of this quantum system. 3. Modulational instability MI is the fundamental topic in nonlinear wave equations, which is depicted a specific range of wave numbers of plane-wave profiles becomes unstable to modulations [25]. In this section, we dis-

JID:PLA AID:25076 /DIS Doctopic: Nonlinear science

[m5G; v1.235; Prn:25/04/2018; 14:56] P.3 (1-9)

W. Qi et al. / Physics Letters A ••• (••••) •••–•••

Fig. 2. Dispersion relation in the absence of nonlinearity for different N, t 1 = t N = 1.

3

Fig. 3. F as a function of unperturbed wave number k, perturbation wave number q, black solid line denotes F = 0. t 1 = t N = 1, N = 2.

cuss the discrete modulational instability of our model under the plane-wave perturbation. First, by factorizing  j = ψ j exp(−i μτ ), we get the stationary counterpart of Eq. (1)

μψ j = −t 1 (ψ j+1 + ψ j−1 ) − t N (ψ j+N + ψ j−N ) − |ψ j |2γ ψ j . (2) The stationary wave function ψ j has plane-wave solutions of the form ψ j = A exp(ikj ) with complex amplitude A and quasimomentum k which can be restricted to first Brillouin zone −π < k < π . Inserting this solution into Eq. (2) we get the dispersion relation

μ(k) = μ0 (k) − | A |2γ ,

(3)

where

μ0 (k) = −2t 1 cos(k) − 2t N cos(kN ).

(4)

According to the definition of phase velocity v p = μ/k and group velocity v g = ∂ μ/∂ k, nonlinearity makes no difference to group velocity, but lowers the phase velocity. μ0 (k) is the dispersion relation in the absence of nonlinearity, as shown in Fig. 2, where we only give the curve for k varying from 0 to π as μ0 (k) is an even function with period 2π . Compared with the NN model where μ0 (k) is a monotonous increase function of k, the helicoidal hopping term actually induces new extrema, and the number of extrema also increases along with the growing N. The stability properties of stationary solution ψ j = A exp(ikj ) can be different due to different amplitude A and quasi-momentum k, this can be identified by means of linear stability analysis. The perturbation introduced into plane-wave solution gives the ansatz

  ∗  j = A +  ae i (qj −ωτ ) +  b∗ e −i (qj −ω τ ) e i [kj −μ(k)τ ] ,

(5)

where q and ω are the perturbation wave number and perturbation frequency, respectively. Since  is small, one can insert this ansatz into Eq. (1) and linearize in  . This yields an eigenvalue problem for the perturbation frequency ω with a pair of solutions

ω± = 2





tl sin(ql) sin(kl) ± 2 F ( F − γ | A |2γ ).

(6)

l =1 , N



ql

Here, F = 2 l tl cos(kl) sin2 ( 2 ). If we set the nonlinearity power γ = 1, Eq. (6) reduces to cubic nonlinear case in Ref. [22]. Real part in the ω± corresponds to the vibration around plane wave solution, whereas imaginary part in the ω± corresponds to the exponential growth of perturbation and thus corresponding to the MI. For a plane wave of amplitude A and quasi-momentum k towards a perturbation of wave number q, the imaginary part of corresponding ω± arises from the square root when F ( F − γ | A |2γ ) < 0.

Fig. 4. Absolute value of imaginary part of ω± as a function of unperturbed wave number k, perturbation wave number q, with different amplitude square | A |2 and nonlinearity γ . Throughout, t 1 = t N = 1, N = 2.

From this analytical expression, we can see clearly that the nonlinearity power γ can dramatically influence the MI condition. In the following, we first consider how the nonlinearity power influence the region where instability might be present. As shown in Fig. 3 where F as a function of k and q with t 1 = t N = 1, N = 2, two kinds of regions can be recognized through the black solid lines denoting F = 0, one represents F < 0 and the other represents F > 0. Noticing that γ > 0, and | A |2 > 0, the instability can only be occurred on the region of F > 0, as clearly see F ( F − γ | A |2γ ) < 0 can only be satisfied when γ | A |2γ is large enough. The other kind of region, however, is impossible to be unstable because F and F − γ | A |2γ are of the same sign for the negative F . The expression of F does not include γ , which means the boundary of these two kinds of regions is not affected by γ . Although γ is irrelevant to the boundary mentioned above, its influence on |Im(ω± )| in “unstable” region is obvious, as shown in Fig. 4 with t 1 = t N = 1, N = 2. Along with the increasing A or γ , the instability is initially appeared near the boundary of F = 0 where the required γ or | A |2 is relatively small here, then other places successively turn unstable until the value of γ | A |2γ reaches the max( F ) = 3.125. We first discuss the instability of k = 0 plane wave, it is noticeable that both nonlinear-power and amplitude can induce two-stage instability as shown in Fig. 4, and MI subse-

JID:PLA

AID:25076 /DIS Doctopic: Nonlinear science

[m5G; v1.235; Prn:25/04/2018; 14:56] P.4 (1-9)

W. Qi et al. / Physics Letters A ••• (••••) •••–•••

4

Fig. 5. Critical amplitude square | A cr |2 of second-stage MI as a function of t 1 = t N = 1, N = 2.

γ , k = 0,

quently occur in two disconnected q-intervals with increasing γ . The first stage of instability occurs in the vicinity of q = 0 at infinitesimal amplitudes or nonlinear-power, the second-stage of instability occurs in the vicinity of q = π with a critical value of γ or amplitude square | A |2 . For a fixed value of γ , following the increasing of the | A |2 , the two-stage MI as will be appeared as shown in the same panel of Fig. 4. However, for a fixed | A |2 , as increasing of γ , the two-stage instability will also be occurred, as depicted in the same column of Fig. 4. Secondly, we find that although the second-stage instability always occurs around q = π , regardless the value of γ , the critical value of amplitude becomes smaller and smaller with the increasing γ . As shown in the Fig. 4 that the critical values are | A cr |2 = 3, 2, 1 for γ = 0.816, γ = 1 and γ = 2, respectively. To understand this one can set k = 0, the expression of F becomes 2

F (k = 0, q) = 2t 1 sin

q 2 1



+ 2t N sin

2

qN



2

(7)

= t 1 + t N + μ0 (q) ≥ 0. 2

The MI of the k = 0 wave requires F − γ | A |2γ < 0, thus

 | A |2 >

2t 1 + 2t N + μ0 (q)

1 γ

2γ

(8)

,

for N = 2, the global minimum of μ0 (q) in the interval of 0 < q < π is at q = 0 as shown in Fig. 2, thus the first stage instability will emerge from there at small amplitudes. Beyond this, the additional local minimum of μ0 (q) is at q = π , thus the second-stage instability will emerge at q = π , with the critical value of amplitude given by

 | A cr |2 =

2t 1 + 2t N + μ0 (q = π ) 2γ

1 γ

=

2t 1

γ

1 γ

.

(9)

Eq. (9) that the critical amplitude is a monotonic decreasing function of γ and the speed of reduction becomes smaller and smaller, as clearly depicted in Fig. 5. For general N, as long as the t N is not very small, the minima of μ0 (q) are mainly determined by cos(kN ). Therefore, the minima appear around qm = m 2Nπ (m = 0, 1, ...,  N2 ) with the corresponding critical amplitudes

 | A cr |2(m)

=

2t 1

γ

2

sin

 π  γ1 m

N

.

(10)

An example of dynamics with γ = 0.816, | A |2 = 4 is shown in Fig. 6(b), seeding a modulationally unstable k = 0 plane wave with

Fig. 6. (a) The value of |Imω± | as a function of k and q; (b) decay dynamics of k = 0 plane wave ψ j with noise; (c) decay dynamics of ϕ (k), the discrete Fourier transformation of ψ j . Throughout, | A |2 = 4, γ = 0.816, t 1 = t N = 1, N = 2.

arbitrary small fluctuation. In the beginning the solution keeps the originate plane-wave form, but quickly localizes after a specific time. In contrast, in Fig. 6(c) we give ϕ (k), the discrete Fourier transform of ψ j . It is seen that at start the ϕ (k) localizes in k = 0, later grows on the most unstable q-values as indicated in corresponding Fig. 6(a), at the time when plane wave localizes in direct space, the ϕ (k) quickly extends to the entire k space. 4. localized breather solutions and its stability properties Time-periodic spatially localized excitations, called discrete localized breathers, are typical solutions for the DNLS equation. In this section we discuss the properties of on-site breather solutions, which are characterized by predominant occupation of a single lattice site. We start from the localized breather solutions of Eq. (2). A key quantity in analyzing the stability properties of breather branch is  2 its norm P = j |ψ j | . First we use numerical simulation to obtain the P (|μ|) curves for different values of γ , as shown in Fig. 7(a)(b).

JID:PLA AID:25076 /DIS Doctopic: Nonlinear science

[m5G; v1.235; Prn:25/04/2018; 14:56] P.5 (1-9)

W. Qi et al. / Physics Letters A ••• (••••) •••–•••

Fig. 7. (a) P (|μ|) curve at different Throughout, t 1 = t N = 1, N = 2.

γ for γ < 1; (b) P (|μ|) curve at different γ for γ > 1; (c) P (a) as a function of a for γ = 0.7; (d) P (a) as a function of a for γ = 1.1.

No matter what nonlinear power is, the breather starts to exist at |μ| = |μ0 | = | − 2(t 1 + t N )| = 4 with zero norm. In the interval of |μ| > 4, solutions exhibit different kinds of stability with a critical value of γ . For example, when γ = 0.6, the norm P is a monotonic increasing function of |μ|, one can value of P corresponds to only one stable solution and there is no bistability. However, when γ ≥ γcr = 0.64, P (|μ|) curve exhibits a local maximum ( P max ) before a local minimum ( P min ), such that while in the interval of P > P max or P < P min one value of P corresponds to exact one stable solution, in the interval of P min < P < P max three different solutions coexist. The two solitons with the smallest and greatest frequency are stable, which indicates that bistability sets in. Accompanied by this bistability, the instability also occur. The third P (intermediate) solution at the interval of negative derivative dd|μ |, according to the Vakhitov–Kolokolov-type [26] stability criterion, is unstable. The frequency (μcr ) at the P max is the critical frequency for the occurrence of instability. These stability features use a theoretical basis in P (a) curve obtained by variational approximation, where a is a real parameter in exponential ansatz ψ j = B exp(−a| j |). We put the detailed derivation of P (a) expression in Appendix A, the result is

Pγ =

2 sech(a) tanhγ +3 (a) coth[(γ + 1)a] csch2 a − coth a csch2 [(γ + 1)a]

× [t 1 + Nt N e

(−aN )

5

(11)

( N + tanh a) cosh(a)].

The P (a) curves are displayed in Fig. 7(c)(d), from which we also see different types of stability for different value of γ . At γ = 0.7, P is a monotonic increase function of a, and all soliton solutions are monostable. The local minimum and local maximum of P (a) curve appear when γ crossing a critical value γcr = 0.85, accompanied by the bistability and instability properties. These extrema and properties can be clearly identified at γ = 1.1. These analyses are in accordance with the results in numerical simulation, the only difference is that in the variational approximation the critical value of γ is 0.85, whereas in the numerical simulation the critical value is 0.64. This is due to the distinction of variational ansatz and breather solution.

Fig. 8. Critical value of nonlinearity γcr as a function of hopping term t N , solid line denotes variational estimate obtained by Eq. (11), markers denote numerical results, t 1 = 1, N = 2.

Numerical simulation demonstrates that when t N = 1, bistability sets in when γ exceeds a critical value of 0.64. Evidently, this kind of critical nonlinear-power also exists in the NN model where t N = 0, but the value is 1.34 [22,18]. The difference implies that critical nonlinear-power changes at different value of t N . This interplay of γ and t N shows that the bistability property is determined by the interplay of nonlinear power and helicoidal hopping. We study this relationship with numerical simulation and variational approximation method, the results are clearly shown in Fig. 8, its reflects that the long-range hopping term dramatically changes the critical nonlinear-power. For investigating the specific breather profile and its stability properties, in Fig. 7(a) we choose three solitons of the same norm P = 10.2 from the γ = 0.8 branch of breathers for example, indicated by the different type of markers. The corresponding profiles are shown in the Fig. 9(a) using the same markers. Since the norm of breathers is a conserved quantity, the unstable solution in negative slope of P (|μ|) curve will rearrange to its coexisting two stable solutions of the same norm. A example of dynamic is shown in Fig. 9(b) where we seed the μ = −5.34 unstable solution with

JID:PLA

6

AID:25076 /DIS Doctopic: Nonlinear science

[m5G; v1.235; Prn:25/04/2018; 14:56] P.6 (1-9)

W. Qi et al. / Physics Letters A ••• (••••) •••–•••

2 Fig. 10. Critical frequencies −μcr , norms P cr and peak amplitudes |ψ0 |cr at local maximum of P (|μ|) curve as the functions of nonlinear power γ . Here, markers denote numerical results. For −μcr and P cr , solid lines denote variational estimate, 2 dashed lines denote continuum estimate results. For |ψ0 |cr , the solid line is the results obtained by Eq. (9). t 1 = t N = 1, N = 2.

ond stage of compare the of breathers. tude for the

the MI for k = 0 plane wave. To verify this link, we amplitude of plane wave with the peak amplitude Remember that in Eq. (9) we obtain critical amplisecond stage plane wave MI, which is given by the

expression | A cr |2(1) =



is reduced to | A cr |2 =

Fig. 9. (a) Profiles of breather indicated in Fig. 7(a) at γ = 0.8; (b) decay dynamics of stationary breather with weak noise at γ = 0.8, μ = −5.34; (c) discrete Fourier transform of stationary breather for different μ at γ = 0.8. Throughout t 1 = t N = 1, N = 2.

weak noise, the final state oscillates around the μ = −6.03 solution. The rearrangement towards μ = −5.13 solution, however a little harder, can also be found. These three solutions showed in Fig. 9(a) also reveal the tendency of localizing to its central site with increasing |μ|. However, coinciding with the critical frequency μcr where the P (|μ|) curve has its maximum and the instability sets in, the wave function of breather becomes a stair-step shape. This can be seen in discrete fourier transform ϕ (k) of breathers shown in Fig. 9(c) as well. For small |μ| the ϕ (k) is localized in the vicinity of k = 0 in quasi-momentum space, and gradually extends with the localization of wave function in direct space. At the critical frequency μcr = −5.22 marked by black solid line, ϕ (k) changes drastically, arises first near k = ±π followed by a quick extending into the full quasi-momentum space. The fact that the instability emerges at k = ±π reminds us of the second stage MI of k = 0 plane wave which also arises from k = π for N = 2. This suggests that the destabilization of breather branch at the maximum of its P (μ) curve is linked to the sec-

2t 1

2 π γ sin m N



2t 1

 γ1

1 γ

γ

, for N = 2, the expression

. Admitting this link, the breather

critical amplitude should be close to the expression. In order to check these, we provide numerically obtained critical breather amplitudes as a function of γ , comparing them to the critical amplitudes of mudulational instability in Fig. 10. We see that these two amplitudes display quite good agreement over a wide intermediate range but some deviation at small and large values of γ due to the distinction of breather and plane wave solutions. This link was first noticed in the model with cubic nonlinearity [22] where γ = 1, here we provide a verification at other values of γ . A remarkable feature is that this link enable us to estimate the critical frequency and norm at the local maximum of P (|μ|) curve. To do this, we use the variational and the continuum approximation to obtain the analytical expressions of peak amplitude and norm as a function of frequency. Now from the above discussion, it is reasonable to assume that the breather turn unstable when the approximatively obtained squared peak amplitude equals



1

 γ 2 π . This yields the critical values of freγ sin N quency μcr and norm P cr . The detailed derivations are given in Appendix A, we only provide results here. The critical frequency and norm obtained by continuum approximation are to | A cr |2(1) =

2t 1

μcr ≈ −2(t 1 + t N ) −  P cr ≈

2t 1 sin2 ( π ) N

2−γ

γ

1  γ

t 1 sin2

π  N

t1 + N 2t N t 1 sin2 ( π ) N

γ

(12)

,

(13)

.

By the variational approximation we get the critical frequency and norm which are slightly different from the continuum approximation

μcr ≈ −2(t 1 + t N ) −  γ

P cr ≈

2t 1

2+γ

γ

2 (2 + 1

γ +γ)

γ (2 + γ1 + γ )  γ

1 + N2

tN t1

t 1 sin2

sin

π  N

π  N

,

(14)

JID:PLA AID:25076 /DIS Doctopic: Nonlinear science

[m5G; v1.235; Prn:25/04/2018; 14:56] P.7 (1-9)

W. Qi et al. / Physics Letters A ••• (••••) •••–•••

 ×

γ (2 + γ1 + γ )(t 1 + N 2t N ) t 1 sin2 ( π ) N

− 2N ( N 2 − 1)  ×

 γ −2 1

t1

tN

γ t1 + N 2t N

1

sin2

γ (2 + γ + γ )(t 1 + N t N )

 t 1 sin2 π N

2

π 

(15)

N

 γ −2 1 .

The estimations of Eqs. (12)–(15) as a function of γ are shown in Fig. 10, comparing to the numerically obtained results at the maxima of P (|μ|) curve respectively. The estimations are in quite good agreement with numerical results in general, whereas at small and large values of γ there are some deviations as predicted above. It shows that γ has a significant impact on the critical frequency and norm. To analyze the impact on the critical frequency, we shall notice that γ and μ are related by amplitude. In the MI analysis, it is shown that critical amplitude of plane wave is a decreasing function of γ . Furthermore, as shown in Appendix A, the amplitude of breather and frequency μ change consistently. The only blind spot is the relation between the amplitude of plane wave and the amplitude of breathers, which is now verified by above reasoning. Taking all these factors into consideration, it is expected that the critical values of μ is a decreasing function of γ , which is fit with the result shown in the figure. The influence of γ on the critical norm is not only more complicated but also significant. In addition, although the results obtained by variational approximation and continuum approximation capture the overall tendency, we also see notable deviation for small and large γ . 5. Conclusion We have investigated the power-law nonlinearity of DNLS equation through the modulational instability analysis and localized breather solutions. As a result, the impact of nonlinearity power γ on the MI is similar to the impact of amplitude, that is, the nonlinear-power induces the additional instability region in the same order as the amplitude, and for quasi-momentum k = 0 both amplitude and nonlinear power induce the two-stage instability. This similarity also means that with increasing γ , the critical values of amplitude for the occurrence of MI become smaller and smaller. For localized breather modes, the numerical simulations reveal that branch of breathers is monostable in weak nonlinearity, but when the nonlinearity power exceeds a critical value, the P (|μ|) curve exhibits a local maximum before a local minimum, so that three different solitons coexist in certain interval of values of P , the two solitons with the smallest and greatest frequency are stable, while the third (intermediate) one with negative derivative dP is unstable, leading to the bistability and instability feature. d |μ| We also find that the long-range hopping term in DNLS equation can dramatically change the dynamics properties of breathers, for example, the critical value of nonlinearity power is 0.64 for t N = 1, much smaller than the critical value 1.34 from the model without long-range hopping, i.e., t N = 0 [18]. Furthermore, we analytically give the bistability conditions by using variational approaches and the continuum estimates, and our analytical results are confirmed by the numerical simulation. Therefore, we get the conclusion that combining the long-range hopping and power-law nonlinearity dramatically change the breather dynamics in the DNLS equation. Our results are valuable not only for theoretical but also for experimental. For theoretical view, we systematically investigate the MI condition and the breathers stability properties in a most general DNLS equation, the results can be wildly implicated in the relevant system if the system has the similar mathematical description. For experimental aspect, (i) we find that the dynamics

7

properties of ultra-cold atoms trapped in the helicoidal optical trap are dramatically influenced by the atomic interactions, due to the fact that every parameter can easily be tuned in ultra-cold atomic experiment, we expect our theoretical results will stimulate the experimental study in this direction. (ii) We also find that the base-pair interaction dramatically influence the DNA nonlinear dynamics, we propose using this mechanism to investigate the energy storage and transport in DNA molecules or cellular gene expression [23] experimentally in the future. Acknowledgements This work is supported by the National Natural Science Foundation of China under Grant Nos. 11647017, 11647049, 11405100 and Science Research Fund of Shaanxi University of Science and Technology, China, under Grant No. BJ16-03. Appendix A. Variational approaches and continuum estimates In the main text we analyzed the relationship between second stage MI of plane wave and the destabilization of localized breather solutions. Admitting this idea, it is reasonable to assume that the critical peak amplitude of breathers is equal to the second stage MI critical amplitude of k = 0 plane wave. Since peak amplitude, norm and frequency of a breather branch are linked to each other, this equivalence enable us to estimate the critical values of frequency and norm. In this appendix, we employ two approximation methods and provide the detailed derivation of approximate relations between peak amplitude, norm and frequency, then give the critical expressions of Eq. (12)–(15). The variational results can also be used to analyze the bistability condition. The first method we give is the variational approximation. The stationary DNLS equation (2) can be obtained from the Lagrangian

L=



− μ|ψ j |2 −

j

1

γ +1

|ψ j |2(γ +1)

∗

∗

− t 1 ψ j (ψ j +1 + ψ j −1 ) − t N ψ j (ψ j + N

+ ψ j−N ) .

(A.1)

In this approximation, instead of the exact breather solution, we use a tractable variational ansatz with a finite number of parameters. For the breather wave packets localized near k = 0 in Fourier space, a commonly used tractable ansatz is ψ j = B exp(−a| j |), where a and B are real variational parameters [22,27,28]. The norm of this ansatz is

P=



|ψ j |2 = B 2 coth(a),

(A.2)

j

to obtain the relation between B and μ, we shall explore the expression of P and a towards μ. Inserting this ansatz into Eq. (A.1) and eliminating the amplitude B via Eq. (A.2), yield

L (a, P ) =

  −t 1 − t N ea(1− N ) − μ P

2P cosh(a)

− 2P t N ( N − 1)e −aN tanh a

(A.3)

P γ +1 coth[(γ + 1)a] − , γ +1 γ +1 coth

a

using the variational equations ∂∂ PL = 0, ∂∂ aL = 0, we get

μ = − 2t 1 sech(a) − 2t N e−aN (1 + N tanh a) γ coth[(γ + 1)a] −P

γ +1

coth

a

,

(A.4)

JID:PLA

AID:25076 /DIS Doctopic: Nonlinear science

[m5G; v1.235; Prn:25/04/2018; 14:56] P.8 (1-9)

W. Qi et al. / Physics Letters A ••• (••••) •••–•••

8

Fig. 11. (a) P , |ψ0 |2 , a of localized breather as the functions of frequency μ for γ = 0.8 obtained by numerical results (solid lines) and variational estimate (markers); (b) wave function of localized breather for γ = 0.8, μ = −5 obtained by numerical results (solid lines) and variational estimate(markers); (c) P , |ψ0 |2 of localized breather as a function of frequency μ for γ = 0.8 obtained by numerical results (solid lines) and continuum estimate (markers); (d) wave function of localized breather for γ = 0.8, μ = −5 obtained by numerical results (solid lines) and continuum estimate (markers). Throughout t 1 = t N = 1, N = 2.

and

Pγ =

2 sech(a) tanhγ +3 (a) coth[(γ + 1)a] csch2 a − coth a csch2 [(γ + 1)a]

× [t 1 + Nt N e

(−aN )

(A.5)

( N + tanh a) cosh(a)].



We put the Eq. (A.5) into Eq. (A.4) to obtain μ(a). Noting that the breathers are broad and the decay parameter a is small in the region we are interested in, we expand μ(a) to second order in a

μ(a) = −2(t 1 + t N ) −

(t 1 + N t N )(2 + γ ) 2

γ

a2 .

(A.6)

By solving this equation for a, yield

 a(μ) =

This is the estimation of peak amplitude given by variational approximation. Now by the reasoning at the beginning of this Appendix, the instability of breather is expected to set in when its peak amplitude crosses

B2 =

2t 1

γ

sin2

(A.7)

which is the expressions of a towards μ. Inserting Eq. (A.7) into Eq. (A.5) and expanding in a once again, we get

1 P γ (a) ≈ γ a

 2(t 1 + N 2 t N )(2 + 3

+ 2( Nt N − N t N )(2 +

1

γ 1

γ

+ γ )a2  + γ )a

3

(A.8)

2γ (2 + γ + γ ) [−μ − 2(t 1 + t N )]. (2 + γ )

(A.9)

,

(A.10)

1 2M

∂x2 − |ψ(x)|2γ ψ(x),

(A.11)

where ψ(x) is the continuum wave function in replace of dis1 crete function ψ j , M is the effective mass with M = μ0 (k)|k=0 = 2 ˜ = μ + 2(t 1 + t N ). For μ ˜ < 0, Eq. (A.11) has 2(t 1 + N t N ) and μ approximate soliton solutions of the form ψ(x) = β sech(α x). By inserting this solution into Eq. (A.11) and expanding it to second order in x we have a set of equations. By solving this set, we get the parameter α and β :



1

B 2γ (μ) ≈

μ˜ ψ(x) = −

.

Combining Eq. (A.7) with Eq. (A.8) gives the expressions of P towards μ. Inserting the P and a expressions into Eq. (A.2), we get

N

from which Eq. (A.9) predicts critical frequency approximately as in Eq. (14), inserting this into Eq. (A.8) yields Eq. (15). Another approach to estimate the critical parameter of breathers is continuum approximation. Here, we study the continuum counterpart of DNLS equation of the form

 [−μ − 2(t 1 + t N )]γ , (t 1 + N 2t N )(2 + γ )

 π  γ1

α=

˜ μγ , (−2 + γ )(t 1 + N 2 t N )

 β=

˜ 2μ

γ −2

(A.12)

 γ1 .

(A.13)

JID:PLA AID:25076 /DIS Doctopic: Nonlinear science

[m5G; v1.235; Prn:25/04/2018; 14:56] P.9 (1-9)

W. Qi et al. / Physics Letters A ••• (••••) •••–•••

The norm of these solutions is thus given by

∞ dx|ψ(x)|2

P= −∞



=2

˜ 2μ γ −2

1  γ

(A.14)

(−2 + γ )(t 1 + ˜ μγ

N 2t

N)

.

[2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13]

Eq. (A.13) is the estimation of the peak amplitude given by the continuum approximation. Equating the amplitude from Eq. (A.10), leads to the continuum approximation of critical frequency as in Eq. (12), inserting back into Eq. (A.14) yields Eq. (13). Fig. 11(a)(b) compare P , |ψ0 |2 , a curve and a localized solution obtained by variational approximation with the results from numerical simulation respectively. In contrast, Fig. 11(c)(d) compare P , |ψ0 |2 curve and a localized solution obtained by continuum approximation with the results from numerical simulation respectively. It is shown that both approximate methods are relatively in good agreement with numerical results, we also find that near the local maximum of P (|μ|) curve with γ = 0.8, the variational approximation tends to underestimate the norm but overestimate the peak amplitude, while the continuum approximation tends to overestimate the norm but underestimate the peak amplitude, which are in agreement with Fig. 10.

[24] [25] [26]

References

[27] [28]

[1] P.G. Kevrekidis, The Discrete Nonlinear Schrödinger Equation: Mathematical Analysis, Numerical Computations and Physical Perspectives, Springer, Berlin, Heidelberg, 2009.

[14] [15] [16] [17] [18] [19] [20] [21] [22] [23]

9

T. Holstein, Ann. Phys. (N.Y.) 8 (1959) 325. O. Morsh, M. Oberthaler, Rev. Mod. Phys. 78 (2006) 179. J.K. Xue, A.X. Zhang, Phys. Rev. Lett. 101 (2008) 180401. M. Salerno, V.V. Konotop, Y.V. Bludov, Phys. Rev. Lett. 101 (2008) 030405. A. Trombettoni, A. Smerzi, Phys. Rev. Lett. 105 (2010) 050405. S.F. Mingaleev, P.L. Christiansen, Y.B. Gaididei, M. Johansson, K.Ø. Rasmussen, J. Biol. Phys. 25 (1999) 41. N. Molkenthin, S. Hu, A.J. Niemi, Phys. Rev. Lett. 106 (2011) 078102. D. Reitz, A. Rauschenbeutel, Opt. Commun. 285 (2012) 4705. T. Dauxois, Phys. Lett. A 159 (1991) 390. J. Stockhofe, P. Schmelcher, Phys. Rev. A 92 (2015) 023605. G. Wang, J.P. Huang, K.W. Yu, Opt. Lett. 35 (2010) 1908. F. Dreisow, G. Wang, M. Heinrich, R. Keil, A. Tünnermann, S. Nolte, A. Szameit, Opt. Lett. 36 (2011) 3963. W. Yan, Commun. Theor. Phys. 59 (2013) 643. E.B. Kolomeisky, J.P. Straley, Phys. Rev. B 46 (1992) 11749. X.D. Bai, J.K. Xue, Phys. Rev. E 86 (2012) 066605. W. Qi, J.K. Xue, Phys. Rev. A 81 (2010) 013608. J. Cuevas, P.G. Kevrekidis, D.J. Frantzeskakis, B.A. Malomed, Physica D 238 (2009) 67. J.K. Pedersen, D.V. Fedorov, A.S. Jensen, N.T. Zinner, J. Phys. B 47 (2014) 165103. J.K. Pedersen, D.V. Fedorov, A.S. Jensen, N.T. Zinner, J. Phys. B 49 (2016) 024002. T. Zhang, G.B. Jo, Sci. Rep. 5 (2016) 16044. J. Stockhofe, P. Schmelcher, Physica D 328 (2015) 9. P. Maniadis, B.S. Alexandrov, A.R. Bishop, K.Ø. Rasmussen, Phys. Rev. E 83 (2011) 011904. S. Zdravkovi, J. Nonlinear Math. Phys. 18 (2011) 463. Z. Rapti, P.G. Kevrekidis, A. Smerizi, A.R. Bishop, Phys. Rev. E 60 (2004) 017601. N.G. Vakhitov, A.A. Kolokolov, Izv. Vysš. Uˇcebn. Zaved., Radiofiz. 16 (1973) 10120; Radiophys. Quantum Electron. 16 (1973) 783. B. Malomed, M.I. Weinstein, Phys. Lett. A 91 (1996) 220. J. Cuevas, P.G. Kevrekidis, D. Frantzeskakis, B. Malomed, Physica D 67 (2009) 238.