Solitary Waves in a Discrete Nonlinear Dirac equation
SSolitary Waves in a Discrete Nonlinear Dirac equation
Jes´us Cuevas–Maraver
Grupo de F´ısica No Lineal, Departamento de F´ısica Aplicada I,Universidad de Sevilla. Escuela Polit´ecnica Superior, C/ Virgen de ´Africa, 7, 41011-Sevilla, SpainInstituto de Matem´aticas de la Universidad de Sevilla (IMUS). Edificio Celestino Mutis. Avda. Reina Mercedes s/n, 41012-Sevilla, Spain
Panayotis G. Kevrekidis
Department of Mathematics and Statistics, University of Massachusetts, Amherst, MA 01003-9305, USA
Avadh Saxena
Center for Nonlinear Studies and Theoretical Division,Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA
In the present work, we introduce a discrete formulation of the nonlinear Dirac equation in the form of adiscretization of the Gross-Neveu model. The motivation for this discrete model proposal is both computational(near the continuum limit) and theoretical (using the understanding of the anti-continuum limit of vanishingcoupling). Numerous unexpected features are identified including a staggered solitary pattern emerging froma single site excitation, as well as two- and three-site excitations playing a role analogous to one- and two-site, respectively, excitations of the discrete nonlinear Schr¨odinger analogue of the model. Stability exchangesbetween the two- and three-site states are identified, as well as instabilities that appear to be persistent overthe coupling strength (cid:15) , for a subcritical value of the propagation constant Λ . Variations of the propagationconstant, coupling parameter and nonlinearity exponent are all examined in terms of their existence and stabilityimplications and long dynamical simulations are used to unravel the evolutionary phenomenology of the system(when unstable). I. INTRODUCTION
Nonlinear dispersive waves in physical systems are often described by the nonlinear Schr¨odinger equation (NLSE), which isboth mathematically and physically studied in a broad range of settings including atomic physics [1], nonlinear optics [2] andmathematical physics [3, 4]. Both the continuum and the discrete [5, 6] installment of the equation have been analyzed in detail.A principal focus of the relevant properties, aside from issues of self-focusing and wave collapse [3] has been the study of theexistence, stability and dynamics of solitary waves in this model, both in lower-dimensional settings (such as one-dimensionalsolitons and multi-solitons) and in higher dimensional settings (vortices, vortex rings, and related structures).However, recent years have seen a gradual increase of interest in the study of near-relativistic settings, where a suitablegeneralization/extension of the NLSE is the so-called nonlinear Dirac equation (NLDE) [7]. In fact, a modified form of the NLSE(with additional terms) is a special case limit of the NLDE at the low-energy limit as has been demonstrated in [8]. Differentrealizations of the NLDE have been proposed in the realm of high energy physics, including the so-called massive Gross-Neveumodel [9], as well as the massive Thirring model [10]. Importantly, the equation has seen a significant volume of studies froma more mathematical perspective. Various aspects have been examined in this context, including the spectral stability and thepotential emergence of point spectrum eigenvalues with nonzero real part (which has been shown to be impossible to happenbeyond the so-called embedded thresholds) [11], the orbital and asymptotic stability under a series of relevant assumptions [12],the nonlinear Schr¨odinger (non-relativistic) limit and its instability for nonlinearities beyond a critical exponent [13], as well asclassical (Vakhitov-Kolokolov) and more suitable to this setting (energy based) criteria [14] for the linear stability of solitarywaves in the NLDE. A series of more computationally/physically oriented studies both in the context of the stability/dynamicsof the NLDE solitary waves [15, 16] (again, in principle for arbitrary nonlinearity powers) and in that of these structures in thepresence of external fields [17] have also recently appeared.It would be relevant to mention, at least in passing, one more framework where Dirac-type equations have received significantattention in recent years in the context of atomic Bose-Einstein condensates. This is, in particular, in the context of artificialgauge fields more broadly, and more specifically spin-orbit coupled Bose-Einstein condensates [18]. There, admittedly, thesetup is somewhat different, as both the Dirac type operator and the Schr¨odinger one co-exist, but the compensating value isthat such settings have already been realized experimentally [19–23]. Moreover, a wide range of coherent structures has beenalready proposed in them including vortices [24, 25], Skyrmions [26], Dirac monopoles [27] and dark solitons [28, 29], as wellas self-trapped states [30], bright solitons [31] and gap-solitons [32].In the present work, we will take a somewhat different path from the above works. In particular, we will consider a relativelystandard example of the NLDE (namely the so-called massive Gross-Neveu model), but motivated by the significant level ofunderstanding and analytical tractability afforded by discrete settings [5, 6], we will instead consider a spatially discrete form ofthe NLDE. A significant part of our motivation for this consideration (and for the particular form of the selected discretization) a r X i v : . [ n li n . PS ] A ug is due to (a) the possibility to deploy the technology of the so-called anti-continuum (AC) limit of MacKay-Aubry [33], inorder to appreciate the stability properties near the limit of uncoupled adjacent sites and (b) the feature that in the continuumlimit of, in principle, infinite coupling, our conclusions are expected to connect to what is known for the corresponding PDEmodels that have been explored in the literature. Admittedly, the discretization that is selected herein is, arguably, not themost natural possible one (in that we utilize next-nearest neighbors in order to discretize the first derivative terms by centereddifferences). Nevertheless, it is identified that it is the most suitable one for the present setting type of stencil and discretesolitary waves are systematically obtained from the AC limit. Moreover, a very recent development worth noting is that spin-orbit Bose-Einstein condensates have recently been considered in the realm of an optical lattice [34], which is often thought (inthe so-called superfluid regime) as being tantamount to a discretization of the original continuum problem, through a suitableWannier function reduction [35]. This suggests that considering discrete variants of Dirac models may be a natural step for nearfuture considerations.While many of our findings are somewhat reminiscent of the corresponding discrete nonlinear Schr¨odinger (DNLS) equationones [5, 6], numerous others are rather unique to the Dirac equation. The single site solution is found to lead to a ratherunexpected waveform which we explain and illustrate to effectively (that is, in its envelope) approach in the continuum limit thesolution of a different homoclinic state problem that will be explicitly discussed below. On the other hand, it is the two-site andthree-site solutions that lead to a continuation all the way to the continuum limit of the Gross-Neveu solitary wave. However,contrary to what is the case for the DNLS, the two-site solution turns out to be stable close to the AC limit, while the three-sitesolution is the unstable one close to that limit. A count of the relevant eigenvalues near the AC limit is systematically given forthese different cases. Subsequently a near-alternation of stability is observed between these two modes (the site- and inter-site-centered ones) that is somewhat reminiscent of the phenomenology identified in the saturable DNLS model [36, 37]. This isexplored systematically, as is the feature of both of these solutions in producing a complex quartet of modes in a suitable bandof the continuous spectrum. This oscillatory instability and its dynamical by-products are traced as a function of the propagationconstant Λ and of the inter-site coupling strength (cid:15) . The dynamical manifestation of the instabilities within the discrete model isshown to lead to different possible features, including the potential mobility of the solitary waves or their splitting into multiplesolitary waves of lower amplitude (and potentially of a different type).Our presentation is structured as follows. In section II, we present an overview of our discrete model and its basic properties.In section III, we examine the different solutions in the vicinity of the AC limit. In section IV we examine the same solutionsfor large (cid:15) , i.e., in the vicinity of the corresponding continuum limit. Finally, in section V, we explore the dynamical instabilitymanifestations of the different solutions. Section VI summarizes our findings and presents our conclusions. II. MODEL AND THEORETICAL SETUP
The NLDE model that we will consider will be the massive Gross-Neveu model with scalar-scalar interactions and a generalpower-law nonlinearity. This is motivated by recent corresponding continuum model explorations both at the level of mathemat-ical analysis [13] and at that of numerical computations [15, 16]. The discrete version of the equation introduced herein will bebased on a centered difference approximation of the first derivative in the form: i ˙ U n = (cid:15) ∇ V n − g ( | U n | − | V n | ) k U n + mU n ,i ˙ V n = − (cid:15) ∇ U n + g ( | U n | − | V n | ) k V n − mV n , (1)with U n and V n being the components of the spinor Ψ n ≡ ( U n , V n ) and ∇ Ψ n ≡ (Ψ n +1 − Ψ n − ) being the discrete gradient,with a centered difference scheme, as indicated above. The connection to the corresponding continuum limit can be assignedby selecting (cid:15) = 1 / (2 h ) with h being the lattice spacing (discretization parameter). It should also be noted in passing that weattempted to discretize by a forward difference scheme, with considerably less promising results. Given also that the centereddifference scheme is a higher order discrete approximation to the corresponding continuum limit, we therefore will only presentresults by means of the centered difference discretization scheme in what follows.Our main focus hereafter will be on stationary solutions and their stability. Such solutions can be found by using U n ( t ) =exp( − i Λ t ) u n , V n ( t ) = exp( − i Λ t ) v n , and satisfy the coupled algebraic equations: (cid:15) ∇ v n − g ( | u n | − | v n | ) k u n + ( m − Λ) u n = 0 ,(cid:15) ∇ u n − g ( | u n | − | v n | ) k v n + ( m + Λ) v n = 0 . (2)Analogously to its continuum counterpart, the dynamical system of Eq. (1) presents a number of conserved quantities, suchas the charge (squared (cid:96) norm): Q = (cid:88) n ρ n , ρ n = | U n | + | V n | , (3)with ρ n being the charge density, and the Hamiltonian: H = 12 (cid:88) n (cid:20) ( U ∗ n ∇ V n − V ∗ n ∇ U n ) − gk + 1 ( | U n | − | V n | ) k +1 + m ( | U n | − | V n | ) (cid:21) . (4)The dynamical equations (1) can be derived from the Hamiltonian (4) by means of the Hamilton’s equations: i ˙ U n = δHδU ∗ n , i ˙ V n = δHδV ∗ n . (5)Once stationary solutions of the algebraic system of Eqs. (2) are calculated (by e.g. fixed points methods), their linearstability is considered by means of a Bogoliubov-de Gennes linearized stability analysis. More specifically, considering smallperturbations [of order O( δ ) , with < δ (cid:28) ] of the stationary solutions, we substitute the ansatz U n ( t ) = e − i Λ t (cid:104) u n, + δ ( a n e iωt + c ∗ n e − iω ∗ t ) (cid:105) , V n ( t ) = e − i Λ t (cid:104) v n, + δ ( b n e iωt + d ∗ n e − iω ∗ t ) (cid:105) (6)into Eqs. (1), and then solve the ensuing [to O ( δ ) ] eigenvalue problem: ω a n b n c n d n = M a n b n c n d n , (7)with M being M = Λ + J n + L n ( u, u ∗ ) −∇ − L n ( u, v ∗ ) L n ( u, u ) − L n ( u, v ) ∇ − L n ( u ∗ , v ) Λ − J n + L n ( v, v ∗ ) − L n ( u, v ) L n ( v, v ) − L n ( u ∗ , u ∗ ) L n ( u ∗ , v ∗ ) Λ − J n − L n ( u, u ∗ ) ∇ + L n ( u ∗ , v ) L n ( u ∗ , v ∗ ) − L n ( v ∗ , v ∗ ) −∇ + L n ( u, v ∗ ) Λ + J n − L n ( v, v ∗ ) (8)for the eigenvalue ω and associated eigenvector { ( a n , b n , c n , d n ) T } . Here, L n ( x, y ) is a function defined as: L n ( x, y ) = kχ k − n x n, y n, , (9)with J n being: J n ≡ gχ kn − m , (10)and χ n ≡ | u n, | − | v n, | . (11)The dispersion relation of the linear excitations corresponds to the continuous spectrum that will be identified in the lin-earization around the trivial u n = v n = 0 ∀ n solution. This relation can be identified by decomposing the perturbations as { a n , b n , c n , d n } = { A, B, C, D } exp( iqn ) in Eqs. (2) and deriving the resulting condition: ω ( q ) = ± Λ ± (cid:113) m + 4 (cid:15) sin q. (12)Consequently, there are two sets of bands in the essential spectrum. The embedded part given by | ω | ∈ [ − Λ + m, − Λ + √ m + 4 (cid:15) ] , and the non-embedded part, | ω | ∈ [Λ + m, Λ + √ m + 4 (cid:15) ] . In what follows, for concreteness we will set m = 1 and vary Λ , as well as (cid:15) , as our relevant parameters. III. RESULTS FROM THE AC LIMIT (SMALL COUPLING REGIME)
In this section, we consider the existence, stability and dynamics of discrete solitons from the AC to the continuum limit. Inthe AC limit, the soliton that can be continued up to the continuum is a three-site soliton, given by v n = 0 ∀ n and u − = u = u = (1 − Λ) / (2 k ) , u n = 0 ∀ | n | ≥ . FIG. 1: Spectrum of the stability matrix (8) for discrete NLDE 3-site solitons with
Λ = 0 . (top) and Λ = 0 . (bottom). The imaginary(left) and real (right) parts of the corresponding eigenfrequencies are shown as a function of the coupling strength (cid:15) . Only the positive realand imaginary parts of the eigenfrequencies are shown. The size of the system is N = 801 . In the top right, the inset is a magnification ofthe relevant Im ( ω ) shown in the figure but at a different scale. The inset in the bottom right panel shows the oscillations of the growth rate fordifferent system sizes when Λ = 0 . . Notice that oscillation amplitude decreases rapidly as the number of lattice nodes increases. Let us explain below the general behavior for Λ > / . Outside this range, the solitary waves are always unstable and hencewe do not consider them further here.It is easy to see from (8) that in the AC limit and for any value of k , this three-site solution possesses 3 pairs of modes at ω = 0 ,3 pairs at ω = ± , ( N − pairs at ω = ± (1 + Λ) and ( N − pairs at ω = ± (1 − Λ) . When the coupling is switched on(see Fig. 1), the wave becomes exponentially unstable because of one among the 3 pairs at ω = 0 that detaches from the originyielding an imaginary eigenfrequency pair in a similar way as occurs e.g. for the two-site structure in the DNLS equation [5].The other two vanishing eigenfrequency pairs remain at the origin. In addition, the eigenmodes at ω = ± detach into threepairs that will subsequently collide with the embedded and essential parts of the spectrum; let us denote those modes as A, B, C(from upper to lower real part of the eigenfrequency). Mode C remains exactly at ω = ± for every coupling. The real partof the eigenfrequency of mode A rapidly increases entering the embedded spectrum at the point where the imaginary part of theeigenmode responsible for the exponential instability reaches its maximum. The exponential instability mentioned previouslydisappears close to (but not at) the point where mode C enters the essential spectrum. However, when the coupling increases,the exponential instability appears again with a similar (non-monotonic) behavior as the previous one, except for the presence ofsmaller growth rates and of a slower decrease in the growth rate (past the point of the maximal growth rate). The most complexparametric dependence is the one experienced by mode B. The latter enters the essential spectrum for a value of (cid:15) higher thanthat for which mode C enters therein. Then, the system becomes oscillatorily unstable and undergoes a Hopf bifurcation [in thecase of finite systems, due to the quantization of the continuous spectrum, this translates into a series of instability bubbles; fora similar scenario in the DNLS see e.g. [38]]. As a consequence, there are many oscillations in the imaginary part of mode Bwhen the coupling is high; the amplitude of those oscillations decreases when the system size increases, as shown in the inset ofbottom right panel of Fig. 1]. When the frequency increases (say Λ (cid:38) . ) the imaginary part of mode B does not asymptoteto a nearly constant value as the coupling strength increases, but, on the contrary, a series of bubbles appears manifesting asoscillations around 0 (see top panels of Fig. 1). The persistence of those oscillations in the continuum limit will be consideredin Section IV (see also Fig. 10 therein).A more complete scenario of the linearization eigenfrequencies is presented in Fig. 2, where the largest imaginary part ofeigenfrequencies with zero and non-zero real part (i.e. responsible for exponential and Hopf bifurcations, respectively), withrespect to (cid:15) ≤ . and / ≤ Λ < is presented. A cut-off for growth rates smaller than − for Hopf bifurcations and − for exponential bifurcations has been introduced. It can be observed that the instability bubbles emerge in the right panel FIG. 2: Logarithm of the largest imaginary part of eigenvalues with zero (left) and non-zero (right) real part with N = 801 . Blank areascorrespond to stable solitons. for Λ (cid:38) . and (cid:15) (cid:38) . . In addition, it is observed that exponential instabilities emerge in several lobes, which suggestsa cascading mechanism of destabilizations and restabilizations that we will return to below, upon examination of the two-sitesolitary wave.A complementary scenario is experienced by the two-site solitary wave, given in the AC limit by v n = 0 ∀ n and u = u =(1 − Λ) / (2 k ) , u n = 0 elsewhere. At this limit, the 2-site structure possesses 2 pairs of modes at ω = 0 , 2 pairs at ω = ± , ( N − pairs at ω = ± (1 + Λ) and ( N − pairs at ω = ± (1 − Λ) . When the coupling is switched on (see Fig. 3), the structureremains stable because of the persistence of both pairs at ω = 0 . Mode A from | ω | = 2Λ does not exist for this case; on theother hand, the oscillatory instabilities caused by mode B also exist for the 2-site case. When increasing the coupling, the solitonexperiences an exponential bifurcation and becomes unstable, contrary to the 3-site soliton (notice that in typical Klein-Gordonand –e.g. saturable– DNLS settings, such stability exchanges take place between 2-site and 1-site breathers or solitons). Here,there are exponential stability exchanges between 2-site and 3-site solitons, although the bifurcations of the two families ofsolutions do not perfectly coincide (nevertheless, in a number of such exchanges, the corresponding stabilization/destabilizationthresholds are fairly proximal). This scenario is summarized in Fig. 4. We should note in passing that these near-exchangesof stability suggest a scenario similar to the ones occurring e.g. in the saturable or cubic-quintic DNLS model where the near-exchange of stability of the 1- and 2-site solitary waves (in that case) is mediated through a series of pitchfork and reversepitchfork bifurcations of asymmetric solution branches [39, 40]. However, we will not pursue the relevant narrow branches ofasymmetric solutions herein.There is an interesting kind of solution that also exists from the AC limit and can be continued to the continuum limit, namelythe one-site soliton. This has the following property which is, in fact, preserved upon continuation for any value of the coupling(see Fig. 5): u n = 0 for odd n and v n = 0 for even n ; however, the charge density of the soliton is qualitatively different fromthat of the three-site solitons. In the AC limit (cid:15) = 0 , u = (1 − Λ) / (2 k ) , and u n = 0 for the rest of sites (with v n = 0 ∀ n ).The form of this solution can be identified as we approach the continuum limit as u n = 0 for odd n and v n = 0 for even n ,by transforming the discrete NLDE equation (1) into the new set of equations: (cid:15) ( v n +1 − v n − ) − gu k +1 n + ( m − Λ) u n = 0 , for even n(cid:15) ( u n +2 − u n ) − ( − k gv k +1 n +1 + ( m + Λ) v n +1 = 0 , for even n (13)which possesses homoclinic solutions in the continuum limit (see Section IV).The spectrum of the one-site solitons at (cid:15) = 0 consists of a single pair of eigenvalues at ω = 0 and another single pair at ω = ± ; apart from these, there are N − pairs at ω = ± (1 + Λ) and ω = ± (1 − Λ) . When the coupling is switched on,as there is only a single pair of eigenmodes at ω = 0 , the soliton does not experience exponential bifurcations; in addition, thenon-existence of mode B prevents the existence of harmful Hopf bifurcations arising in 3-site and 2-site solitons (see Fig. 6).The only observed instability is an exponential one arising for a finite value of coupling and caused by a mode that bifurcatesfrom the essential spectrum as the coupling strength increases; the growth rate which has a non-monotonic dependence on thecoupling and tends asymptotically to zero when reaching the continuum limit, and its maximum value decreasing with Λ (forfixed coupling) are shown in Fig. 7 for more details. Similar to the 3-site structures, there is a complementary family of solitonsconsisting of 2-site structures with a hole in between, characterized by u = u . FIG. 3: Same as Fig.1 but for the discrete NLDE 2-site solitons with
Λ = 0 . (top) and Λ = 0 . (bottom).FIG. 4: (cid:15) vs Λ plane where different unstable regimes for 3-site and 2-site solitons are displayed. 3-site solitons are unstable inside the fullline regions, whereas unstable 2-site solitons are inside dashed lines. Above the dashed-dotted lines, both solutions are oscillatorily unstable. IV. RESULTS NEAR THE CONTINUUM LIMIT (LARGE COUPLING REGIME)
In this section, on the one hand, we will connect the findings of our model with some previous results about the stability ofthe continuous NLDE. On the other hand, we will perform the Bogoliubov-de Gennes (BdG) spectral stability analysis of thediscrete NLDE for a large coupling such as (cid:15) = 5 , which corresponds to a spatial discretization parameter h = 0 . .Previous results from Comech (see Refs. [11, 41–43]) show that close to the non-relativistic limit ( Λ (cid:46) ), the Vakhitov–Kolokolov criterion should hold. Based on it [42], it is concluded that no unstable eigefrequency should emerge from ω = 0 closeto this limit for k = 1 or k = 2 (contrary to the k ≥ , k ∈ N case where a pair of eigenfrequencies with a nonzero imaginarypart and a zero real part are present) and, consequently, no exponential instability should exist in that limit. Additionally, forany k , the existence of an eigenfrequency | ω | = 2Λ is also predicted. This mode enters the linear mode band at Λ = 1 / (whenincreasing Λ ).The work of [17], based on the so-called Bogolubsky criterion, as well as that of [11] suggest that solitary waves are alwaysunstable for Λ < Λ c . [It is worth noting here that neither of the two criteria mentioned above is able to give a necessary condition FIG. 5: (Left) The two-component profiles for a 1-site soliton with
Λ = 0 . at (cid:15) = 5 . (Right) Charge density ρ n = | u n | + | v n | for thesolitary wave at the left (blue line) and a 3-site soliton with the same parameters (red line). It is clear that the former does not asymptote to thelatter, but rather to a different envelope that will be revealed in section IV below.FIG. 6: Same as Fig.1 but for discrete NLDE 1-site solitons with Λ = 0 . (top) and Λ = 0 . (bottom). and, consequently, the minimum value for which solitons are stable must be determined numerically]. In the cubic case ( k = 1 )it is predicted in [17] that Λ c = 0 . . However, in a recent paper [16], further numerical simulations have suggested thatsolitons may be dynamically stable for Λ ≥ . .The analytical form of the profile of solitons in the continuum limit is given by [15, 17]: u ( x ) = (cid:115) (1 + Λ) cosh ( kβx )1 + Λ cosh(2 kβx ) (cid:20) ( k + 1) β kβx ) (cid:21) / k , v ( x ) = (cid:115) (1 − Λ) sinh ( kβx )1 + Λ cosh(2 kβx ) (cid:20) ( k + 1) β kβx ) (cid:21) / k , (14)with β = √ − Λ . When k = 1 the equations above can be simplified to: u ( x ) = (cid:112) − Λ)[1 − µ tanh ( βx )] cosh( βx ) , v ( x ) = (cid:112) µ (1 − Λ) tanh( βx )[1 − µ tanh ( βx )] cosh( βx ) , (15) FIG. 7: Exponential instability loci (left) for a 1-site soliton are shown in the (cid:15) - Λ plane in the left panel. The right panel shows the maximum(over the considered (cid:15) variations) growth rate at each value of Λ . For Λ > . , the growth rates for the 1-site soliton are smaller than − and cannot be accurately traced because of machine precision.FIG. 8: (Left) Solitary wave profiles and (right) spectral planes for h = 0 . , L = 80 and Λ = 0 . (top) and Λ = 0 . (bottom). Linesin the right figures indicate the limits of the embedded and non-embedded parts of the essential spectrum. Instabilities in the non-embeddedspectrum should vanish when L → ∞ . with µ = (1 − Λ) / (1 + Λ) . As demonstrated in [15], continuous solitons become double-humped for Λ smaller than a criticalvalue for every k . Fig. 8 shows the profile and spectral planes for two different examples of solitons close to the continuum limitwith k = 1 .We show in Fig. 9 the stability eigenvalues for k = 1 in a domain x ∈ [ − L/ , L/ , with L = 80 and a discretization step h = 0 . . Although there are instabilities caused by eigenvalue collisions in the non-embedded spectrum, we have neglectedthem, as they disappear in the limit of h → and L → ∞ . The waves are found to be unstable for small Λ , with a growthrate that decreases when Λ is increased. The source of instabilities is a localized mode (with non-zero imaginary part of itseigenfrequency even when Λ → ) that enters the embedded spectrum at Λ ≈ . . Once inside the linear modes band,this localized mode causes multiple bubbles, yet at Λ ≈ . , it returns to the real eigenfrequency axis and the solitary wavebecomes stable. Nevertheless, this stability is ephemeral, as the soliton becomes unstable again at Λ ≈ . . From this point,there is a succession of instability bubbles, whose amplitude (i.e., the maximal growth rate associated with them) decreases with Λ . Notice also the existence of the eigenvalue with ω = 2Λ , which enters the embedded spectrum at Λ = 1 / . FIG. 9: Spectrum of the stability matrix (8) for solitons with h = 0 . and L = 80 . For the sake of simplicity, only the positive real andimaginary parts of the eigenvalues are shown.FIG. 10: Dependence of the growth rates of the solitary waves shown in Fig. 9 for different domain lengths. As L increases, we progressivelycan discern the envelope of the infinite domain limit. In order to observe the behavior of bubbles when the domain is enlarged, we have included Fig. 10 where the growth rateis plotted for L = 80 , and . It is observed that the number of bubbles increases with L , but their width decreases. Inany case, the envelope of the bubbles tends to zero asymptotically when Λ approaches 1, in a similar way as it was observed fordark solitons in DNLS settings [38]. Unfortunately, the convex nature of the relevant (apparent) envelope curve is inconclusivein connection to the stability aspect. In particular, it is unclear, based on the present computations, whether the curve, as h → , still intersects the axis and no longer features an unstable mode past a critical Λ c , as is the case with our finite h , finite domaincomputations in Fig. 10, The alternative scenario is that the approach to the stable NLS limit of Λ → (a glimpse of which isillustrated in Fig. 10) is merely asymptotic. It would be especially interesting to pursue this intriguing aspect further, pushingthe envelope of the currently available numerical tools.As mentioned in Section III, the 1-site solitons can also exist in the continuum limit. There, by neglecting the irrelevant (inthis setting) inactive odd sites for one of the fields, and the even ones for the other, the envelope of the solitary waves can be seento approach the homoclinic orbits of the following system of ODEs that is found by obtaining the continuum limit of (13): ∂ x u = ( − k gv k +1 − ( m + Λ) v,∂ x v = gu k +1 − ( m − Λ) u . (16)Using phase plane numerical analysis (not shown here), we have confirmed that Eqs. (16) possess a homoclinic orbit g = m = 1 and k = 1 , for a wide range of Λ ’s. We have also confirmed that it is at these very homoclinic orbits that the envelope of ourNLDE 1-site solution tends as the coupling strength is increased.As a final comment regarding the stability analysis, we note that our approach allows to examine not only the variations as afunction of the propagation constant Λ , as well as the coupling strength (cid:15) , but additionally also with respect to the nonlinearityexponent parameter k . Fig. 11 shows some typical examples of this variation for h = 0 . and values of Λ = 0 . (top) and Λ = 0 . (bottom). The parametric variation of k reveals both the Hopf and exponential instabilities of the system. As regardsthe latter, we note that for sufficiently high values of k , an eigenfrequency bifurcating from the continuous spectrum crossesthe spectral plane origin becoming imaginary, in accordance with the expectation that for sufficiently high k a blow-up type0 FIG. 11: Spectrum of the stability matrix (8) for solitons with h = 0 . and L = 80 . Λ = 0 . in upper panels and Λ = 0 . in bottom ones, inthis case as a function of the nonlinearity exponent k . In the former case, exponential instabilities emerge for k < . and k > . , whereasin the latter case, those instabilities take place for k < . and k > . .FIG. 12: k vs Λ plane where the behavior of solitons with respect to exponential instabilities is displayed. Notice that the critical value k = 2 is retrieved at the non-relativistic (NLSE) limit. instability (which for Λ → , i.e., the NLSE limit, should occur for k > ) should emerge. The relevant critical points are k = 3 . and k = 2 . , respectively for the considered values of Λ of . and . . On the other hand, another interestingobservation is that a similar bifurcation to exponential instability appears to emerge in the small positive k i.e., the weaklynonlinear limit. This instability arising for k < . and k < . in the top and bottom, respectively, panel of Fig. 11 is worthexamining further in its own right, possibly through a perturbative calculation.A systematic exploration of the k - Λ plane of the relevant exponential instability is identified in Fig. 12. As expected in thenon-relativistic NLSE limit of Λ → , no exponential instabilities are observed when k < , whereas for k > , the solitons areunstable i.e., amenable to collapse. We can see that as Λ decreases from that limit, the corresponding critical k for the instabilitymonotonically increases.1 FIG. 13: Evolution of the soliton with
Λ = 0 . and (cid:15) = 0 . . Top panel shows the evolution of the charge density ρ n . Right panels depict thefields at t = 0 (dashed line) and at t = 1000 (solid line). Bottom left panel displays the spectral plane of the solitary wave. V. DYNAMICAL EVOLUTION OF INSTABILITIES
We have also briefly analyzed the dynamics of unstable 3-site solutions in different regimes. Below, we give a numberof selected results in connection to the relevant numerical evolution, although admittedly a systematic classification of thedynamical implications of the different identified instabilities and of the various possible configurations identified herein is aseparate numerical project in its own right.Fig. 13 shows the evolution of a solitary wave with
Λ = 0 . and (cid:15) = 0 . , i.e. inside the lower lobe of exponential instabilities.We observe that the structure emits linear wave “radiation” and subsequently deforms towards a more compact configurationwith fewer high-amplitude excited sites (more specifically one in each component). If a solution within the intermediate lobe istaken (as e.g. that of Fig. 14, where Λ = 0 . and (cid:15) = 0 . ), it is observed that the soliton moves along the lattice. In this regime,both an exponential and an oscillatory instability are present. Generally, for cases of larger coupling, we find that the solutionsare more prone to becoming mobile, upon the manifestation of the dynamical instability.Fig. 15 depicts an oscillatorily unstable wave with Λ = 0 . and (cid:15) = 2 . Interestingly, the latter splits into two daughter-wavesas a result of the oscillatory growth. Once the original structure is split, the charge density at even sites is close to zero (a statesimilar to the 1-site soliton). That is, the offspring in this case belong to the same class of solitary waves as the 1-site solutionexamined above. If an oscillatory unstable soliton is taken within the region of oscillatory instability bubbles (see Fig. 16, where Λ = 0 . and (cid:15) = 2 ), the soliton is put into motion. Once again, this is a relatively common feature of case examples with largevalues of (cid:15) that are more proximal to the continuum limit of the problem. Notice, however, additionally that in the process ofshedding away radiative wavepackets that manifests the dynamical instability and sets the solitary wave in motion, the amplitudeof the structure decreases, which indicates that its effective Λ increases and hence renders it more robust.Lastly, we should note that we have also examined dynamical instabilities of other structures such as 1- and 2-site solitons.These often, too, result in mobile coherent structures, especially for large values of (cid:15) , although some states, such as the staggered1-site wave are less amenable to extensive traveling throughout the lattice, perhaps partly due to their special spatial structure. VI. CONCLUSIONS AND FUTURE CHALLENGES
In the present work we have examined a lattice analogue of the nonlinear Dirac equation. Motivated by the considerablevolume of both mathematical and computational investigations of the continuum limit of the corresponding problem, we havedeveloped a prototypical discretization scheme whose continuum limit is the Gross-Neveu model. This is a model that while2
FIG. 14: Evolution of the soliton with
Λ = 0 . and (cid:15) = 0 . . Top panel shows the evolution of the charge density ρ n . Right panels depict thefields at t = 0 (dashed line) and at t = 5500 (solid line). Bottom left panel displays the spectral plane of the solitary wave.FIG. 15: Evolution of the soliton with Λ = 0 . and (cid:15) = 2 . Top panel shows the evolution of the charge density ρ n . Right panels depict thefields at t = 0 (dashed line) and at t = 4000 (solid line). Bottom left panel displays the spectral plane of the solitary wave. FIG. 16: Evolution of the soliton with
Λ = 0 . and (cid:15) = 2 . Top panel shows the evolution of the charge density ρ n . Right panels depict thefields at t = 0 (dashed line) and at t = 8000 (solid line). Bottom left panel displays the spectral plane of the solitary wave. it does not presently possess a straightforward physical realization (e.g. analogous to what is the case for its DNLS analogueand, say, optical waveguide arrays [44]), it nevertheless is of substantial interest in its own right for a number of reasons. Itis useful (and relevant to understand), on the one hand, as a numerical scheme and a computational tool for approximatingthe corresponding continuum limit (in regimes of large coupling strength (cid:15) ). On the other hand, its analytical tractability in thevicinity of the anti-continuum limit of uncoupled sites makes it a useful starting point for the exploration of the spectral propertiesof solitary waves. In the AC limit, there is a complete control over these spectral properties and corresponding eigenvalues, andit then remains to appreciate the continuation of these over the coupling strength (cid:15) , in order to understand both the features ofthe discrete model and those of its continuum limit. Moreover, the physical realization of quasi-discrete systems possessingDirac-like dynamics such as spin-orbit Bose-Einstein condensates in the presence of optical lattices very recently [34], seems tostrongly suggest the potential experimentally-relevant realization of models of this class in the near future.In light of the above motivations, here we have shown a multitude of unexpected properties that merit further studies notonly from a computational but also importantly from a rigorously mathematical point of view. In particular, we showed thata single site excitation does not continue, as might be expected, to a continuum solitary wave of the Gross-Neveu model.Instead, it forms a remarkable staggered structure that approaches in the limit of (cid:15) large (while being preserved as a state) theenvelope of the homoclinic orbit of a different dynamical model. This appears to be a fairly robust structure in its parametricdependence over (cid:15) and Λ . On the other hand, the two- and three-site initial excitations play the role, respectively, of the one-and two-site excitations of the DNLS. Yet, here a situation more akin to the saturable analogue of the DNLS occurs [36, 37],whereby exchanges of stability between the on-site and inter-site solutions arise. Likely, and analogously to corresponding DNLScubic-quintic or saturable settings, these exchanges are mediated by pitchfork bifurcations (and reverse pitchforks) generatingasymmetric waveforms, a topic potentially worthy of further investigation in the future. Additionally, these states appear topossess quartet of eigenfrequencies chiefly responsible for their instability. While this instability appears to reach an asymptoticgrowth rate (over (cid:15) variations) for values of Λ below a critical one, it is an open problem whether indeed this instability isexpected to be present in the continuum limit of the problem. In that connection, it is relevant to point out that we have observedthe manifestation of the instability to potentially lead to traveling and mobility of the structure, while in other cases, we haveobserved it to lead to a fragmentation of the solitary wave into the staggered structures, a feature which would not be “accessible”in the continuum limit. It should also be pointed out that despite our computation for different domain sizes as a function of Λ forlarge (approaching the continuum) values of (cid:15) , the concavity of the relevant eigenvalue dependence precludes a straightforwarddetermination of the associated critical value of Λ . It can be safely inferred that the instability (at least for sub-critical exponents k < /n ) is absent in the nonrelativistic Schr¨odinger limit. Nevertheless, whether this occurs asymptotically as Λ → or at a4finite Λ c (the latter being observed in the case of our finite -but large- coupling and domain size) remains yet another importantopen question.Naturally, the present investigation, as a primary one of its kind, raises a considerable volume of additional questions mer-iting future examination both at the discrete and at the continuum limit. In particular, a key issue is how the discrete modelasymptotes to the actual corresponding continuum. It is especially important to understand how the spectral properties may bemodified in the limit. Another very interesting avenue of research would be to develop and utilize the solvability conditions thatwere especially handy in the DNLS case to understand the unusual existence and stability features in the corresponding Diraccase. Understanding also better the role (especially in the dynamics) of the unusual staggered structure would be especiallyrelevant. Other themes, such as a classification of the dynamical instability scenaria for different states or the identification ofthe exponential instability in the near-linear limit of small k have also emerged. Beyond the realm of “single pulses” focusedupon herein, the cases of multi-pulses, pulse interactions and related themes are entirely open, to the best of our knowledge, notonly in the discrete case but largely also in the continuum one. Finally, all these investigations could naturally be generalized tohigher dimensions, where also vortical and related structures could potentially arise [5, 6]. Some of these topics are currentlyunder investigation and will be reported in future publications. Acknowledgements
This work was supported in part by the U.S. Department of Energy (A.S.). P.G.K. acknowledges support from the NationalScience Foundation under grants CMMI-1000337, DMS-1312856, from FP7-People under grant IRSES-606096 from the Bina-tional (US-Israel) Science Foundation through grant 2010239, and from the US-AFOSR under grant FA9550-12-10332. We areindebted to Faustino Palmero for technical assistance with some parts of the manuscript. [1] P.G. Kevrekidis, D.J. Frantzeskakis, and R. Carretero-Gonz´alez,
Emergent Nonlinear Phenomena in Bose-Einstein Condensates ,Springer-Verlag, Berlin, 2008.[2] Yu.S. Kivshar and G.P. Agrawal,
Optical solitons: from fibers to photonic crystals , Academic Press (San Diego, 2003).[3] C. Sulem and P.L. Sulem,
The Nonlinear Schr¨odinger Equation , Springer-Verlag (New York, 1999).[4] M.J. Ablowitz, B. Prinari and A.D. Trubatch,
Discrete and Continuous Nonlinear Schr¨odinger Systems , Cambridge University Press(Cambridge, 2004).[5] P. G. Kevrekidis,
The Discrete Nonlinear Schr¨odinger Equation: Mathematical Analysis, Numerical Computations and Physical Per-spectives , Springer-Verlag (Heidelberg, 2009).[6] D.E. Pelinovsky,
Localization in Periodic Potentials: From Schr¨odinger Operators to the Gross-Pitaevskii Equation , Cambridge Univer-sity Press (Cambridge, 2011).[7] S.Y. Lee, T. K. Kuo, and A Gavrielides, Phys. Rev. D , 2249 (1975).[8] F.M. Toyama, Y. Hosono, B. Ilyas, Y. Nogami, J. Phys. A , 3139 (1994).[9] D. J. Gross and A. Neveu, Phys. Rev. D , 3235 (1974).[10] W. Thirring, Annals Phys. , 91 (1958).[11] N. Boussa¨ıd and A. Comech. On spectral stability of the nonlinear Dirac equation . Preprint. arXiv:1211.3336 [math.AP].[12] N. Boussa¨ıd, S. Cuccagna, Comm. PDE , 1001 (2012).[13] A. Comech, M. Guan, S. Gustafson, arXiv:1209.1146.[14] A. Comech, G. Berkolaiko, A. Sukhtayev, arXiv:1306.5150.[15] F. Cooper, A. Khare, B. Mihaila and A. Saxena. Phys. Rev. E , 036604 (2010).[16] S. Shao, N.R. Quintero, F.G. Mertens, F. Cooper, A. Khare, A. Saxena, arXiv:1405.5547.[17] F.G. Mertens, N.R. Quintero, F. Cooper, A. Khare and A. Saxena. Phys. Rev. E , 046602 (2012).[18] J. Dalibard, F. Gerbier, G. Juzeli¨unas, and P. ¨Ohberg, Rev. Mod. Phys. , 1523 (2011).[19] Y.-J. Lin, R. L. Compton, K. Jimenez-Garcia, J. V. Porto, and I. B. Spielman, Nature , 628 (2009).[20] Y.-J. Lin, K. Jimenez-Garcia, and I. B. Spielman, Nature, , 83 (2011).[21] C. Qu, C. Hamner, M. Gong, C. Zhang, and P. Engels, Phys. Rev. A , 021604(R) (2013);[22] L.J. LeBlanc, M. C. Beeler, K Jimenez-Garcia, A. R. Perry, S. Sugawa, R.A. Williams and I. B. Spielman, New J. Phys. , 073011(2013).[23] Jin-Yi Zhang, Si-Cong Ji, Zhu Chen, Long Zhang, Zhi-Dong Du, Bo Yan, Ge-Sheng Pan, Bo Zhao, You-Jin Deng, Hui Zhai, Shuai Chen,and Jian-Wei Pan Phys. Rev. Lett. , 115301 (2012).[24] X.-Q. Xu and J. H. Han, Phys. Rev. Lett. , 200401 (2011).[25] J. Radi´c, T. A. Sedrakyan, I. B. Spielman, and V. Galitski, Phys. Rev. A , 063604 (2011); B. Ramachandhran, B. Opanchuk, X-J. Liu,H Pu, P. D. Drummond, and H. Hu, Phys. Rev. A , 023606 (2012).[26] T. Kawakami, T. Mizushima, M. Nitta, and K. Machida, Phys. Rev. Lett. , 015301 (2012).[27] G. J. Conduit, Phys. Rev. A , 021605(R) (2012).[28] O. Fialko J. Brand, and U. Z¨ulicke, Phys. Rev. A , 051605(R) (2012) . [29] V. Achilleos, J. Stockhofe, P. G. Kevrekidis, D. J. Frantzeskakis, and P. Schmelcher, EPL , 20002 (2013)[30] M. Merkl, A. Jacob, F. E. Zimmer, P. ¨Ohberg, and L. Santos, Phys. Rev. Lett. , 073603 (2010).[31] V. Achilleos, D. J. Frantzeskakis, P. G. Kevrekidis, and D. E. Pelinovsky, Phys. Rev. Lett. , 264101 (2013).[32] Y. V. Kartashov, V. V. Konotop, and F. Kh. Abdullaev, Phys. Rev. Lett. , 060402 (2013).[33] R.S. MacKay and S. Aubry, Nonlinearity (1994) 1623-1643.[34] C. Hamner, Y. Zhang, M.A. Khamehchi, M.J. Davis, P. Engels, arXiv:1405.4048.[35] G.L. Alfimov, P.G. Kevrekidis, V.V. Konotop and M. Salerno, Phys. Rev. E , 046608, 5 pages (2002).[36] L. Hadzievski, A. Maluckov, M. Stepi´c and D. Kip, Phys. Rev. Lett. , 033901 (2004).[37] T. R. O. Melvin, A. R. Champneys, P. G. Kevrekidis, and J. Cuevas, Phys. Rev. Lett. , 124101 (2006).[38] M. Johansson and Yu.S. Kivshar. Phys. Rev. Lett. (1999) 85.[39] R. Carretero-Gonzlez, J.D. Talley, C. Chong and B.A. Malomed. Physica D
77 (2006).[40] see e.g. R.A. Vicencio and M. Johansson, Phys. Rev. E , 046602 (2006) and references therein.[41] G. Berkolaiko and A. Comech. Math. Model Nat. Phenom. (2012) 13.[42] A. Comech. On the meaning of the Vakhitov-Kolokolov stability criterion for the nonlinear Dirac equation . Preprint. arXiv:1107.1762v2[math.AP].[43] A. Comech.
Linear instability of nonlinear Dirac equation in 1D with higher order nonlinearity . Preprint. arXiv:1203.3859v2 [math.AP].[44] F. Lederer, G.I. Stegeman, D.N. Christodoulides, G. Assanto, M. Segev and Y. Silberberg, Phys. Rep.463