Modulational Instability and Discrete Localized Modes in Two Coupled Atomic Chains with Next-Nearest-Neighbor Interactions

A pair of one dimensional atomic chains which are coupled via the Klein-Gordon potential is considered in this study, with each chain experiencing both nearest and next-nearest-neighbor interactions. The discrete nonlinear Schrödinger amplitude equation with next-nearest-neighbor interactions is thus derived from the out-phase equation of motion of the coupled chains. This is achieved by using the rotating wave approximations perturbation method, in which both the carrier wave and envelope are explicitly treated in the discrete regime. It is shown that the next-nearest-neighbor interactions greatly modifies the region of observation of modulational instability in the atomic chain. By exploring the discrete Hirota-Bilinear method, we obtain the discrete one-soliton solution which is localized around the origin and structurally stable because it conserves it form as time evolves. However when the atomic chain is purely subjected to a symmetric coupling potential, we observe a structurally unstable discrete excitation that changes into an up-and-down asymmetric localized modes; both in the presence and absence of next-nearest-neighbor interactions. Results of numerical simulations clearly depicts the long term evolution of these discrete nonlinear excitations, that evolve from symmetric to asymmetric localized modes in the atomic chain.


Introduction
Modulational instability (MI) is an ubiquitous phenomena in nonlinear dynamics that have been observed in a broad spectrum of systems ranging from dual-core nonlinear optical fiber [1], nonlinear LC self-dual circuits [2], fibers with random polarization mode dispersion [3], negative refractive index medium [4], discrete optical nonlinear arrays [5], Bose-Einstein condensates [6,7], deoxyribonucleic acid (DNA) chains [8], plasmas [9], among other nonlinear physical systems. In fact, MI is generally observed in systems where there is a dynamic balance between nonlinearity and dispersion, and characterized by small amplitude and phase perturbations growing exponentially.
Since the discovery of solitary waves by John Scott Russell, solitons have been experimentally observed in many spatially discrete physical systems, such as granular materials, nonlinear Hamiltonian lattices [10], electrical transmission lines, and neural networks [11][12][13][14]. The recurrence phenomenon was originally observed in the Fermi-Pasta-Ulam lattice [15,16]. The eventual formulation of the soliton concept by Zabusky and Kruskal to elucidate on the recurrence phenomenon, triggered enormous interest in scientists around the world to investigate on this fascinating soliton entity [16]. Of particular interest is in discrete regimes, where the role and consequences of MI have been extensively studied within the frame work of discrete nonlinear Schrödinger (DNLS) amplitude equations [17][18][19][20][21][22][23][24][25]. This generally leads to the observation of discrete breathers, that were first identified by Takeno and Sievers in a perfect anharmonic lattice [26,27]. Physically, the discrete MI emanates from the interplay between nonlinearity and spatial discreteness, and a potential mechanism for the formation of nonlinear localized modes. This fully accounts for energy localization in discrete nonlinear lattices [28].
Numerous research activities in the past have so far concentrated on the properties of MI and discrete breathers in the DNLS equation, with nearest-neighbor (NN) interactions. However to effectively study the dynamics of most real physical systems, necessary measures must be taken to incorporate long-range coupling interactions, which may play an important role in explaining some important physical processes in the discrete lattice. In this light, we see that charged groups of molecules which are mediated by Coulomb interactions in DNA systems are governed by long range interactions [29][30][31][32][33]. The dipole-dipole interaction characterized by a typical long range power law dependence within molecular subunits, mimics a long-range interaction that is responsible for exciton and vibron energy transfers in biomolecules and crystals [34]. Also, in the dynamics of excitons and photons interaction in semiconductors, the polariton effects naturally leads to long range interaction radius in molecular crystals [29,30,35]. Furthermore an optical zigzag waveguide array [36], depicts a classical zigzag lattice model which have been successfully applied to improving the next-nearest-neighbor (NNN) coupling in some planar lattices. The NNN coupling is thus enhanced, and eventually leads to a modification of the region of linear dispersion curve. Consequently, there is an increase in the range of allowed frequencies for the effective propagation of electromagnetic waves in photonic crystal waveguide and coupled resonator optical waveguide in the form of robust soliton entity with minimal distortions [36]. The NNN interactions of peptide groups along a polypeptide chain enhances stability in the secondary structure of −helix and − sheet regions of proteins [31,34]. Finally results from studies carried out, strongly suggests that the NNN coupling accounts for interesting phenomena in discrete lattices. These include the identification of the staggering and unstaggering stationary localized breather modes [20], observation of pulse propagation in nonlinear photonic crystal wave guides [36], existence of bright and dark discrete breathers in 4 nonlinear lattice [37], linear increase of the input power required for soliton formation in lattices [38], and the reversal phenomenon of the propagation direction of the Brillouin zone boundary mode [39]. This study mainly captures a generic model of an atomic lattice and theoretically investigate on the effects of the NNN coupling, with potential applications in quantum-information storage [40] and DNA systems [41] serving like a driving force behind our current investigation.
In this work, we consider two atomic chains which are nonlinearly coupled via Klein-Gordon potential, and experiencing both nearest and NNN interactions. Consequently in Sect. 2, we obtain the Hamiltonian of the coupled chains and explore the Hamilton-Jacobi formalism to obtain the appropriate equations of motion. By using the rotating wave approximation scheme in which both the amplitude and phase of the wave are constraint in the discrete regime, the DNLS amplitude equation with NNN interaction is obtained. In Sect. 3 we discuss the stability of carrier plane wave solutions of the DNLS equation with NNN coupling. Regions of occurrence of MI in the parameter space are identified, with the influence of NNN interaction on this parameter space demonstrated. Analysis in Sect. 4 focuses on analytically obtaining discrete breather mode solutions of the amplitude equation, with the help of the Hirota's binilearization method. We also utilize a numerical integration based on the standard fourth order runge-Kutta scheme, to investigate the effects of the NNN interactions on the evolution of discrete localized modes in the atomic chain. Finally we give a brief summary in Sect. 5. Figure 1 depicts the mechanical model of the two one-dimensional atomic chains under consideration. We assume that each atom on each chain is subjected to linear interactions with its nearest and NNN couplings on the same chain. The atoms are equally spaced along each chain in the rest position. The chains are lying parallel to each other and nonlinearly coupled by nonlinear springs represented by the hollow rectangular blocks, and the coils represent linear springs. In forming kinetic equations for the motion of such a lattice, we assume that all atoms have identical mass m, the spring constant for NN interaction is k 1 , while that for NNN interaction is k 2 . The energy stored in the nonlinear springs is denoted by the potential V(x n − n ) , where the atomic displacements of atoms in the chains are x n and n . Physically, the system can be characterized by the following Hamiltonian with the Klein-Gordon coupling potential given by where 1 , 2 , 3 are constants. The Klein-Gordon potential emanates from a more general potential known as Morse-potential. It is given by V( ) = D 0 (exp(a ) − 1) 2 , and can be approximated as

