Standing and traveling waves in a model of periodically modulated one-dimensional waveguide arrays

PHYSICAL REVIEW E 108, 024214 (2023) Standing and traveling waves in a model of periodically modulated one-dimensional waveguide arrays Ross Parker *and Alejandro Aceves Department of Mathematics, Southern Methodist University, Dallas, Texas 75275, USA Jesús Cuevas-Maraver Grupo de Física No Lineal, Departamento de Física Aplicada I, Universidad de Sevilla, Escuela Politécnica Superior, C/Virgen de Africa 7, 41011 Sevilla, Spain and Instituto de Matemáticas de la Universidad de Sevilla, Edificio Celestino Mutis, Avenida Reina Mercedes s/n, 41012 Sevilla, Spain P. G. Kevrekidis † Department of Mathematics and Statistics, University of Massachusetts Amherst, Amherst, Massachusetts 01003, USA (Received 26 January 2023; accepted 11 July 2023; published 16 August 2023) In the present work we study coherent structures in a one-dimensional discrete nonlinear Schrödinger lattice in which the coupling between waveguides is periodically modulated. Numerical experiments with single-site initial conditions show that, depending on the power, the system exhibits two fundamentally different behaviors. At low power, initial conditions with intensity concentrated in a single site give rise to transport, with the energy moving unidirectionally along the lattice, whereas high-power initial conditions yield stationary solutions. We explain these two behaviors, as well as the nature of the transition between the two regimes, by analyzing a simpler model where the couplings between waveguides are given by step functions. For the original model, we numerically construct both stationary and moving coherent structures, which are solutions reproducing themselves exactly after an integer multiple of the coupling period. For the stationary solutions, which are true periodic orbits, we use Floquet analysis to determine the parameter regime for which they are spectrally stable. Typically, the traveling solutions are characterized by having small-amplitude oscillatory tails, although we identify a set of parameters for which these tails disappear. These parameters turn out to be independent of the lattice size, and our simulations suggest that for these parameters, numerically exact traveling solutions are stable. DOI: 10.1103/PhysRevE.108.024214 I. INTRODUCTION The study of nonlinear lattice dynamics has been fundamental in advancing our understanding of, for example, light propagation in nonlinear optics [1] and the wave-function properties of atomic condensates [2]. In the former realm, the relevant models consider the propagation of light in coupled arrays of optical waveguides, while in the latter setting, they explore the evolution of the mean-field wave function in the context of deep optical lattices. In both scenarios, the universal model of interest (also considered as an envelope wave model in other discrete settings, including mechanical and electrical lattices [3,4]) has been the prototypical discrete nonlinear Schrödinger (DNLS) lattice [5]. Progressively, over the past few years, a topic that has been gaining significant traction has been the exploration of topological features in both linear and nonlinear systems exhibiting wave dynamics. Indeed, recent studies in a diverse host of fields including but not limited to photonics [6], phononics [7,8], metamaterials [9], and atomic physics *[email protected] †ke[email protected] [10] highlight unique dynamical properties resulting from the interplay of nonlinearity and topology. Relevant realizations of, e.g., Su-Schrieffer-Heeger lattice systems and associated anomalous edge states have also been recently proposed in the work of [11] with the potential for application in the context of topoelectrical metamaterials. Notably, their interplay has been leveraged to produce solitonic excitations and domain walls [12–19] and to generate robust states propagating on domain edges [20–24] that defy discreteness-induced barriers such as the famous Peierls-Nabarro barrier [25]. The resulting topologically protected states achieve unidirectional, uninhibited propagation around lattice defects in topological lattices [26]. These intense recent efforts have been summarized, e.g., in [27,28] and also in the very recent and detailed review of [29]. Among the many ongoing efforts in the field of topological photonics, we single out here a series of highly influential recent experiments of Rechtsman and co-workers [16,30–32]. Topological photonics has its roots in two seminal papers [33,34], where the authors delineated in detail the one-on-one correspondence between condensed-matter physics (CMP) and photonics. In particular, the propagation direction zplays the role of time in the original CMP setting and hence leads to the notion of pseudotime. For the same reason, the conjugate 2470-0045/2023/108(2)/024214(13) 024214-1 ©2023 American Physical Society ROSS PARKER et al. PHYSICAL REVIEW E 108, 024214 (2023) wave number is referred to as a pseudofrequency (t↔z;ω↔ k). The first of these works [16] showcased the experimental realization of Floquet solitons in a topological band gap, the numerical existence and stability of which we subsequently explored in [35]. More recently, such dispersive nonlinear systems with a coupling dependent on the evolution variable were proposed as a suitable realization of nonlinear Thouless pumps [31], and the topological properties of the bands such as the Chern number were argued to govern the resulting soliton motion. In [32] the analogy with the quantum Hall effect and the original proposal of the Thouless pump [36]was taken further by studying how nonlinearity acts to quantize transport via soliton formation and spontaneous symmetrybreaking bifurcations. In the present work, influenced by these studies, we consider the system analyzed in [32], but we depart from the adiabatic regime of focus in that work. By doing so, we are able to capture topologically induced stationary and dynamic states beyond the adiabatic approximation. We do so by enabling the computation of numerically exact stationary solutions but importantly also traveling solutions. Not only do we generate such waveforms by generic dynamical evolution experiments, but we also study a simple variant of the model which considers piecewise-constant coupling strengths (in a way reminiscent of the celebrated Kronig-Penney model [37]). There it becomes evident that at a qualitative level, the transition between standing and traveling waves mirrors the self-trapping transition of the DNLS dimer [38]. This may provide quite relevant insight towards understanding the symmetry-breaking transitions and dynamics within the intensely studied topic of nonlinear Thouless pumps. Our work is structured as follows. In Sec. II we present the theoretical setup and our quantitative diagnostics used in the model of interest. In Sec. II B we rescale the model so that we can use the propagation distance Las a parameter. We subsequently turn to numerical computations in Sec. III, starting with the evolution of single-site initial conditions and then gaining insights from the simplified piecewise-constant coupling model. In Sec. III A we perform evolution experiments which demonstrate that there are two fundamental behaviors to the system. For low initial power, the initial intensity moves either to the left or to the right in the lattice. The direction of motion depends on the lattice site chosen for the initial condition. For high initial power, the intensity remains confined to a single lattice site. In addition, there does not appear to be a sharp transition between these two behaviors when the starting intensity of the single site is continuously varied. In Sec. III B we consider a simplification of the model in which the coupling between waveguides is given by step functions. An analysis of this simplified model for an optical dimer explains these two observed behaviors, as well as the lack of a sharp transition between them. In Sec. III C we numerically construct both stationary and traveling coherent structures. As opposed to what occurs with single-site initial conditions, these coherent structures reproduce themselves exactly after an integer multiple of the coupling period. We use Floquet theory to determine the spectral stability of the stationary coherent structures, which are periodic orbits of the system. In Sec. IV we summarize our findings and present our conclusions, including a number of directions for future study. II. THEORETICAL ANALYSIS A. Mathematical model As discussed above, as well as motivated by experiments such as those of [31,32], we study light propagation in an array of coupled optical waveguides, where the coupling is periodically modulated along the axis of light propagation. Mathematically, this is described by the nonautonomous variant of the DNLS model of the form idun dz + m Hn,m(z)um+g|un|2un=0,(1) where un(z) is the complex amplitude of light propagating at the waveguide in the lattice site indexed n,zis the propagation distance (in the direction along the waveguides), and His the linear z-dependent (i.e., dependent on the evolution variable) tight-binding Hamiltonian or equivalently the lattice coupling profile. It is important to clarify here that the nonautonomous nature of the system under study is in connection with the notion of pseudotime (corresponding to the propagation distance), as indicated above. It is with that sense of nonautonomy in mind that we will proceed hereafter. The parameter gquantifies the strength of the cubic Kerr nonlinearity. For H,asin[32], we use an off-diagonal implementation of the Aubry-André-Harper model [39,40] with three sites per unit cell, resulting in the model idun dz +Jn(z)un+1+Jn−1(z)un−1+g|un|2un=0.(2) The z-dependent coupling functions Jn(z) are periodic in zwith spatial period L, which we will refer to as the coupling period. We note that [32] considers this model in the adiabatic regime, i.e., for very large L(see, for example, Fig. 2 in [32], where L=8000), in which case the system is approximately at a frozen equilibrium for every z. This is a central point to the analysis presented therein, which is explicitly geared towards (and limited to) such an adiabatic regime. By contrast, the parameter regime we consider herein is that of relatively small L, e.g., L=2π. In this case, stationary solutions are (genuine) periodic orbits of the system, which in turn enables us to use the tools of Floquet analysis to determine their spectral stability. While we consider larger Lbelow (see, in particular, Fig. 14), we note that our methods do not allow us to compute exact periodic orbits for which Lis very large, i.e., ones which would approach the adiabatic regime (see the end of Sec. III C 1 for further details). The choice of Jn(z), Jn(z)=J0+Ccos2π Lz+4π 3n+π 6,(3) groups the lattice sites into unit cells comprising three waveguides each (Fig. 1; see also Fig. 1 in [32]), since Jn(z)=Jm(z) for m≡n(mod 3). This choice of Jn(z) is slightly modified from Eq. (3) in the Supplemental Material of [32]; squaring the cosine function ensures that the coupling is always positive, which would be the case in a physical realization of the model. We note that if we do not square the cosine (as in the Supplemental Material of [32]), thus allowing for negative coupling values, we obtain the same qualitative behavior. When C=0, the nearest-neighbor coupling is constant and 024214-2 STANDING AND TRAVELING WAVES IN A MODEL OF … PHYSICAL REVIEW E 108, 024214 (2023) FIG. 1. Coupling functions J0(z), J1(z), and J2(z)ofz-dependent nearest-neighbor couplings over one spatial period L. Eq. (2) reduces to the ordinary DNLS equation with coupling via the discrete second difference operator Eq. (4), written in a corotating frame with frequency 2J0, which is given by idun dz +J0(un+1−2un+un−1)+g|un|2un+2J0un=0.(4) B. Model rescaling and density evolution In order to use the spatial period Las a parameter, we rescale the propagation distance using the change of variables z=LZ so that the coupling period is always 1, i1 L dun dZ +Jn(Z)un+1+Jn−1(Z)un−1+g|un|2un=0, Jn(Z)=J0+Ccos2πZ+4π 3n+π 6. (5) At any propagation distance Z, the power of the solution is its squared 2norm P(un)= n|un(Z)|2,(6) where the sum is taken over the entire lattice. The optical intensity at lattice site nis the square amplitude |un|2. The power is conserved, i.e., P(un) is independent of Z.UsingEq.(5) and its complex conjugate, we derive the flux equations for the density matrix elements ρmn =unum, dρmn dZ =iL[Jm(Z)ρm+1,n+Jm−1(Z)ρm−1,n −Jn(Z)ρm,n+1−Jn−1(Z)ρm,n−1] +iLg(ρmm −ρnn )ρmn.(7) The evolution of the optical intensity (or density) of the solution at lattice site n, which is given by ρnn =|un|2,is dρnn dZ =iL[Jn(Z)ρn+1,n+Jn−1(Z)ρn−1,n −Jn(Z)ρn,n+1−Jn−1(Z)ρn,n−1] =−2LIm[Jn(Z)ρn+1,n+Jn−1(Z)ρn−1,n],(8) 0 0.2 0.4 0 0.5 1 0 1 2 0 1 2 (b)(a) (d)(c) FIG. 2. Colormap showing the intensity of the solution of Eq. (5) evolving in Z, starting with a single excited site at n=0 with intensity (a) P=0.5, (b) P=1, (c) P=2.25, and (d) P=2.5. where we used the fact that ρn,m=ρm,n. We can split the last line of Eq. (8) into two components QL n(Z)=−2LJn−1(Z)Imρn−1,n, QR n(Z)=−2LJn(Z)Imρn+1,n, (9) where QL n(Z) and QR n(Z) are the flow of intensity into site n from the left and the right, respectively, and a positive sign indicates that intensity is flowing into site nfrom the labeled neighboring site. III. NUMERICAL COMPUTATIONS A. Single-site evolution As an initial experiment, we consider dynamical simulations starting with a single excited site at Z=0. Unless otherwise specified, the parameters in the section are g=1, L=2π,J0=0.05, and C=0.4 and the simulations are run on a finite lattice with m=30 lattice points, with periodic boundary conditions on the ends of the lattice (i.e., a ring, which allows for waves to loop around when they reach the boundaries). Evolution in Zis performed with ode45 in MATLAB using a tolerance of 10−9. For a single-site initial condition of sufficiently high power (above a threshold between P=2.25 and 2.5 for input intensity at n=0 and between P=2.15 and 2.25 for input intensity at n=−1), the energy remains concentrated at a single site and the resulting excitation appears to be stable [see Figs. 2(d) and 3(d)]. We will address the associated slight difference in the power threshold between initial excitations at n=0 and −1 below. As the initial power is lowered, this single-site solution becomes prone to mobility; in both cases, this leads the initially concentrated intensity to leak to the right within the lattice before dispersing throughout the lattice [Figs. 2(c) and 3(c)]. For lower-power initial conditions (between P=0.5 and 1), the initial intensity moves either to the left in the 024214-3 ROSS PARKER et al. PHYSICAL REVIEW E 108, 024214 (2023) 0 0.2 0.4 0 0.5 1 0 1 2 0 1 2 (a) (b) (c) (d) FIG. 3. Colormap showing the intensity of the solution of Eq. (5) evolving in Z, starting with a single excited site at n=−1 with intensity (a) P=0.5, (b) P=1, (c) P=2.15, and (d) P=2.25. lattice [for the initial excited site at n=0, see Figs. 2(a) and 2(b)] or to the right [for the initial excited site at n=−1, see Figs. 3(a) and 3(b)] before dispersing through the lattice. One explanation for this observation is as follows: For the first third of the period (Z∈[0,1 3]), the strongest coupling is between sites n=0 and n=−1viaJ−1=J2[see Eq. (3) and Fig. 1], so intensity can flow to the left from n=0to n=−1, which occurs when the initial power is sufficiently low. For Z∈[1 3,2 3], the strongest coupling is between n=−1 and n=−2, and for Z∈[1 3,2 3], the strongest coupling is between n=−2 and n=−3; thus we expect the intensity to travel three sites to the left over one period. Similarly, for the rightward-moving solution starting at n=−1, the rightward coupling is strongest for a rightward sequence of sites on the Zintervals [0,1 3], [2 3,1], and [4 3,5 3]; thus this solution moves to the right three sites in two periods. A similar rightward-moving behavior occurs when the initial excited site is n=1 (not shown). The first coupling for n=1isto the right on the interval [1 3,2 3]; the behavior is then similar to that of the rightward-moving solution for n=−1. We thus conclude that this fundamental difference between the leftward- and rightward-moving solutions and the associated speeds is a direct consequence of the form of the periodic coupling function, together with the lattice site at which the initial intensity is placed; this is suggestive also towards the difference in power thresholds noted above. In addition, we believe that this intuitive explanation is quite straightforward and thus appealing towards understanding solitary-wave motion in such nonautonomous discrete nonlinear systems. B. Simplified model A further heuristic and qualitative (yet in our view informative and intuitive) explanation of these different behaviors can be obtained by considering the simplification of the system of Eq. (2) obtained by approximating the coupling functions Jn(Z) with step functions, as is done, e.g., in the setting of FIG. 4. Simplified coupling functions Eq. (10)forZ∈[0,1]. the Kronig-Penney model [37]. Note that such a perspective has also been beneficial in a quantitative fashion in the case of nonlinearity (rather than dispersion) management in works such as [41,42]. Specifically, we define Jn(Z)=⎧ ⎨ ⎩ Cχ[0,1/3](Z),n≡2(mod3) Cχ[1/3,2/3](Z),n≡1(mod3) Cχ[2/3,1](Z),n≡0(mod3) (10) for Z∈[0,1] and extend periodically for all Z(Fig. 4). The function χ[a,b](Z) is the characteristic function of the interval [a,b], defined by χ[a,b](Z)=1forZ∈[a,b] 0 otherwise. Using this approximation, on the interval Z∈[0,1 3], the only active coupling is between the sites n≡0 (mod 3) and their left neighbors, which effectively creates a string of independent dimers. Our analysis follows Kenkre and Campbell’s study of self-trapping in a DNLS dimer [38]. Similar to what was done in that work, we will derive a second-order differential equation for the difference in intensity between the two sites of the dimer. This ordinary differential equation (ODE) will have solutions in terms of Jacobi elliptic functions and these will be used to explain the observed transition between mobile and trapped solutions as the coupling strength is increased. We note that this allows us to write the intensities |un|2(not the complex amplitudes un) in terms of Jacobi elliptic functions. Looking only at the sites n=0 and n=−1, i.e., one of these dimers, the system of equations Eq. (8) reduces to the four equations dρ0,0 dZ =iLC(ρ−1,0−ρ0,−1), dρ−1,−1 dZ =iLC(ρ0,−1−ρ−1,0), dρ−1,0 dZ =iL[C(ρ0,0−ρ−1,−1)+g(ρ−1,−1−ρ0,0)ρ−1,0], dρ0,−1 dZ =iL[C(ρ−1,−1−ρ0,0)+g(ρ0,0−ρ−1,−1)ρ0,−1]. (11) 024214-4 STANDING AND TRAVELING WAVES IN A MODEL OF … PHYSICAL REVIEW E 108, 024214 (2023) Let p=ρ0,0−ρ−1,−1be the difference between the intensities of site n=0 and site n=1. As in [38], we will derive a second-order ODE for p. Following the analysis in [43], we define q=i(ρ−1,0−ρ0,−1), r=ρ−1,0−ρ0,−1.(12) Using Eq. (11), we obtain the system of first-order ODEs dp dZ =2LCq,(13) dq dZ =−L(2Cp−gpr),(14) dr dZ =−Lgpq.(15) Since d dZ p2=2pdp dZ =4Lcpq, Eq. (15) becomes dr dZ =−g 4C d dZ p2, which has the solution r=r0−g 4Cp2−p2 0,(16) where r0and p0are the initial conditions for rand p,respectively, at Z=0. Differentiating Eq. (13) and substituting Eq. (14) and Eq. (16), we obtain the second-order differential equation for p, d2p dZ2=L2(Ap−Bp3), A=−4C2+2Cgr0+g2 2p2 0,B=g2 2.(17) This second-order autonomous differential equation with a linear term and a cubic nonlinearity has an exact solution in terms of Jacobi elliptic functions [see Sec. 22.13(iii) of [44] as well as [45]]. We are interested in the case of a single-site initial condition with intensity Pat site n=0, for which p0=Pand r0=0. Following [38], for fixed Cand g,Eq.(17) has an exact solution p(Z)=Pcn2CLZ;k=gP 4C,P<P∗ PdngPL 2Z;k=4C gP ,P>P∗,(18) where the critical intensity P∗is given by P∗=4C g.(19) The functions cn(z;k) and dn(z;k) are the Jacobi elliptic functions with elliptic modulus k(we note that these are often written in terms of the elliptic parameter m, where m=k2). Since the power of the solution is conserved, i.e., |u0(Z)|2+ |u−1(Z)|2=P, the intensity at sites n=0 and n=−1onthe interval Z∈[0,1 3] is given by |u0(Z)|2=1 2[P+p(Z)], |u−1(Z)|2=1 2[P−p(Z)].(20) (a) (b) FIG. 5. Plot of |u0|2in simplified model, given by Eq. (20), vs Z for initial intensity (a) P<P∗and (b) P>P∗. Although this solution holds only for Z∈[0,1 3](Z=1 3is marked with a solid vertical line), it is continued to Z=1 for illustrative purposes. The other parameters are C=0.5, g=1, and P∗=2. There are two fundamental behaviors of the single-site initial condition in the simplified model, depending on whether P< P∗or P>P∗. In the dimer, a sharp transition between the two behaviors occurs at P=P∗(see [38]). We note that this transition is somewhat blurred in this model, even with the simplified coupling function, since the initial dimer coupling is broken [Jn(Z)=0] at Z=1 3. Below we will discuss P< P∗and P>P∗separately. 1. Case 1: P <P∗ For P<P∗, the solution p(Z) involves the Jacobi cn function, which oscillates about 0 with period 2K(k)/CL, where K(k)=π/2 0 dθ 1−k2sin2θ (21) is the complete elliptic integral of the first kind. (We note that the period of oscillation becomes infinite as Papproaches P∗from below.) The intensity |u0|2, given by Eq. (20), exhibits large-amplitude oscillations with this period from 0 to P, centered at P/2 [Fig. 5(a)]. Intensity initially flows to the left; if the coupling is not cut off at Z=1 3(and no other couplings are activated), there is a critical value Z∗ 1of Zat which the intensity from site n=0 has been completely transferred to site n=−1; after this point, intensity starts flowing back in the other direction. This critical value Z∗ 1is larger for larger starting intensity P[see Fig. 5(a)]. For most configurations, including all of the examples in Fig. 5,the critical value Z∗ 1>1 3; thus there is still intensity remaining in site 0 when the coupling switches off at Z=1 3.Fig.6(a) plots the fraction of intensity that has been transferred from site n=0 to site n=−1atZ=1 3as the starting intensity varies. If this fraction is close to 1, numerical evolution experiments find leftward-moving solutions starting from a single-site initial condition [see Figs. 7(a)–7(c)]. As the initial intensity P approaches P∗from below, the fraction of intensity transferred at Z=1 3decreases to approximately 0.5 and we do not expect to see pure leftward-moving solutions. Numerical evolution experiments in this case show that the initial intensity splits into a leftward- and a rightward-moving solution [Fig. 7(d)]. We note that it is possible to choose parameters so that, for the simplified model, the period of p(Z)is2 3, i.e., Z∗ 1=1 3, so that the initial intensity has been completely transferred from n=0ton=−1 when the coupling switches off. The 024214-5 ROSS PARKER et al. PHYSICAL REVIEW E 108, 024214 (2023) (a) (b) FIG. 6. (a) Fraction of intensity transferred from site n=0 to site n=−1atZ=1 3for P<P∗for varying C. (b) Fraction of intensity remaining at site n=0atZ=1 3for P>P∗for varying C. For each curve, P∗is indicated with a star. The other parameter is g=1. next coupling is then between n=−1 and −2forZ∈[1 3,2 3]. This pattern continues, and so for this choice of parameters, the simplified model supports an exact leftward-moving solution which persists for a large interval in Z[Figs. 8(a) and 8(b)]. A comparison between the evolution of the original and simplified systems from single-site initial conditions for P<P∗, illustrating the similarity of the model dynamics, is shown in Figs. 9(a)–9(c). We will see below in Sec. III C 2 that left-moving coherent structures exist in the full model, but these do not have pure single-site initial conditions. We also note that by symmetry of Eq. (11), a similar analysis holds for rightward-moving solutions starting with single-site initial conditions at n=−1 on the interval Z∈[0,1 3]. 2. Case 2: P >P∗ For P>P∗, the solution p(Z) involves the Jacobi dn function, which oscillates about 1 with period 4K(k)/gPL, where K(k) is defined by Eq. (21). (We again note that the pe- (a) (b) (c) (d) FIG. 7. Colormap showing the intensity of the solution of Eq. (5) with the simplified coupling function Eq. (10)forP<P∗evolving in Z, starting with a single excited site at n=0 with intensity (a) P= 0.5, (b) P=1, (c) P=2, and (d) P=3. The fraction of intensity transferred from site n=0 to site n=−1atZ=1 3is (a) 0.9910, (b) 0.9959, (c) 0.9909, and (d) 0.6636. The other parameters are P∗= 3.2, C=0.8, and g=1. (a) (c) (b) (d) FIG. 8. (a) and (c) Colormap of evolution in Zof the single-site initial condition in the simplified model and (b) and (d) intensity |u0|2of site n=0 (right) on Z∈[0,1 3] (right). The parameters are chosen so that (a) and (b) Z∗ 1=1 3and (c) and (d) Z∗ 2=1 3. The starting intensity (a) and (b) P=0.0625 and (c) and (d) P=1.9453. The other parameters are C=0.75 and g=1. riod of oscillation becomes infinite as Papproaches P∗from above.) The intensity |u0|2exhibits small-amplitude oscillations (which become progressively smaller as Pis increased) with this period about the initial intensity P[Fig. 5(b)]. As in case 1, the intensity initially flows to the left. If the coupling is not cut off at Z=1 3(and no other couplings are activated), there is a critical value Z∗ 2of Zat which point the system has returned to its initial state, i.e., the intensities at sites n=0 and n=−1 are once again Pand 0, respectively. For most configurations, including all of the examples in Fig. 5,the critical value Z∗ 2= 1 3; thus when the coupling switches off, there is some net transfer of intensity to the neighboring site (a) (b) (c) (d) (e) (f) FIG. 9. Intensity |u0|2of site n=0onZ∈[0,1 3]forP<P∗ [(a) P=0.5, (b) P=1, and (c) P=1.5] and P>P∗[(d) P=2.5, (e) P=3.5, and (f) P=4.5] for the original system (blue solid line) and the simplified system (orange dotted line). The other parameters are C=0.5, g=1, and P∗=2. 024214-6 STANDING AND TRAVELING WAVES IN A MODEL OF … PHYSICAL REVIEW E 108, 024214 (2023) (a) (b) (c) (d) FIG. 10. Colormap showing the intensity of the solution of Eq. (5) with the simplified coupling function Eq. (10)forP>P∗ evolving in Z, starting with a single excited site at n=0 with intensity (a) P=4, (b) P=5, (c) P=6, and (d) P=7. The fraction of intensity remaining at site n=0atZ=1 3is (a) 0.9938, (b) 0.8853, (c) 0.9815, and (d) 0.9797. The other parameters are P∗=3.2, C=0.8, and g=1. n=−1. Fig. 6(b) plots the fraction of intensity remaining at site n=0atZ=1 3for varying C. If this fraction is close to 1, numerical evolution simulations show stationary solutions starting from a single-site initial condition; the closer this fraction is to 1, the longer these stationary solutions will persist before breaking up (see the case examples in Fig. 10). As in case 1, we note that it is possible to choose parameters so that, for the simplified model, the system returns exactly to its starting condition at Z=1 3, i.e., Z∗ 2=1 3 [Fig. 5(b)]. In this case, for a single-site initial condition with a specific starting intensity, the simplified model supports a localized in space, time-periodic solution which persists for a large interval in Z. A comparison between the evolution of the original and simplified systems for P>P∗is shown in Figs. 9(d)–9(f). We discuss genuinely time-periodic (numerically exact) coherent structures in the full model in Sec. III C 1. C. Coherent structures in the full model These evolution experiments are strongly suggestive of the fact that the system Eq. (5) supports two classes of coherent structures: localized in space, time-periodic solutions, which are centered at a particular lattice site, and moving solutions, which reproduce themselves exactly a specific number of sites to the left or to the right. Recall that for the rescaled system Eq. (5), the coupling period is 1. The stationary coherent structures will be periodic orbits whose period is a multiple of the coupling period, i.e., a positive integer N. We note that while it may be possible to find such solutions which have a noninteger period, Floquet analysis requires that the period of the solutions be commensurate with that of the coupling. Similarly, we will look for moving solutions that reproduce themselves, shifted left or right, after an integer period. For appropriate choices of system parameters, we can compute both types of solutions numerically. For both localized (i.e., nonmoving) and moving solutions, we use a shooting method with periodic boundary conditions imposed on Z, starting with a single-site initial guess. In addition, for the former case, we validate this method by using numerical parameter continuation with AUTO [46] to solve a periodic boundary value problem. Unless otherwise specified, the parameters in the section are the same as in the preceding one. 1. Stationary (nonmoving) solutions First, we look at the stationary solutions. At the anticontinuum limit (J0=0 and C=0), the lattice sites are decoupled and an initial intensity Pat lattice site nwill yield a standing-wave solution of frequency P, i.e., of the form un(Z)=√Pe2πiPZ . Since such a solution has period 1/Pand stationary solutions must have an integer period, these solutions will exist in a discrete family for every integer period N, i.e., approximately P=k/Nfor sufficiently large positive integer k. For period N=1 and the parameters in the preceding section, for example, we expect to have time-periodic nonmoving solutions for approximate integer intensities P⩾2 (see Figs. 11 and 12 for the first two of these solutions). By looking at the intensity [Figs. 11(c) and 12(c)] and the real part [Figs. 11(d) and 12(d)] of the central site n=0, we see that they are approximately standing waves with frequencies 2 and 3, respectively. We note that the stationary solutions do not decay to 0 with increasing |n|, but rather the tails exhibit small-amplitude oscillatory patterns [see Figs. 11(b) and 12(b)]; the specific pattern of oscillations depends on the lattice size (not shown). Looking at the sites adjacent to the central one, the left neighbor u−1peaks on the interval [0,1 3] when the coupling J2(z) is most active, and the leftward flow QL 0<0, indicating flow of intensity out of site 0 to the left. The right neighbor peaks on the interval [2 3,1] when the coupling J0(z) is most active, and the rightward flow QR 0<0, indicating flow of intensity out of site 0 to the right. Both QL 0 and QR 0are close to 0 on the interval [1 3,2 3], when neither has strong nearest-neighbor coupling. Since these nonmobile solutions are true periodic orbits with integer period, so their period is equal to or commensurate with that of the coupling, their spectral stability can be determined by Floquet theory. Numerical computation of the Floquet multipliers of the stationary solutions is shown in the insets in Figs. 11(a) and 12(a). The lower-power solution has two pairs of Floquet multipliers off the unit circle, which is characteristic of an oscillatory instability. Long-term evolution in Z[Fig. 11(f)] shows that this solution remains coherent until approximately Z=130. By contrast, the Floquet spectrum of the higher-power solution lies on the unit circle, indicating spectral stability. Long-term evolution in Z[Fig. 12(f)] shows that this solution is still coherent at Z=500. Using numerical parameter continuation, we start with the DNLS soliton at C=0 and slowly vary C, tracing the curves in Fig. 13 (J0=0.05 throughout). Since we are looking for solutions with period 1, the starting single-site intensity must take integer values so that the frequency of the DNLS standing wave at C=0 is commensurate with this period. In all cases, 024214-7 ROSS PARKER et al. PHYSICAL REVIEW E 108, 024214 (2023) (a) (b) (c) (d) (e) (f) FIG. 11. (a) Initial intensity |un(0)|2, with the inset showing Floquet multipliers, and (b) logarithm of the initial intensity for the stationary solution with an approximate power of 2. (c) Intensity and (d) real part of three central sites over one period. (e) Leftward flow QL 0and rightward flow QR 0of the intensity at the central site n=0 and (f) long-term evolution in Z. The other parameters are C=0.4 and J0=0.05. a turning point is reached, as Cis increased, at which point the parameter continuation in the coupling parameter Creverses direction. This turning point occurs at a larger value of Cfor solutions which start at a higher initial intensities at C=0. All stationary solutions initially have their Floquet spectrum confined to the unit circle and thus are spectrally stable. Spectral stability is lost at some point before the turning point observed in the graph, when Floquet multipliers collide and leave the unit circle, creating an oscillatory instability. Solutions on the upper branches of the bifurcation diagram are periodic solutions to the DNLS equation which are not pure standing waves. To leading order, these upper solutions are the sum of two Fourier modes, as opposed to standing waves, which comprise a single Fourier mode. Substituting the finite Fourier ansatz un(Z)= N  k=−N an,ke2πikz (a) (b) (c) (d) (e) (f) FIG. 12. Same as Fig. 11 but for the nonmoving solution with an approximate power of 3. into Eq. (5) and projecting onto each of the Fourier basis functions, we can obtain expressions for the coefficients an,k for each wave number k. A fast Fourier transform of the numerical solution on the upper branches suggests that the solutions at each site are composed predominantly of the modes with wave numbers 0 and 1. Thus, to leading order, these solutions are of the form un(Z)=an,0+an,1e2πiωZ, where the coefficients an,0and an,1satisfy J0(an+1,0+an−1,0)+a3 n,0+2a2 n,1an,0=0, J0(an+1,1+an−1,1)+a3 n,1−wan,1+2an,1a2 n,0=0. We note that if an,0=0 for all n, the second equation reduces to the DNLS equation, in which case the solution is a standing wave. We can also continue solutions in the coupling period L (Fig. 14). The intensity of the central peak, hence the overall power of the solution, decreases with increasing L; thus solutions with greater starting power at L=2πpersist for higher L. For example, solutions with (approximate) starting powers of 2, 3, and 4 at L=2πpersist up to L=11.7, 15.7, and 16.6, respectively (see Fig. 14). Again, there is a turning point where the continuation reverses directions, which occurs at larger Lfor higher-power branches. The central site for the 024214-8 STANDING AND TRAVELING WAVES IN A MODEL OF … PHYSICAL REVIEW E 108, 024214 (2023) (a) (b) FIG. 13. (a) Branches of stationary solutions with spatial period (in Z) of 1, obtained from numerical parameter continuation starting with the DNLS soliton at C=0. (b) Intensity (top row) and logarithm of the intensity (bottom row) of the initial condition, corresponding to C=0.05 at the labeled points on the bifurcation diagram in (a). Solid lines correspond to solutions with the Floquet spectrum contained in the unit circle and dotted lines correspond to solutions with some of the Floquet spectrum outside the unit circle. Solutions on other branches at these values of Care qualitatively similar. upper and lower branches of each loop has approximately the same intensity; the higher power of the upper branches is due to higher intensity in the tails of the solutions. For contrast, the spatial period of the solutions in Fig. 2 in [32] is L=8000, which simulates the adiabatic regime; since the power of the solution decreases as Lis increased by parameter continuation, we would have to start with a solution with extremely high power at L=2πto be able to reach such a large Lusing this method. Obtaining such solutions with a very large spatial period Ldirectly by a shooting method is similarly computationally impractical. Finally, we note that, while we have only considered nonmoving solutions with period of 1, stationary solutions do exist for other positive integer periods. For example, if we start with a single-site initial condition with intensity k/2 for positive odd integer k, we expect that we will obtain a stationary solution with period 2. (We have verified that this is the case for single-site initial conditions with intensities 3 2and FIG. 14. Parameter continuation in coupling period L. The plot shows the power of the solution vs L, starting with solutions of approximate power 2, 3, and 4 at L=2π. Solid lines correspond to solutions with the Floquet spectrum contained in the unit circle and dotted lines correspond to solutions with some of the Floquet spectrum outside the unit circle. The diagram is shown only for L⩾2π. The other parameters are J0=0.05 and C=0.25. 5 2.) While in principle this can be done for any integer period, it becomes computationally intractable for larger periods. 2. Moving solutions Next we look for moving solutions. For a given lattice size m, we find that leftward-moving solutions exist (Fig. 15) for all values of Cwithin an interval [CL(m),CR(m)] [see Figs. 16(b) and 16(c)]. These are true coherent structures, in that the entire solution reproduces itself exactly after one period, shifted three sites to the left. (In the numerical simulation, where we are using periodic boundary conditions on (a) (b) (c) (d) FIG. 15. (a) Initial intensity |un(0)|2and (b) logarithm of the initial intensity for the left-moving solution. (c) Intensity of the solution evolved in Zover a period for a few select sites and (d) space-time contour plot evolution of the intensity for the traveling wave. 024214-9