The Two-Dimensional Fractional Discrete Nonlinear Schrodinger Equation
TThe Two-Dimensional Fractional Discrete Nonlinear Schr¨odinger Equation
Mario I. Molina
Departamento de F´ısica, Facultad de Ciencias, Universidad de Chile, Casilla 653, Santiago, Chile (Dated: July 6, 2020)We study a fractional version of the two-dimensional discrete nonlinear Schr¨odinger (DNLS)equation, where the usual discrete Laplacian is replaced by its fractional form that depends on afractional exponent s that interpolates between the case of an identity operator ( s = 0) and that ofthe usual discrete 2D Laplacian ( s = 1). This replacement leads to a long-range coupling among sitesthat, at low values of s , decreases the bandwidth and leads to quasi-degenerate states. The meansquare displacement of an initially-localized excitation is shown to be ballistic at all times with a‘speed’ that increases monotonically with the fractional exponent s . We also compute the nonlinearmodes and their stability for both, bulk and surface modes. The modulational stability is seen toincrease with an increase in the fractional exponent. The trapping of an initially localized excitationshows a selftrapping transition as a function of nonlinearity strength, whose threshold increases withthe value of s . In the linear limit, there persists a linear trapping at small s values. This behavioris connected with the decrease of the bandwidth and its associated increase in quasi-degeneracy. Introduction . Let us consider the discrete nonlinearSchr¨odinger (DNLS) equation in d -dimensions[1–3]: i dC n dτ + (cid:88) m C m + γ | C n | C n = 0 (1)where the sum is over the nearest-neighbor sites. Pa-rameters V and γ represent the coupling between near-est neighbor sites and the nonlinear coefficient, respec-tively. The DNLS equation has proven useful in describ-ing a variety of phenomena in nonlinear Physics, such aspropagation of excitations in a deformable medium[4, 5],dynamics of Bose-Einstein condensates inside coupledmagneto- optical traps[6, 7], transversal propagation oflight in waveguide arrays[8–11], self-focusing and collapseof Langmuir waves in plasma physics[12, 13] and descrip-tion of rogue waves in the ocean[14], among others. Itsmain features include the existence of localized nonlinearsolutions with families of stable and unstable modes, theexistence of a selftrapping transition[15, 16] of an initiallylocalized excitation, and a degree of excitation mobilityin 1D[17]. All these characteristics have made the DNLSinto a paradigmatic equation that describes the propaga-tion of excitations in a nonlinear medium under a varietyof different physical scenarios.Recently, the topic of fractional derivatives has gainedincreased attention. It started with the observation thata usual integer-order derivative could be extended to afractional-order derivative, that is, ( d n /dx n ) → ( d s /dx s ),for real s , which is known as the fractional exponent. Thetopic has a long history, dating back to letters exchangedbetween Leibnitz and L’Hopital, and later contributionsby Euler, Laplace, Riemann, Liouville, and Caputo toname some. In the Riemann-Liouville formalism[18–20],the s-th derivative of a function f ( x ) can be formallyexpressed as (cid:18) d s dx s (cid:19) f ( x ) = 1Γ(1 − s ) ddx (cid:90) x f ( x (cid:48) )( x − x (cid:48) ) s dx (cid:48) (2)for 0 < s <
1. For the case of the laplacian operator∆ = ∂ /∂ r , its fractional form ( − ∆) s can be expressed as[21] ( − ∆) s f ( x ) = L d,s (cid:90) f ( x ) − f ( y ) | x − y | d +2 s (3)where, L d,s = 4 s Γ[( d/
2) + s ] π d/ | Γ( − s ) | (4)where Γ( x ) is the Gamma function, d is the dimension,and 0 < s < s decreases.This causes an increase in the system’s degeneracy andaffects its capacity to selftrap excitations. The model . For a square lattice ( d = 2), the kineticenergy term (cid:80) m C m in Eq.(1), can be written as 4 C n +∆ n C n where ∆ n is the discretized Laplacian∆ n C n = C p +1 ,q + C p − ,q − C p,q + C p,q +1 + C p,q − , (5)where n = ( p, q ). Equation (1) can then be written indimensionless form as i dC n dt + 4 C n + ∆ n C n + χ | C n | C n = 0 . (6) a r X i v : . [ n li n . PS ] J u l m y m y m y s=0.1 s=0.4 s=0.9 Figure 1. Effective coupling K ( m , m ) for several fractionalexponents. N = 11 × where t ≡ V τ and χ ≡ γ/V . We proceed now to replace∆ n by its fractional form (∆ n ) s in Eq.(6). The form of this fractional discrete Laplacian for d = 2 is givenby[35, 36](∆ n ) s C n = (cid:88) m (cid:54) = n ( C m − C n ) K s ( n − m ) (7)where, K s ( m ) = 1 | Γ( − s ) | (cid:90) ∞ e − t I m (2 t ) I m (2 t ) t − − s dt (8)with m = ( m , m ) and I m ( x ) is the modified specialBessel function. An alternative expression for (∆ n ) s is(∆ n ) s C j = L ,s (cid:88) m (cid:54) = j ( C m − C j ) G , , (cid:16) / , − ( j − m s,j − m s )1 / s,j − m , − ( j − m (cid:12)(cid:12)(cid:12) (cid:17) (9)where j = ( j , j ) and m = ( m , m ), and G ( ... ) is theMeijer G-function. The symmetric kernel K s ( m ) playsthe role of a long-ranged coupling. Near s = 1, K ( m ) → δ m , u where u = (1 ,
0) or u = (0 , I ν ( z ) ∼ πν (cid:16) ez ν (cid:17) ν → ∞ . (10)plus Γ( n ) ∼ n / ( n/e ) n and Γ( n + s ) ∼ n s Γ( n ) for n → ∞ , one obtains K s ( m ) ∼ (cid:18) m + m m m (cid:19) − ( m + m ) | Γ( − s ) | ( m + m ) m + m m m m m (11)where m = ( m , m ). Thus, for the ‘diagonal’ case m = m = m , one obtains K s ( m ) ∼ | Γ( − s ) | − m √ m ( m → ∞ ) . (12)This decay is faster than in the one-dimensional case[34].We consider now stationary modes defined by C n ( t ) = e iλt φ n , which obey( − λ + 4) φ n + (cid:88) m (cid:54) = n ( φ m − φ n ) K s ( m − n ) + χ | φ n | φ n = 0(13)It should be mentioned that in expressions (6) and (13)the term 4 is to be replaced by 3 (2) for sites at theedge (corner), when dealing with a finite square lattice.Figure 1 shows K ( m ), where we see how the range of thecoupling increases as s decrease. This has the effect ofincreasing the coupling between distant sites, leading todeep consequences, as we will show below. Plane waves . Let us set χ = 0 and look for plane wavesolutions, φ n = A exp( i k . n ), where we are assuming an band w i d t h k=(1,0), (0,1)k=(1,1) ba lli s t i cs peed Figure 2. Left: Bandwidth along three k -space directions, asa function of the fractional exponent s . Right: Character-istic ballistic speed (square) as a function of the fractionalexponent. Note that as s →
1, this speed approaches 2, asexpected. ( N = 11 × infinite square lattice. After a short algebra, one obtainsthe dispersion relation, λ ( k ) = 4 + (cid:88) m (cid:0) e i k . m − (cid:1) K s ( m ) (14)Unfortunately, in this case it is not possible to rewriteEq.(14) in closed form in terms of special functions, as inthe one-dimensional case. Figure 2 shows the bandwidthalong two different directions in k -space. For both casesthe bandwidth decreases monotonically as s decreases.This flattening of the band increases the degeneracy ofthe modes, and the system becomes closer to an idealmodel known as the ‘simplex’[37, 38] where every siteis coupled to every other site with equal strength. Thisleads to a strong localization of an initially localized ex-citation. In our case, this will become evident when welook at selftrapping. Root mean square (RMS) displacement . One of theways to quantify the propagation of an excitation across a - (cid:1) (cid:1) - (cid:1) (cid:1) k x k y - (cid:1) (cid:1) - (cid:1) (cid:1) k x k y - (cid:1) (cid:1) - (cid:1) (cid:1) k x k y - (cid:1) (cid:1) - (cid:1) (cid:1) k x k y s=0.1 s=0.4 s=0.1s=0.7 s=1 Figure 3. Contour plot of λ ( k x , k y ) for several fractional ex-ponents. N = 11 × lattice, is by the Root Mean Square (RMS) displacement,defined as (cid:104) n (cid:105) = (cid:88) n n | φ n | / (cid:88) n | φ n | . (15)A general result concerning RMS in lattices is that if λ ( − k ) = λ ( k ), λ ( k ) = λ ( k + 2 π q ) and ∇ λ ( k ) | F BZ =0, where FBZ is the first Brillouin zone, then it can beproven that the RMS for an initially localized excitation( φ n (0) = δ n , ) is always ballistic and given (in 2D) by[39] (cid:104) n (cid:105) = (cid:20) π (cid:90) FBZ ( ∇ λ ( k )) d k (cid:21) t (16)where FBZ is the first Brillouin zone. Using the form of λ ( k ) given by Eq.(14), we obtain (cid:104) n (cid:105) = (cid:88) m m K s ( m ) K s ( − m ) (17)The parity properties of the modified Bessel functions,imply that K s ( − m ) = K s ( m ). Therefore, (cid:104) n (cid:105) = (cid:34) (cid:88) m m ( K s ( m )) (cid:35) t (18)where m = ( m , m ). The quantity inside the squarebrackets can be interpreted as the square of a char-acteristic ballistic speed. Given the rapid decrease of K s ( m ) with distance, (cid:104) n (cid:105) becomes well-defined for all0 < s <
1. Figure 2b shows this speed square as a func-tion of the fractional exponent s . The speed rises fromzero at s = 0 up to 4 at s = 1, which is the usual ballisticvalue. The vanishing of the speed at s = 0 implies that,in this limit, the excitation is unable to move and remainlocalized at the initial site. In fact, at s = 0 any initial m x m y m x m y m x m y m x m y s=0.1s=0.2s=0.4s=1 Figure 4. Fundamental linear mode for several values ofthe fractional exponent. Left: Contour plots. Right: Threedimensional plots. N = 11 × condition, which is a combination of all plane waves, willbe unable to diffuse way. This result is in consonancewith the observation that, at s →
0, the band becomescompletely degenerate, λ ( k ) → ∇ λ ( k ) to vanish. Nonlinear modes . Let us consider now nonlinear modes( χ (cid:54) = 0), i.e., solutions to Eqs.(13). They constitute asystem of N × N nonlinear algebraic equations, wherethe form of the nonlinear term adopted here correspondsto Kerr nonlinearity found in coupled waveguide arrays,as well as in the semi classical description of an elec-tron propagating in a deformable lattice. Numerical so-lutions are obtained by the use of a multidimensionalNewton-Raphson scheme, using as a seed the solutionobtained from the decoupled limit, also known as theanticontinuous limit. We take a finite N × N latticewith open boundary conditions and examine two-modefamilies, “bulk” modes, located far from the boundariesand “surface” modes located near the beginning (or end)of the lattice. The stability of these nonlinear modesis carried out in the standard manner, which we sketchhere for completeness: We perturb our stationary solu-tion C n ( t ) = ( φ n + δ n ( t )) exp( iλt ), where | δ n ( t ) | (cid:28) | φ n | .We replace this in the evolution equation (6) [ with ∆ n re-placed by (∆ n ) s ]. After a linearization procedure, wherewe neglect any powers of δ n ( t ) beyond the linear one,we obtain a linear evolution equation for δ n ( t ). Next,we decompose δ n ( t ) into its real and imaginary parts: δ n ( t ) = x n ( t ) + i y n ( t ). To simplify the notation wemap the sites of the two-dimensional square lattice intothose of an open one-dimensional chain: ( n , n ) → n ≡ n +( n − N , with 1 ≤ n , n ≤ N . Then, the equationsfor x n ( t ) and y n ( t ) can be written in the form: d dt (cid:126)x − A B (cid:126)x = 0 , d dt (cid:126)y − B A (cid:126)y = 0 , (19)where (cid:126)x = ( x , x , · · · x N × N ) and (cid:126)y = ( y , y , · · · y N × N ).Matrices A and B are given by A nm = [ − λ − K s (0) + (cid:15) n + 2 | φ n | − φ n ] δ nm + K s ( n − m )(20) B nm = [ λ + K s (0) − (cid:15) n − | φ n | − φ n ] δ nm − K s ( n − m )(21)where (cid:15) n = α − (cid:80) j K s ( n − j ), with α = 2 for a corner site, α = 3 for an edge site, and α = 4 for a bulk site. Linearstability is determined by the eigenvalue spectra of thematrices AB (or BA ). When all eigenvalues are real andnegative, the system is stable, otherwise it is unstable.In the more general case that considers possible complexeigenvalues, one defines the instability gain G as: G = Max of (cid:26) (cid:16) Re( g ) + (cid:112) Re( g ) + Im( g ) (cid:17) (cid:27) / (22)for all g , where g is an eigenvalue of AB ( BA ). Thus,when G = 0, the mode under inspection is stable; other-wise it is unstable. Results from the above procedure areshown in figure 5 which shows power versus eigenvaluebifurcation diagrams for some bulk and surface modes,and for several values of the fractional exponents. Alsoshown are generic shapes of the two-dimensional modes(in this case we have taken a much smaller lattice to re-duce computation time) . Since we are dealing with afinite square lattice with open boundaries, there are twosurface modes: the ‘edge’ one, and the ‘corner’ one. It isobserved (not shown) that, after some few layers belowthe boundary, the surface modes become almost indistin-guishable from the bulk ones. Thus, there is a continuoustransition from surface to bulk modes. As for the bulkmodes, we have focused on the ‘odd’ mode, centered ona single site, the ‘even’ mode, centered on two nearest-neighbor sites, and the ‘ring’ mode, centered around aclosed loop of 4 sites (the two first names originate fromthe usage employed for one-dimensional lattices). Forthe bulk modes, we observe that the fundamental modeis always stable for all s values while the even and ring P O W E R (a) oddevenrings=0.2 P O W E R (b) oddevenrings=0.6 P O W E R (c) oddevenrings=1 P O W E R (d) edgecorner s=0.2 P O W E R (e) corner edges=0.6 P O W E R (f) corneredge s=1 Figure 5. Power content versus eigenvalue for some bulk andsurface modes, for different fractional exponents. Continuous(dashed) curves denote stable (unstable) modes. N = 5 × modes are unstable. All the bulk curves seem to touchthe edge of the linear band, at low powers. As s is de-creased, the modes become wider, but there are not otherdramatic changes on the shape of the modes. The surfacemodes decay quickly away from the boundary, but theirbifurcations curves look rather similar as s is varied.An interesting special case to examine is that ofthe stability of a nonlinear uniform solution. In one-dimension, the instability of this mode has been observedto give rise to discrete solitons and has, in fact, beenproposed as a practical way to produce them[40]. Letus consider a solution of the form C n = φ exp( iλt ).After replacing into the evolution equation one obtains λ = 4 + χφ . Thus, C n ( t ) = φ exp( i (4 + χφ )). Afterinserting this form into Eqs.(20) and (21), one obtains A nm = [ − (cid:88) j (cid:54) = n K s ( n − j )] δ nm + K s ( n − m ) (23) B nm = [ − (cid:88) j (cid:54) = n K s ( n − j )+3 χ φ ] δ nm + K s ( n − m ) . (24)As before, we look at the eigenvalues of AB (or BA )and record the instability gain G . Figure 6a shows G asa function of the nonlinearity strength, χ | φ | , for severalfractional exponents. We notice that, the stable region( G = 0) increases with an increase in the fractional ex-ponent s . A possible explanation could be that, as s isdecreased form a large value (i.e., near s = 1), the rangeof the coupling among sites increases, causing that anyperturbation on a given site is instantly felt on distantsites. This situation of mutual, long-range perturbationsinhibits the stability of the uniform front, as opposed thatcase when a perturbation of a given site only affects itsimmediate neighbors. Selftrapping . One the well-known effects of the Kerrnonlinearity is the onset of selftrapping where, for a non-linearity strength above a threshold value, an initially lo-calized excitation does not diffuse away completely whenplaced on a given site of the lattice. After some time,part of the excitation remains localized in the immediatevicinity of the initial site, while the rest diffuse away in aballistic manner. We want to examine the possible effectof a fractional exponent on this trapping phenomenon.To ascertain the presence of a selftrapping transition, oneexamines the long-time average probability at the initialsite (site ‘zero’) (cid:104) P (cid:105) = 1 T (cid:90) T | C ( t ) | dt (25)We have computed (cid:104) P (cid:105) for several s values, comparingthe different selftrapping curves. Results are shown inFig.6b. We see that, as s is decreased from unity, theposition of the selftrapping transition decreases as welland, in the linear limit we notice a degree of linear trap-ping that increases for lower values of s . These resultscan be explained as follows: as s is decreased, so does thewidth of the band (Fig.2). This flattening of the bandsoriginate a smaller group velocity of the modes. On theother hand, the local nonlinearity is roughly equivalentto a linear impurity of strength χ | φ | . Thus, we havean effective linear impurity embedded in a lattice whosemodes have low group velocity (because of relatively flatband). The combination of these two effects, facilitatesthe trapping of the excitation and thus, decreases thenonlinearity threshold needed. There is yet another ef-fect we can see in Fig.6: As s decreases, the amount oftrapping in the limit of zero nonlinearity increases. Thislinear trapping approaches unity for s → (cid:104)| C | (cid:105) = ( N × N − + 1( N × N ) (26)Therefore, at large N the trapped fraction approachesunity. Conclusions . We have examined the consequences ofemploying a fractional version of the usual discrete Lapla-cian, parametrized by a fractional exponent, on the exis-tence and stability of nonlinear modes, the free propaga-tion of localized linear and nonlinear excitations, and the - χ | ϕ G A I N po w e r Figure 6. (a) Modulational instability gain versus nonlin-earity strength. From the leftmost to the rightmost curve s = 0 . , . , . , . , . , . , . , . N × N = 25) (b) Time-averaged probability at the initial site as a function of non-linearity, for several different values of the fractional expo-nent, s . From the leftmost to the rightmost curve s =0 . , . , . , . , . T = 20 , N × N = 81). selftrapping of initially localized excitations, on a two-dimensional square lattice. We found that the main ef-fect of a fractional exponent is to introduce a long-rangecoupling interaction among the sites of the lattice. Themean square displacement is always ballistic with a speedthat decreases with a decrease in the fractional exponent.At small values of s , one observes a decrease in the band-width with a corresponding increase in degeneracy. Thestability of the low-lying excitations is not dissimilar tothe one found in the one-dimensional case, while the mod-ulational stability increases with an increase in s . Finally,the selftrapping of an initially localized excitation showsa selftrapping transition with a threshold that increaseswith s , and shows a degree of linear selftrapping at low s values. This was explained to be a consequence of theincrease in degeneracy as the band gets flatter and flatteras s → different mathematical Laplacians .Also, the fact that the effect of the fractional Laplacianis to introduce a long-range coupling means that one canreproduce experimentally the effect of a fractional Lapla-cian in an optical context, by setting up an appropriatedistribution of refractive indices and inter site distancesbetween waveguides in a waveguide array. Thus, it ispossible in principle, to explore the fractional dynamicsvia optical experiments. ACKNOWLEDGMENTS
The author is grateful to Luz Roncal for valuable dis-cussions. This work was supported by Fondecyt Grants1160177 and 1200120. [1] P. G. Kevrekidis, The Discrete Nonlinear SchrodingerEquation (Springer, Berlin Heidelberg 2009).[2] J. C. Eilbeck, P. S. Lomdahl, and A. C. Scott. The dis-crete selftrapping equation, Physica D 16 (1985) 318-338.[3] J. C. Eilbeck, M. Johansson, The discrete nonlinearSchrodinger equation-20 years on, Proceedings of theConference on Localization and Energy Transfer in Non-linear Systems, Madrid, Spain (2002). (World Scientific,2003).[4] A. S. Davydov, Solitons and energy transfer along proteinmolecules, J. Theor. Biology 66 (1977) 379-387.[5] P. L. Christiansen and A. C. Scott (Eds). Davydov’sSoliton Revisited: Self-trapping of Vibrational Energyin Protein (Plenum Press, New York, 1990).[6] O. Morsh, M. Oberthaler, Dynamics of Bose-Einsteincondensates in optical lattices, Rev. Mod. Phys. 78 (2006)179-215.[7] V. A. Brazhniy, V. V. Konotop, Theory of nonlinear mat-ter waves in optical lattices, Mod. Phys. Lett. B 18 (2004)627-651.[8] D. N. Christodoulides, R. I. Joseph, Discrete self-focusingin nonlinear arrays of coupled waveguides, Optics letters13 (1988) 794-796.[9] F. Lederer, G. I. Stegeman, D. N. Christodoulides, G.Assanto, M. Segev, Y. Silberberg, Discrete solitons inoptics, Physics Reports 463 (2008) 1-126.[10] Rodrigo A. Vicencio, Mario I. Molina, and Yuri S.Kivshar, Phys. Rev. E 70 (2004) 026602.[11] J. W. Fleischer, M. Segev, N.K. Efremidis, D.N.Christodoulides, Observation of two-dimensional discretesolitons in optically induced nonlinear photonic lattices,Nature 422 (2003) 147-150.[12] V.E. Zakharov, Collapse and self-focusing of Langmuirwaves, in: Handbook of Plasma Physics, Vol. 2, BasicPlasma Physics, eds. A.A. Galeev, R.N. Sudan, (ElsevierNorth-Holland 1984), pp. 81-121.[13] V. E. Zakharov, Collapse of Langmuir Waves, Sov. Phys.JETP 35 (1972) 908-914.[14] M. Onorato, A. R. Osborne, M. Serio, S. Bertone, FreakWaves in Random Oceanic Sea States, Phys. Rev. Lett.86 (2001) 5831.[15] M. I. Molina, G. P. Tsironis, Dynamics of self-trappingin the discrete nonlinear Schr¨odinger equation, PhysicaD 65 (1993) 267-273.[16] G. P. Tsironis, W. D. Deering, M. I. Molina, Applicationsof self-trapping in optically coupled devices, Physica D68 (1993) 135-137.[17] Rodrigo A. Vicencio, Mario I. Molina, and Yuri S.Kivshar, Phys. Rev. E 70 (2004) 026602.[18] R. Herrmann, Fractional Calculus - An Introduction forPhysicists (World Scientific Singapore 2014).[19] Bruce West, Mauro Bologna, Paolo Grigolini, Physics ofFractal Operators, (Springer 2003).[20] An Introduction to the Fractional Calculus and Frac-tional Differential Equations, por Kenneth S. Miller,Bertram Ross (Ed.) (John Wiley & Sons 1993). [21] N.S. Landkof, Foundations of Modern Potential Theory(Translated from the Russian by A.P. Doohovskoy), DieGrundlehren der mathematischen Wissenschaften, vol.180, (Springer-Verlag New York 1972.[22] R. Metzler, J. Klafter, The random walk’s guide toanomalous diffusion: a fractional dynamics approach,Phys. Rep. 339 (2000) 1-77.[23] I. M. Sokolov, J. Klafter and A. Blumen, Fractional ki-netics, Physics Today , 44-58 (2002).[24] G. M. Zaslavsky, Chaos, fractional kinetics, and anoma-lous transport, Phys. Rep. 371 (2002), 461-580.[25] N. C. Petroni and M. Pusterla, Levy processes andSchrodinger equation, Physica A , 824 (2009).[26] L. A. Caffarelli, A. Vasseur, Drift diffusion equations withfractional diffusion and the quasi-geostrophic equation,Ann. of Math. 171 (2010) 19031930.[27] P. Constantin and M. Ignatova, Critical SQG in boundeddomains, Ann. PDE 2 (2016) 1-42.[28] M. F. Shlesinger, G. M. Zaslavsky and J. Klafter, Strangekinetics, Nature 363 (1993), 31-37.[29] N. Laskin, Fractional quantum mechanics, Phys. Rev. E (2000) 3135. (2000).[30] N. Laskin, Fractional Schrodinger equation, Phys. Rev.E ,056108 (2002).[31] M. Allen, A fractional free boundary problem related toa plasma problem, Comm. Anal. Geom. , 1665 (2019).[32] H. Berestycki, J.-M. Roquejoffre and L. Rossi, The influ-ence of a line with fast diffusion on Fisher-KPP propa-gation, J. Math. Biol. 66 (2013) 743.[33] A. Bueno-Orovio, D. Kay, V. Grau, B. Rodriguez and K.Burrage, Fractional diffusion models of cardiac electri-cal propagation: role of structural heterogeneity in dis-persion of repolarization, Journal of The Royal SocietyInterface , 20140352 (2014).[34] M. I. Molina, “The Fractional Discrete NonlinearSchr¨odinger Equation”, M. I. Molina, Phys. Lett. A ,126180 (2020).[35] Oscar Ciaurri, Luz Roncal, Pablo Raul Stinga, Jose L.Torrea, Juan Luis Varona, Nonlocal discrete diffusionequations and the fractional discrete Laplacian, regu-larity and applications, Advances in Mathematics 330(2018) 688.[36] Luz Roncal, private communication.[37] Danilo Rivas, Mario I. Molina, Seltrapping in flat bandlattices with nonlinear disorder, Sci. Rep. , 5229(2020).[38] J. D. Andersen, V. M. Kenkre, Self-trapping and timeevolution in some spatially extended quantum nonlinearsystems: Exact solutions, Phys. Rev. B (1993), 11134.[39] A. J. Martinez, and M. I. Molina, Diffusion in infinite andsemi-infinite lattices with long-range coupling, J. Phys.A: Math. Theor. , 275204 (2012).[40] Y. S. Kivshar, and M. Peyrard, Modulational instabilitiesin discrete lattices. Phys. Rev. A46