The Model
The terms in equation (3) are derived by using the Taylor series expansion, with 1,2,3,4,5 being the quadratic, cubic, quartic, ... coefficients respectively. The potential is symmetric for 2 = 4 = 0; hence termed 4 potential for 1 ≠ 0, 3 ≠ 0, 5 = 0 as in reference [37], or a more general case of 6 potential for coefficients [42]. However for non-varnishing values of 2 and 4 , the potential wells are asymmetric as with equation (2). The dynamics of periodic base pairs opening in a finite stacking enthalpy DNA was recently investigated for non zero values of all the coefficients 1 , 2 , 3 , 4 , 5 [43].
For k 2 = 0 , the Hamiltonian (1) of the coupled chains is physically similar to the Peyrard-Bishop model that governs the dynamics of localized large amplitude excitations in the DNA molecule [44,45]. By using the Hamilton-Jacobi formalism, the equation of motion for x n and n thus reads  It is convenient for us to use the in-phase motion u n = x n + n , and out-phase motion v n = x n − n , in order to transform equation (4) into is a linear differential-difference equation with the usual plane wave solutions (phonons), while (5b) is the equation that contains the nonlinear terms and can only be solved using perturbation techniques. By considering linear oscillations of the out-phase motion v n , at frequency and wavenumber q, with the chain of lattice spacing a, we can easily obtain the dispersion relation given by

