Solitary waves in a two-dimensional nonlinear Dirac equation: from discrete to continuum

Solitary waves in a two-dimensional nonlinear Dirac equation: from discrete to continuum J. 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, Spain Instituto de Matem´ aticas de la Universidad de Sevilla (IMUS). Edificio Celestino Mutis. Avda. Reina Mercedes s/n, 41012-Sevilla, Spain P. G. Kevrekidis Department of Mathematics and Statistics, University of Massachusetts Amherst, Amherst, MA 01003-4515, USA A.B. Aceves Department of Mathematics, Southern Methodist University, Dallas, Texas 75275, 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 explore a nonlinear Dirac equation motivated as the continuum limit of a binary waveguide array model. We approach the problem both from a near-continuum perspective as well as from a highly discrete one. Starting from the former, we see that the continuum Dirac solitons can be continued for all values of the discretization (coupling) parameter, down to the uncoupled (so-called anti-continuum) limit where they result in a 9-site configuration. We also consider configurations with 1- or 2-sites at the anti-continuum limit and continue them to large couplings, finding that they also persist. For all the obtained solutions, we examine not only the existence, but also the spectral stability through a linearization analysis and finally consider prototypical examples of the dynamics for a selected number of cases for which the solutions are found to be unstable. I. INTRODUCTION Optical waveguide arrays [1,2] constitute one of the settings that have led to numerous recent experimental and theoretical developments as regards the analysis of wave phenomena in Hamiltonian lattices. Both in this and in the related context of photorefractive crystals, features such as discrete diffraction [3] and its management [4], Talbot revivals [5], PT -symmetry and its breaking [6], as well as discrete solitons [3,7] and vortices [8,9] were not only theoretically predicted but also experimentally observed. Variants of the theme of optical waveguide arrays have involved multi-component models bearing multiple polarizations [10,11], waveguides featuring quadratic (so-called χ2) nonlinearities [12,13], and the examination of dark-solitonic states [14,15]. Another theme where extensive related studies have been conducted is the atomic physics realm of Bose-Einstein condensates (BECs) in optical lattices [16,17]. Recently, research on binary waveguide arrays has been gaining momentum [18–22]. Part of the reason for this is that under suitable limiting conditions, this system can lead to Dirac-like nonlinear equations that are of increasing prominence and wide relevance in diverse physical contexts. These include, among others, bosonic evolution in honeycomb lattices [23,24] and a growing class of atomically thin two-dimensional (2D) Dirac materials [25] such as graphene, silicene, germanene and transition metal dichalcogenides [26]. They also arise when studying light propagation in honeycomb photorefractive lattices (the so-called photonic graphene) [27–29]. These Dirac settings have been argued to present fundamental differences, e.g., with respect to their more well known, non-relativistic limits of the nonlinear Schr¨ odinger equation, such as the absence of collapse for an extended interval of frequencies in two spatial dimensions; see, e.g., the recent work of [30] and the discussion therein. Our aim in the present work is to revisit the context of binary waveguide arrays, motivated by the realistic models studied in [18,19] for two-dimensional geometries. In an earlier work, we considered the phenomenology of a 1D discrete nonlinear Dirac equation [31]. However, this was done for a discretization of the less straightforwardly applicable (from a physical perspective) Soler model. Here, we turn our attention to the more realistic setting of the waveguide arrays with onsite nonlinearity, aiming to explore existence, spectral stability and dynamics of nonlinear modes. We do this in a complementary way between the (highly) discrete and the continuum limits. On the one hand, we explore the configurations at strong coupling (the continuum soliton) and subsequently reduce the coupling going towards the discrete limit. Here, we eventually find that a configuration bearing nine sites turns out to be the limiting highly discrete analogue of the continuum solitary wave. On the other hand, we also start our search for model solutions from the highly discrete limit of vanishing coupling (the so called anti-continuum limit of [32]) with one- or two-site configurations and continue them into the strong coupling regime. Utilizing a spectral stability analysis, we identify the regimes of lattice coupling as well as of solution frequency for which the relevant waveforms are dynamically stable. When instability is identified, some prototypical examples of the configuration’s unstable evolution are given. Page 1 of 13 AUTHOR SUBMITTED MANUSCRIPT - JPhysA-107720.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript 2 Our exposition will be structured as follows: in Sec. II we provide an overview of the theoretical properties of both the discrete and the continuous model. In Sec. III, we numerically explore the existence, stability and dynamical features of the models, while in Sec. IV we summarize our findings and present our conclusions as well as some interesting directions for future work. II. THEORETICAL SETUP Following the setup of [18,19], the two-dimensional (continuum) Dirac model of relevance to the binary waveguide problem is of the form: i∂tψ1=−(i∂x+∂y)ψ2+ (m−g|ψ1|2)ψ1, i∂tψ2=−(i∂x−∂y)ψ1−(m+g|ψ2|2)ψ2,(1) where (ψ1, ψ2)denotes the spinor field (mode amplitude), while mrepresents the mass associated with the propagation mismatch between the two different types of waveguides. The cubic nonlinearity stems from the Kerr effect and breaks the Lorentz symmetry (contrary, e.g., to what is the case with the Soler model; cf. [30,31]). We consider this model both at the discrete and at the continuum level. A. Continuous model The analysis of the continuum model of Eq. (1) can be performed in radial coordinates, following a procedure similar to the one proposed in [30]. This results in the form: i∂tψ1=−e−iθ (i∂r+∂θ r)ψ2+ (m−g|ψ1|2)ψ1, i∂tψ2=−eiθ (i∂r−∂θ r)ψ1−(m+g|ψ2|2)ψ2.(2) The form of this equation suggests that we look for stationary solutions as ψ(r, t) = exp(−iωt)ϕ(r)with ϕ(r) = [u(r)eiSθ i v(r)ei(S+1)θ],(3) where u(r)and v(r)are real-valued. The value S∈Zcan be cast as the vorticity of the first spinor component (the second component in this formulation has vorticity S+ 1). The equation fulfilled by stationary profiles then only depends on r, casting the problem into a 1D one: (∂r+S+ 1 r)v+ (m−ω−gu2)u= 0, −(∂r−S r)u−(m+ω+gv2)v= 0,(4) with r > 0. As we discuss in the next section, stationary solutions are sought by numerical means. In order to capture the linear stability of the stationary solutions, we introduce the following ansatz into (2): ψ(r, t) =  {u(r) + δ[a1(r)eiqθeλt +b∗ 1(r)e−iqθeλ∗t]}eiSθ i{v(r) + δ[a2(r)eiqθeλt +b∗ 2(r)e−iqθeλ∗t]}ei(S+1)θ e−iωt,(5) and subsequently solve the ensuing [to O(δ)] eigenvalue problem: λ(a1, a2, b1, b2)T=Mq(a1, a2, b1, b2)Twith Mqbeing Mq=i  L1L2 −L∗ 2−L∗ 1  −iq r  σ10 0σ1  ,(6) Page 2 of 13AUTHOR SUBMITTED MANUSCRIPT - JPhysA-107720.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript 3 and L1=  ω−m+ 2gu2−(∂r+S+1 r) (∂r−S r)ω+m+ 2gv2 , L2=  gu20 0gv2 , σ1=  0I I0  .(7) The key observation which facilitates a computation of the spectrum is that the explicit form of Eq. (6) for Mqcontains rand ∂r, but not θ. This allows us to compute the full 2D stability spectrum Mas the union of spectra of the one-dimensional spectral problems: σ(M) = ∪ q∈Z σ(Mq).(8) In what follows, for concreteness we will set m=g= 1 (as this choice can be made by renormalizing the time and the wave function). Now, the only free parameter that will be considered in the continuum limit is the frequency ω. B. Discrete model The discrete version of Eq. (1) will be based on a discretization similar to that considered, e.g., in [19] (cf. the discussion around Eq. (2) therein). From a numerical approximation perspective, this is tantamount to a centered difference discretization of the first derivative in Eq. (1) and leads to: i∂tVn,m =−C[i∇xUn,m +∇yUn,m]+(m−g|Vn,m|2)Vn,m, i∂tUn,m =−C[i∇xVn,m −∇yVn,m]−(m+g|Un,m|2)Un,m,(9) with Un,m and Vn,m (−N/2 + 1 ≤(n, m)≤N/2) being the components of the spinor (amplitude modes) Ψn,m ≡ (Un,m, Vn,m), and ∇xΨn,m ≡(Ψn+1,m −Ψn−1,m)+ν(Ψn+2,m −Ψn−2,m),∇yΨn,m ≡(Ψn,m+1 −Ψn,m−1)+ ν(Ψn,m+2 − Ψn,m−2)being the xand ycomponents of the discrete gradient and νa parameter that accounts for next-nearest-neighbour coupling; ν= 0 corresponds to the case of solely nearest-neighbour coupling when C= 0. The connection to the corresponding continuum limit can be assigned by selecting C= 1/(2h)with hbeing the lattice spacing (discretization parameter). The dynamical system of Eq. (9) presents a number of conserved quantities, such as the charge (squared ℓ2norm): Q=∑ n∑ m ρn,m, ρn,m =|Un,m|2+|Vn,m|2,(10) with ρn,m being the charge density, and the Hamiltonian: H=−1 2∑ n[CV ∗ n,m(i∇xUn,m +∇yUn,m) + CU∗ n,m(i∇xVn,m −∇yVn,m)−g k+ 1(|Un,m|4+|Vn,m|4) + m(|Un,m|2−|Vn,m|2)]. (11) Equations (9) can be derived from the Hamiltonian (11) by means of Hamilton’s equations: i˙ Un,m =δH δU∗ n,m ,i˙ Vn,m =δH δV ∗ n,m .(12) Our main focus hereafter will be on stationary solutions and their stability as well as their dynamics. Such solutions can be found by using Un,m(t) = exp(−iωt)un,m,Vn,m(t) = exp(−iωt)vn,m, when they possess frequency ωand satisfy the coupled algebraic equations: (ω−m+g|vn,m|2)vn,m +C[i∇xun,m +∇yun,m]=0, (ω+m+g|un,m|2)un,m +C[i∇xvn,m −∇yvn,m]=0.(13) Once stationary solutions of the algebraic system of Eqs. (13) are calculated (by, e.g., fixed point methods as discussed below), their linear stability is considered by means of a linearized stability analysis. More specifically, considering small perturbations [of order O(δ), with 0< δ ≪1] of the stationary solutions, we substitute the ansatz Page 3 of 13 AUTHOR SUBMITTED MANUSCRIPT - JPhysA-107720.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript 4 Un,m(t) = e−iωt [un,m +δ(an,meλt +c∗ n,meλ∗t)], Vn,m(t) = e−iωt [vn,m +δ(bn,meλt +d∗ n,meλ∗t)](14) into Eqs. (9), and then solve the ensuing [to O(δ)] eigenvalue problem: λ(an,m, bn,m, cn,m, dn,m)T= M(an,m, bn,m, cn,m, dn,m)Twith Mbeing M=i  L1L2 −L∗ 2−L∗ 1  (15) and L1=  ω−m+ 2g|u|2C(i∇x+∇y) C(i∇x−∇y)ω+m+ 2g|v|2 , L2=  u20 0v2 .(16) The potential existence of an eigenvalue with non-vanishing real part suggests the existence of a dynamical instability. If all the eigenvalues are imaginary, then the solution is spectrally (neutrally) stable. As in the continuum limit, we will set m=g= 1 and vary ωas well as Cas our relevant parameters in order to characterize the behavior of the solution and the variation of its stability properties. III. NUMERICAL RESULTS A. Continuous model We start our exposition by showing the numerical results regarding fundamental solutions (S= 0 solitary waves, for which we nevertheless note that their second component ψ2bears a vortex of charge 1) and S= 1 vortices in the continuous setting. Numerical analysis has been performed in a similar fashion as in Ref. [30], using spectral methods for dealing with spatial derivatives. Figure 1shows several examples of the profiles for S= 0 and S= 1 stationary solutions. To assess stability, Fig. 2shows the dependence on the frequency of the real and imaginary parts of the eigenvalues for S= 0 and S= 1 solitary waves. We observe that, similar to the Soler model [30], the continuum S= 0 solitons are unstable below a critical frequency (ω= 0.388) for this model. However, it is interesting to observe that contrary to the Soler model, only q= 0 instabilities are present for the S= 0 case; these instabilities are of exponential nature and, consequently, can be predicted by the Vakhitov– Kolokolov criterion, as the curve representing charge versus frequency presents a maximum at the bifurcation point (see bottom panel of Figs. 2). In addition, solutions get more localized with decreasing frequency and tend to be localized at r= 0 as ω→0. It is important to also note that the continuous spectrum of this problem extends throughout the imaginary axis, except for the interval (−m+|ω|, m −|ω|)[33]. Moreover, ω= 0, belonging to the point spectrum due to symmetries (such as translation and phase invariance) is at the middle of the spectral gap. We can see from the stability analysis of Fig. 2that for the fundamental solution, there are wide intervals of stability/instability for S= 0, while the S= 1 waveform is unstable for all values of ω; cf. with the similar results of [30]. This prompts us to turn to the dynamical evolution of the instability. The study of the dynamics of the unstable solutions for S= 0 shows that they feature collapse, as is depicted in Fig. 3. As a numerical diagnostic for the accuracy of the simulation, we have monitored the relative norm error εdefined as ε(t) = |Q(t)−Q(0)| Q(0) ,(17) with Q(t)being the soliton’s charge. We have observed that the norm is preserved within a factor ∼10−7, despite an obvious departure of the norm from its initial value as collapse occurs. In [35] one can find a movie with the soliton evolution. In every movie, the top panels correspond to the density of each spinor component and bottom ones to their phase. B. Discrete model We now turn to the existence, stability and dynamics of solitons in the discrete 2D nonlinear Dirac equation (NLDE) in the form of Eq. (9). Given the generic instability of S= 1 states in the continuum, in our discrete studies, we will restrict Page 4 of 13AUTHOR SUBMITTED MANUSCRIPT - JPhysA-107720.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript 5 0 5 10 15 0 0.5 1 1.5 2 2.5 r u(r) ω=0.9 ω=0.7 ω=0.5 ω=0.3 0 5 10 15 0 0.2 0.4 0.6 0.8 1 r v(r) ω=0.9 ω=0.7 ω=0.5 ω=0.3 0 5 10 15 0 0.5 1 1.5 2 r u(r) ω=0.9 ω=0.7 ω=0.5 ω=0.3 0 5 10 15 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 r v(r) ω=0.9 ω=0.7 ω=0.5 ω=0.3 FIG. 1: Radial profiles of the spinor components for (upper panels) S= 0 solitary waves and (lower panels) S= 1 vortices for different values of ω. considerations to states without vorticity in the first spinor component. The properties of these states are to be compared with the S= 0 continuum counterparts. We start by the case when only nearest neighbour coupling (i.e. ν= 0) is considered. The solutions are obtained by making use of fixed point methods based on the anti-continuous (AC) limit [32], that allows to solve Eq. (13). One of the main difficulties in this case is to identify a suitable solution in the AC limit given that there are many solutions at C= 0 that can be extended to the continuous C→ ∞limit. In fact, all the solutions we have analyzed can be traced upon increasing the coupling strength Ctowards the continuum limit. Among all the solutions in the AC limit, remarkably we find that the one that leads to the S= 0 solitary waves of the continuum limit – discussed in the previous section – is the 9-site soliton, i.e. u0,0=u0,±1=u±1,0=u±1,±1=√1−ω. The vfield must be vanishing at C= 0, as can be seen from Eq. (13). The left panel of Fig. 4shows the profile of a typical such solution at finite coupling. Notice that in the left panels of Fig. 4, one can observe that, contrary to the soliton in the continuum, the imaginary part of un,m is not null; however, ||Im(un,m)||∞tends asymptotically to 0 when C→ ∞. An additional advantage of the AC limit is that the decoupled nature of the lattice enables us to analytically calculate the spectrum at C= 0. This consists of 9 pairs of eigenvalues at 0, N2−9pairs at λ=±i(1 + ω)and N2pairs at λ=±i(1 −ω). As in the 1D case, some of the eigenvalues become real when the coupling is switched on and the soliton is unstable for finite coupling. In particular, 7 eigenvalue pairs detach from zero whereas 2 additional pairs remain at zero for every coupling. Of the seven remaining pairs, two become imaginary and five pairs acquire a real part; for very low coupling C.0.01, four among these real pairs experience Hamiltonian Hopf bifurcations yielding complex quartets. This scenario persists for C < 0.80. Beyond this value, a rather complicated scenario arises, as can be seen in Fig. 5, where the dependence of the stability eigenvalues as a function of the coupling constant Cfor ω= 0.7is shown. One can observe that for large Cthere are only two sources of oscillatory instabilities: (1) for C&2.7, one of the quartets that exists for small Cpersists even for high C; (2) at C≈8.8, an additional Hamiltonian Hopf bifurcation takes place. Contrary to the Soler model [31], the 2iω eigenvalue is not present either in the discrete or in the continuum limit. It is intriguing to note that for the discretization considered at least one of the Hamiltonian Hopf bifurcations persists for large values of C. Such a feature has been previously encountered in the finite difference analysis of the Soler model in [31]. Another interesting feature that occurs below a certain critical frequency is the existence of two branches of solitons, one starting at C= 0 and another one finishing at C→ ∞. These branches are connected by an intermediate one through two saddle-center bifurcations i.e., the dependence of Qon Chas an S-shaped form. This type of bifurcation point arises for higher values of C when ωdecreases. Figure 6shows the relevant dependence of the charge with respect to Cfor ω= 0.7and ω= 0.6. We can observe in the latter case the existence of this intermediate branch, as well as the associated fold points. In a similar fashion to the 1D case [31], 1-site solitons also exist and can be continued to the continuum limit. In the AC limit, Page 5 of 13 AUTHOR SUBMITTED MANUSCRIPT - JPhysA-107720.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript 6 0.3 0.4 0.5 0.6 0.7 0.8 0.9 −1 −0.5 0 0.5 1 Im(λ) ω 0.3 0.4 0.5 0.6 0.7 0.8 0.9 −1 −0.5 0 0.5 1 Im(λ) ω 0.3 0.4 0.5 0.6 0.7 0.8 0.9 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 ω Re(λ) q=0 q=±1 q=±2 q=±3 q=±4 0.3 0.4 0.5 0.6 0.7 0.8 0.9 −1 −0.5 0 0.5 1 Re(λ) ω 0.3 0.4 0.5 0.6 0.7 0.8 0.9 5.5 6 6.5 7 7.5 8 8.5 9 Q ω 0.3 0.4 0.5 0.6 0.7 0.8 0.9 24 26 28 30 32 34 Q ω FIG. 2: Dependence of the (top) imaginary and (middle) real part of the eigenvalues with respect to ω. Left (respectively, right) panels correspond to S= 0 solitary waves (S= 1 vortices). The panels at the left illustrate that the S= 0 solution is only unstable for ω < 0.388. The panels at the right indicate that the S= 1 solution is unstable for all values of the frequency ω. The color that corresponds to each qis indicated in the middle left panel. Bottom panels show the dependence of the charge with respect to the frequency. According to the Vakhitov–Kolokolov criterion, the maximum of the curve indicates the occurrence of an exponential instability caused by radially symmetric perturbations (q= 0). they are given by u0,0=√1−ω, while v0,0= 0. The spectrum at C= 0 comprises N2−1pairs at λ=±i(1+ω),N2pairs at λ=±i(1 −ω)and a sole pair at λ= 0. As the single pair at 0 must remain there because of the U(1) symmetry, no eigenvalue can depart from 0 and the soliton is stable, at least for small values of the coupling parameter C. The middle panel of Fig. 4shows the profile of a 1-site soliton for ω= 0.7and C= 1. Notice that such solitons, for every coupling, possess the following properties: •Re(un,m) = 0 if nand mare odd, •Im(un,m) = 0 if nand mare even, •Re(vn,m) = 0 if nis odd and mis even, •Im(vn,m) = 0 if nis even and mis odd, Page 6 of 13AUTHOR SUBMITTED MANUSCRIPT - JPhysA-107720.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript 7 t=1.0 −2 0 2 −2 0 2 t=13.0 −2 0 2 −2 0 2 t=15.0 −2 0 2 −2 0 2 t=16.0 −2 0 2 −2 0 2 t=16.2 −2 0 2 −2 0 2 t=16.4 −2 0 2 −2 0 2 t=16.6 −2 0 2 −2 0 2 t=16.8 −2 0 2 −2 0 2 t=17.0 −2 0 2 −2 0 2 0 5 10 15 −3.5 −3 −2.5 −2 −1.5 −1 −0.5 0x 10−7 ε(t) t FIG. 3: Evolution of the unstable S= 0 soliton with ω= 0.3. Dynamics towards the collapse with the conservation of charge is observed. The left panels illustrate the evolution of the total density at different times in a two-dimensional contour plot. The right panel shows the evolution of the power based diagnostic of Eq. (17) which is well conserved during the simulation. Re(un,m) −20 −10 0 10 20 −20 −10 0 10 20 0.2 0.4 0.6 0.8 1 1.2 Im(un,m)×103 −20 −10 0 10 20 −20 −10 0 10 20 −4 −2 0 2 4 Re(vn,m) −20 −10 0 10 20 −20 −10 0 10 20 −0.4 −0.2 0 0.2 0.4 Im(vn,m) −20 −10 0 10 20 −20 −10 0 10 20 −0.4 −0.2 0 0.2 0.4 Re(un,m) −20 −10 0 10 20 −20 −10 0 10 20 Im(un,m)×103 −20 −10 0 10 20 −20 −10 0 10 20 Re(vn,m) −20 −10 0 10 20 −20 −10 0 10 20 Im(vn,m) −20 −10 0 10 20 −20 −10 0 10 20 0 0.5 1 −0.5 0 0.5 −0.4 −0.2 0 0.2 0.4 −0.4 −0.2 0 0.2 0.4 Re(un,m) −20 −10 0 10 20 −20 −10 0 10 20 Im(un,m)×103 −20 −10 0 10 20 −20 −10 0 10 20 Re(vn,m) −20 −10 0 10 20 −20 −10 0 10 20 Im(vn,m) −20 −10 0 10 20 −20 −10 0 10 20 0 0.5 1 −0.5 0 0.5 −0.4 −0.2 0 0.2 0.4 −0.4 −0.2 0 0.2 0.4 FIG. 4: Profiles of the 9-site (left set of 2×2panels), 1-site (middle set of panels) and 2-site (right set of panels) soliton with ν= 0,ω= 0.7 and C= 1. resembling the properties of the soliton of Fig. 2 of [18] in the large Climit. This also endows the solitons’ real and imaginary parts with a “staggered” structure with alternating rows missing, as seen in the middle panels of Fig. 4. Figure 7shows the dependence of the complex eigenvalues on C. One can see that the soliton is stable (see the real part in the bottom left panel of the figure) for Cbelow 1.04. At this critical point there is a bifurcation caused by a mode that destabilizes, after bifurcating from the linear modes band. Above this point, the soliton is exponentially unstable, becoming stable again at C= 1.631. However, at C= 1.73 it experiences a similar instability anew, while the structure finally stabilizes ∀C≥3.0in this case of ω= 0.7. Another interesting structure is the 2-site soliton, which, in the continuum limit is reminiscent of the soliton of Fig. 5 in [18]. In the AC limit, such a wave structure is given by u0,0=u0,1=√1−ω, while once again the vfield is vanishing. The spectrum at C= 0 is composed by N2−2pairs at λ=±i(1 + ω),N2pairs at λ=±i(1 −ω)and two pairs at λ= 0. As in the 9-site case, the two pairs remain invariant at 0, and thus no eigenvalue departs from 0, making the structure spectrally stable for small C. The right panel of Fig. 4shows the profile of a 2-site soliton for ω= 0.7and C= 1. Here too, similarly to the 1-site soliton, we observe a staggered pattern with the following features: •Re(un,m) = 0 if nis odd, •Im(un,m) = 0 if nis even, •Re(vn,m) = 0 if nis even, •Im(vn,m) = 0 if nis odd. Figure 7shows the dependence of the linearization eigenvalues with C. One can see that the behavior is essentially the same as in the case of the 1-site soliton, with the bifurcations taking place at the same points as in the 1-site case. In fact, the spectrum is almost identical in both 1-site and 2-site case, as the panels of Fig. 7show. Page 7 of 13 AUTHOR SUBMITTED MANUSCRIPT - JPhysA-107720.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript 8 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 Im(λ) C 0 0.5 1 1.5 2 2.5 3 0 0.05 0.1 0.15 0.2 0.25 Im(λ) C 0 2 4 6 8 10 0 0.05 0.1 0.15 0.2 0.25 Im(λ) C 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 Re(λ) C 0 0.5 1 1.5 2 2.5 3 0 0.1 0.2 0.3 0.4 0.5 Re(λ) C 0 2 4 6 8 10 0 0.1 0.2 0.3 0.4 0.5 Re(λ) C FIG. 5: Dependence of the stability eigenvalues for the 9-site soliton with respect to Cfor ν= 0 and ω= 0.7. The left panels show the spectrum for a lattice with N= 20 sites in each dimension. As the system is too small, continuation for couplings higher than C= 1 is not pursued. In addition, computation of the full spectrum for large lattices is not possible. For this reason, the middle and right panels show a partial spectrum (i.e., the imaginary and real part of the eigenvalues with the lowest imaginary part) up to C= 3 and C= 10, respectively for N= 200. We also indicate with red (black) dots, the eigenvalues corresponding to oscillatory (exponential) instabilities, whereas blue dots correspond to stable eigenmodes. In each case, the top panels show the imaginary parts and the bottom panels the real parts (corresponding to unstable growth) of the relevant eigenvalues. 0 0.5 1 1.5 2 2.5 0 50 100 150 200 Q C 0 0.5 1 1.5 2 2.5 0 50 100 150 200 Q C 1.82 1.825 1.83 1.835 1.84 88 90 92 94 FIG. 6: Dependence of the charge (power) Qon the coupling constant for the 9-site solitons with ν= 0 and ω= 0.7(left) and ω= 0.6 (right). The inset in the latter zooms in the region where three branches coexist for a narrow interval of values of C. Following the ideas introduced in [34], one way of lifting the degeneracy between 1-site and 2-site solitons in the vicinity of the anti-continuum limit could be to introduce next-nearest-neighbour coupling. In order to establish this, we have analyzed the stability of such solutions for ν= 0.1when Cis varied. The outcome for ω= 0.7is shown in Fig. 8, where it can be observed that the change on the stability eigenvalues of the 1-site soliton is only quantitative whereas this change is qualitative for 2-site solitons. In fact, in the latter case, an eigenvalue bifurcates from λ= 0 towards the imaginary axis and leads to a Hopf bifurcation observed at C= 0.51 (through its collision with another pair of imaginary eigenvalues). This renders the solution unstable past this point; at C= 1.00 the solution experiences an exponential bifurcation similar to that of the ν= 0 case (and of the single site solution); finally, there is another exponential bifurcation at C= 1.22 so that above this point, 2-site solitons do not exist. Notice that there are two eigenvalues pairs at λ= 0 for ν= 0; when ν= 0, only one of the eigenvalues pairs remains Page 8 of 13AUTHOR SUBMITTED MANUSCRIPT - JPhysA-107720.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript 9 0 0.2 0.4 0.6 0.8 1 1.2 0 0.5 1 1.5 2 2.5 3 Im(λ) C 0 1 2 3 4 5 0 0.05 0.1 0.15 0.2 0.25 C Im(λ) 1−site 2−site 0 0.2 0.4 0.6 0.8 1 1.2 0 0.1 0.2 0.3 0.4 0.5 Re(λ) C 0 1 2 3 4 5 0 0.1 0.2 0.3 0.4 0.5 C Re(λ) 1−site 2−site FIG. 7: Dependence of the stability eigenvalues for the 1-site soliton (full blue line) and the 2-site soliton (dashed red line) with respect to Cfor ν= 0 and ω= 0.7. The left panels show the spectrum for a lattice with N= 20 sites in each dimension. The right panels show the relevant eigenvalues [i.e., only part of the spectrum as concerns especially Im(λ)] for N= 200. Notice that the two spectra are almost identical and, consequently, full blue and dashed red lines match almost perfectly. In every case, the top panels show the imaginary parts and the bottom panels the real parts (corresponding to unstable growth) of the relevant eigenvalues. at 0 when C&0, whereas the other one detaches leading to the Hopf bifurcation; the imaginary part of the departing eigenvalue grows faster with Cwhen νincreases (while it does not move if ν= 0). This is basically the mechanism that explains how the degeneracy is lifted: as the 1-site soliton only possesses one eigenvalue pair at λ= 0 at the AC limit, increasing νcannot have an effect of degeneracy breaking, contrary to what occurs in the 2-site case. Let us also mention that although staggered patterns for 1-site and 2-site solitons are also observed for ν= 0, there are no sites where the real/imaginary part of un,m or vn,m becomes exactly 0. We have also considered the stability of the above mentioned configurations for fixed coupling near the continuum limit and variable frequency and ν= 0. Figure 9shows the profile of such solutions for ω= 0.5and C= 5. Figure 10 shows, for C= 2.5and C= 5, the dependence of the complex eigenvalues with ωfor the 9-site solitons. Notice that solitons only exist above a critical value of ω(C); this critical value decreases with C. For instance, solitons only exist for ω≥0.47 (ω≥0.22) if C= 2.5(C= 5). This is caused by a bifurcation similar to that shown in Fig. 6. This bifurcation is not found (at least for the considered coupling C= 5) in the 1-site and 2-site solitons (see Fig. 11). In the latter case, we observe that while generally dynamical instabilities may arise for small values of ωas well as in a narrow interval in the vicinity of ω= 1, a wide parametric interval of frequencies exists where the solitary waves are dynamically stable. Finally, we have considered the dynamics of prototypical unstable solutions with ν= 0. To this aim, we have included the spectral plane of the corresponding solution and a link to the movie with the evolution. We start with the typical evolution for the 9-site soliton at small coupling, where it possesses several coexisting instabilities (see Fig. 12a). The evolution for C= 0.5and ω= 0.7, shown in [36], leads to the destruction of the structure. In the case of 1-site solitons, we consider two cases corresponding to each exponential instability (Figs. 12b-c); the first one leads to the expansion i.e. dispersion of the solitary wave ([37], for ω= 0.7and C= 1.45), and the second one to slow soliton motion ([38], for ω= 0.7and C= 1.9). For the 2-site soliton, where the spectral dependence on Cis qualitatively similar to the 1-site soliton, the dynamics for the first instability is similar to the 1-site case (i.e. the soliton expands with time); however, for the second instability (see spectrum in Fig. 12d) we observe that the 2-site soliton breaks up ([39], for ω= 0.7and C= 1.9), repartitioning its mass into predominantly single site structures. Another interesting regime for the dynamical observation of the solutions’ instability regards the setting close to the continuum limit. To this end, we have fixed C= 5 and observed the 9-site soliton dynamics for three cases (see Figs. 12e-g), which are traced in [40] (ω= 0.7), [41] (ω= 0.6) and [42] (ω= 0.3). In the first case, there is only an oscillatory instability that leads Page 9 of 13 AUTHOR SUBMITTED MANUSCRIPT - JPhysA-107720.R1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 Accepted Manuscript