Fig. 2
Linear dispersion curve of the chain according to Eq. (6) for Ω 2 0 = 0.032, K 1 = 1.00, and the values of K 2 are varied modify the linear dispersion curve of wave propagation in the atomic chain. For higher values of K 2 , the maximum value of wave frequency is no longer at the edge of the Brillouin zone q = ∕a , but is now shifted to a new position with wavenumber q o given by [37] The system of equation (5b) for K 2 = 0 has been extensively studied within the framework of nonlinear Schrödinger amplitude equation, and myriad of breather and envelope-soliton solutions obtained in the low-amplitude limit [18,19,42]. For example, Kivshar et al. exploited the rotating-wave approximation method to derive the DNLS equation [18,19], with solutions in the form of bright and dark-profile localized structures.
In a highly discrete regime where the nearest and NNN coupling forces between particles are weak (i.e. Ω 0 >> K 1 and Ω 0 >> K 2 ), the equation that governs the outphase motion (5b) can now be re-written as where << 1 is a small perturbation parameter. The rotating-wave approximation scheme where both the carrier wave and envelope are strictly treated in the discrete regime [26,46,47], seems to be the most suitable method to explore. Consequently we now analyze slow temporal variations of the wave envelope by considering solution of equation (8) in the form [18] where the discrete amplitudes A n , B n , C n and their complex conjugates A * n , C * n , are functions of T = 2 t , and we are strictly operating at the lower cut-off mode frequency Ω 0 . Upon substituting the ansatz (9) into equation (8) and simplifying, we now obtain the following equations in their respective harmonics: Finally, upon substituting relations (12) and (11) into equation (10), leads us to where Equation (13) is the DNLS equation with NNN interactions. For P 2 = 0 , equation (13) can equally govern the propagation of laser wave packets in optical media, dynamics of matter wave solitons in Bose-Einstein condensates, nonlinear excitations in DNA molecules, among other potential applications [48][49][50]. We will investigate on the MI of the amplitude equation (13) in the next section, in order to effectively map out the appropriate regions for the observation of localized modes in the atomic chain.

Modulational Instability Analysis
MI generally precedes the generation of localized wave modes in discrete lattices, which has been experimentally identified in electrical transmission lines [51]. The generation of higher order-rogue wave solutions for the generalized discrete Hirota equation was equally shown to be driven by the process of MI [52]. Consequently, we now explore the stability properties of equation (13) which possesses exact plane wave solutions where A 0 is the amplitude, is the angular frequency, and q is the wavenumber. By substituting (15) into the discrete amplitude equation (13), leads to the nonlinear dispersion relation where 0 (qa) = 4P 1 sin 2 qa 2 + 4P 2 sin 2 (qa). To examine the linear stability of the solution (15), we consider small perturbations B n (.i.e |B n | << A 0 ) on the amplitude of the plane wave such that Using Equations (13), (16), and (17), we obtain the following linearized equation in B n (T) as where the asterisk denotes complex conjugation.
We now assume a general modulation solution to equation (18) of the form where Q and Ω respectively represent the wave vector and the frequency of the linear modulation. From the systems (18) and (19), we can easily deduce the dispersion relation for the perturbation wave as where Γ = P 1 cos(qa)sin 2 (Qa∕2) + P 2 cos(2qa)sin 2 (Qa) . The imaginary part of Ω can only be realized when Γ(Γ − |A 0 | 2 ) < 0 , hence leading to the exponential growth in the perturbation B n (T) which triggers the generation of localized modes in the coupled chains. Figure 3 depicts the regions in the (qa, Qa) plane for identifying such modes, which is realized by computing the magnitude of the imaginary part of equation (20). It should be noted that for the special case P 2 = 0 , equation (20) reduces to result obtained by Kivshar [18].
As shown in Fig. 3, the exponential growth in small perturbations B n (T) generally occurs at small amplitudes plane waves with the NNN coupling parameter P 2 greatly modifying their area of occurrence. This is inextricably linked to the dispersion relation (16). To comprehensively analyze this, we consider small values of Qa and simplifying Γ in order to have It is clear that for small values of Qa and positive values of the curvature of the dispersion ( i.e. �� 0 (qa) > 0 ), then Γ > 0 . Hence MI is only possible when > 0 (i.e. 10 2 > 9Ω 2 0 ) ; since Γ ∝ (Qa) 2 , a minimal amplitude A 0 equally triggers infinitesimal-Qa instability. On the other hand, if �� 0 (qa) < 0 , then minimal-Qa instability is only possible at small amplitudes provided 10 2 < 9Ω 2 0 . We can now confidently say that the regions of the Brillouin zone in which �� 0 (qa) > 0 , are more vulnerable to MI provided 10 2 > 9Ω 2 0 . For the special case in Fig. 3(a) where P 1 = 1.00, P 2 = 0.00, and |A 2 0 | = 1.00 , MI is rare to observe for any small values of Qa probably because �� 0 (qa) < 0 . However by progressively increasing the values of P 2 from 0.20 to 1.00 as in Fig. 3b-f, we increase our chances of observing MI in the chain because the yellowish area(s) also progressively stretches downward. We also note that for these small values of Qa; two major distinct regions of MI progressively emerges with boundaries around qa = 1.00 , as the values of P 2 is progressively increased. However, the MI phenomena for the region where qa < 1.00 , is weaker compared to their counterpart region qa > 1.00.
We further investigate a more interesting MI phenomena at wavenumber qa = 0.00 and varying values of Qa. At this stage, we have since P 1 , P 2 are assumed to be positive constants. Consequently, for MI to be observed at qa = 0 , it requires that 10 2 > 9Ω 2 0 , with the amplitude of plane wave constraint in the interval Values of Qa with smaller 0 (Qa) become unstable for smaller plane wave amplitudes, and the first instability always emerging at the global minimum where 0 (Qa) = 0 for Qa = 0.00 . There is also an additional local minimum in the dispersion curve that will induce another region of MI in the chain. This will occur at Qa = , which corresponds to another critical plane wave amplitude given by A careful examination of Fig. 3c-f, clearly indicates two regions of MI for qa = 0.00 , which is consistent with the analytical predictions.

Discrete Localized Modes
Enormous effort has been made to obtain analytical solutions to a number of discrete nonlinear partial differential equations. Some prominent analytical methods includes: quadrature method, ansatz method, method of separation of variables, Painlevé method, Hirota method, IST method, B äcklund transform method, Darboux transformation method among others [11,18,51,[53][54][55][56]. Furthermore, a more powerful analytical method like the discrete generalized (m, N − m)-fold Darboux transformation technique; have already been developed to obtain two-component localized wave solutions for a reduced semi-discrete two-component coupled system on a two-chain lattice [57]. By exploiting this latter method, breathing-soliton and singular rogue wave solutions were obtained for a discrete nonlocal coupled Ablowitz-Ladik equation of reverse-space type [58]. Some of these methods are more appropriate with some mathematical problems under specific conditions. A variety of intrinsically localized stationary solutions which are also termed discrete breathers or discrete bright solitons, can be obtained for the DNLS amplitude equation (13), by exploiting the Bilinear Hirota analytical technique as elaborated below.

Bilinear Hirota Method
Concretely, the Hirota's bilinear method was proposed by Hirota in 1971 to solve the Korteweg-de Vries equation [53]. It does not involve any sophisticated technique, but only differentiation. This method can construct multiple-soliton solutions as well as one-soliton solutions effectively using exponentials. Hence it also helps to study discrete equations as will be demonstrated here. The principle of the Hirota bilinear method lies on the implementation of the Hirota bilinear operator, which transforms the nonlinear evolution equations into several coupled bilinear equations. This act breaks down the original complicated equation into a series of relatively simple equations. A detail study on the Hirota bilinear transformation to different nonlinear evolution equations can be found in Ref. [56]. The Hirota's bilinear operator (or D-operator) is defined as where s and l are positive integers. Depending on the nature of the physical problem in question, several modifications and improvements are made on the D-operator in system (25), to obtain appropriate solutions.
The most appropriate discrete localized modes of equation (8) is now obtained by using the ansatz (9), in order to have with Φ = bn + Ω 0 t + 2 2(P 1 + 2P 2 )cosh(a)cos(b) − 2(P 1 + P 2 ) t. We have set = a + ib and = 1.00 , without loss of generality. We thus obtain the profile of the solution given in Fig. 6 for the asymmetric coupling potential (i.e. for = −1.00 ), while Fig. 7 depicts the profile in the case of symmetric coupling potential (i.e. for = 0.00 ). For another case of asymmetric coupling potential in which = 1.00 ; we obtain the profiles in Fig. 8 that generally shows dark-like localized modes because of the constant dip. Concretely for = −1.00 and b = 1 + 0.001i as shown in Fig. 6, we observe a structurally stable stationary discrete bright solitary wave centered around the origin. It is a dynamic structure that mimics a localized standing wave in a breathing mode, owing to the slight variation in the wave amplitude as time evolves. In fact, the wave amplitudes in Fig. 6d-f for P 2 = 1.00 supersedes those of their counterparts in Fig. 6a-c for P 2 = 0.00 ; indicating that the NNN coupling enhances the amplitude of the localized modes. The nature of the coupling potential greatly influence the form of the localized mode in the atomic chain. This is reflected in Fig. 7 for = 0.00 , where the nonlinear excitation evolve into a discrete localized mode with up-and-down asymmetry [13,60]. For P 2 = 1.00 as in Fig. 7d-f, there is a rapid transformation of the up-and-down asymmetric localized mode with minimal change in amplitude. As time increases, we observe that the amplitude of the asymmetric discrete solitary wave drastically reduces as shown in Fig. 7a-c for P 2 = 0.00 . This loss may be as a result of thermal vibration in the lattice with the energy of the localized modes carried away as small amplitude linear waves in the form of radiation [18]. Lastly, the wave profile observed in Fig. 8 for = 1.00 is just (33) A n (T) = 8(P 1 + P 2 )sinh 2 k+k * 2 + 8P 2 cosh(k + k * ) exp J the reflection about the spatial axis of the wave pattern identified in Fig. 6. The darklike localized modes identified in the atomic chain as in Fig. 8, may signal the absence of MI. This therefore underscores the crucial role played by the asymmetric parameter , in the dynamics of discrete soliton of the two coupled chains.

Numerical Analysis
All the solution profiles so far obtained are approximations because of the assumptions made in the analytical calculations. We now conduct direct numerical simulations of equation (8) to determine the outcome of the long term evolution of the localized discrete solutions derived in the previous subsection. Consequently, we initially split the second order discrete equation (8) into two first order discrete ordinary differential equations. We then define the differential equations as a function that takes one vector of variables in the differential equation plus a time vector as an argument. Lastly, we return the derivative of the vector and use the built-in MAT-LAB solver ode45; that is based on the fourth order Runge-Kutta scheme. The initial displacement is given by while the initial velocity is simply obtained by differentiating solution (34) with respect to time and set t = 0.0 . This leads to the long term evolution of discrete localized modes as shown in Fig. 9. At the initial stage (i.e. t = 0 ), we observe a bright discrete pulse in Fig. 9 that is symmetric about the origin. As time changes, the localized discrete bright pulse becomes structurally unstable as it evolves from a symmetric to asymmetric wave profile around the origin. A close examination between Figs. 9 and 6d-f reveal the same structural features of the discrete localized wave profiles, as they share the same parameter values. In fact the amplitudes of the localized modes in both figures fluctuates around the numerical value of 80, clearly showing the consistency between the analytical curves and long term numerical profiles of the discrete localized modes.

Conclusion
For integrable systems, the word "soliton" has a mathematical meaning of maintaining constant amplitude wave profile with it form unchanged even after collision. For non-integrable systems like the DNLS equation with NNN coupling earlier derived, the equation merely support but solitary waves (localized modes) which can loosely be referred to as soliton because of the localized and long-lived nonlinear energy packet. Due to non-integrability, a soliton may lose energy due to radiation [18]. It can equally acquire an internal mode as a result of a degree of freedom corresponding to shape deformation [61,62], which is a clear indication that solitons may carry with them an internal store of energy. We have investigated on the MI properties and discrete localized modes in two chains which are non-linearly coupled via Klein-Gordon potential. Each chain independently experiences both nearest and NNN interactions, with the model mimicking the Peyrard-Bishop DNA model. The DNLS equation with NNN interactions was derived from the out-phase equation of motion by exploring the rotating wave approximation, in which the wave amplitude and phase are constraint in the discrete regime. Within the frame work of the linear stability analysis, we calculated the growth rate of the MI, and graphically captured the areas of their existence. It was established that the introduction of the NNN coupling parameter P 2 , significantly modified the shape of the areas of MI. We further explored the Hirota-Bilinear method to obtain discrete one-soliton of the DNLS equation with NNN interactions. Our results showed that for symmetric coupling potential (i.e. = 0.0 ) and in the absence of NNN interaction parameter P 2 , our system supports a discrete asymmetric localized mode which is structurally unstable as time (35) v n (0) = 2 8(P 1 + P 2 )sinh 2 (a) + 8P 2 cosh(2a) exp(an) 8(P 1 + P 2 )sinh 2 (a) + 8P 2 cosh(2a) + exp(2an) cos(bn) + 2 2 Ω 2 0 8(P 1 + P 2 )sinh 2 (a) + 8P 2 cosh(2a) exp(an) 8(P 1 + P 2 )sinh 2 (a) + 8P 2 cosh(2a) + exp(2an) 2 cos(2bn) unfolds. However by introducing the P 2 parameter, the discrete localized modes evolve into standing wave profile with up-and-down asymmetry. We equally identified structurally stable discrete localized mode in the form of a dark soliton when the asymmetric parameter assumes a positive value. Results of numerical investigations mainly confirmed the long term evolution of the discrete symmetric pulse around the origin to asymmetric pulse structures. This study provides an elegant mathematical approach that can be used to rigorously study highly localized discrete modes in a DNA molecule. The dynamics of DNA is increasingly being investigated in order to adequately comprehend the formation of genetic codes, and intrinsic processes of replication and transcription [41,43,63,64]. Discrete solitons have also been shown to be very robust entities. This is the case with elastic soliton interactions investigated on a triangular-lattice ribbon, in which the dynamics of the system is governed by the discrete reduced integrable nonlinear Schrödinger coupled equations [65]. In the same light, the propagation and elastic interactions of discrete solitons have equally been investigated with higher-order coupled Ablowitz-Ladik equation [66]. Consequently, our next study will focus on using the discrete Hirota-Bilinear method to obtain multi-soliton (i.e. two, three, and four soliton solutions) for the DNLS equation with NNN interactions. This will give us the opportunity to investigate on the effects of helicoidal interactions on the propagation and elastic collision of discrete solitons in DNA structures.