{'fraction_non_alphanumeric': 0.045508335582452715, 'fraction_numerical': 0.013079494259555298, 'mean_word_length': 4.331857427496126, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Tight-binding methods for general longitudinally driven photonic lattices -- edge states and solitons

A systematic approach for deriving tight-binding approximations in general longitudinally driven lattices is presented. As prototypes, honeycomb and staggered square lattices are considered. Time-reversal symmetry is broken by varying/rotating the waveguides, longitudinally, along the direction of propagation. Different sublattice rotation and structure are allowed. Linear Floquet bands are constructed for intricate sublattice rotation patterns such as counter rotation, phase offset rotation, as well as different lattice sizes and frequencies. An asymptotic analysis of the edge modes, valid in a rapid-spiraling regime, reveals linear and nonlinear envelopes which are governed by linear and nonlinear Schrodinger equations, respectively. Nonlinear states, referred to as topologically protected edge solitons are unidirectional edge modes. Direct numerical simulations for both the linear and nonlinear edge states agree with the asymptotic theory. Topologically protected modes are found; they possess unidirectionality and do not scatter at lattice defect boundaries.

A systematic approach for deriving tight-binding approximations in general longitudinally driven lattices is presented. As prototypes, honeycomb and staggered square lattices are considered. Timereversal symmetry is broken by varying/rotating the waveguides, longitudinally, along the direction of propagation. Different sublattice rotation and structure are allowed. Linear Floquet bands are constructed for intricate sublattice rotation patterns such as counter rotation, phase offset rotation, as well as different lattice sizes and frequencies. An asymptotic analysis of the edge modes, valid in a rapid-spiraling regime, reveals linear and nonlinear envelopes which are governed by linear and nonlinear Schrödinger equations, respectively. Nonlinear states, referred to as topologically protected edge solitons are unidirectional edge modes. Direct numerical simulations for both the linear and nonlinear edge states agree with the asymptotic theory. Topologically protected modes are found; they possess unidirectionality and do not scatter at lattice defect boundaries.

I. INTRODUCTION
In recent years the study of topological edge/interface and surface modes has been a very active area of research. Fundamental research in systems where timereversal symmetry is broken have been found to support unidirectional edge modes that remain intact over long distance and against defects; these systems include magneto-optics [1,2], photonics [3][4][5] and acoustics [6]. While these systems have very different underlying physics they all share some common properties; they have a boundary that separates two distinct regions. A boundary discontinuity alone can be enough to support edge states, but to generate topologically protected modes in these systems typically time-reversal symmetry must be broken. The manner in which this is accomplished is different for each system; here however we focus on optical beams propagating through longitudinally driven waveguide arrays.
In photonic lattices, regions of high refractive index can be etched into bulk silica using femtosecond laser etching techniques [15,16] and act as waveguides for the beam. One way to break time-reversal symmetry is to rotate these waveguides longitudinally, creating an array configuration that changes along along the direction of beam propagation. Such lattices, which are orthogonal to the direction of the beam propagation, have been constructed at optical frequencies and topologically protected edge modes have been experimentally observed [5].
One of the special features topological edge modes exhibit is their one-way scatter-free motion even when lat-tice defects are encountered. Such traveling modes are said to be topologically protected meaning that unidirectionality is preserved even under significant deformations of the lattice. Topological invariants, such as Chern number [17] or Zak phase [18], can be associated with these topological modes. The robust nature of topologically protected edge modes suggests that they will be useful in many applications where small imperfections are always present.
Photonic lattices that can support linear edge modes include honeycomb [5] and staggered square lattices [19,20] each of which has two lattice sites per unit cell. We refer to any lattice with two or more sites in a unit cell as non-simple. In this work we provide a direct route for deriving tight-binding equations that describe beam propagation in general longitudinally varying lattices with either simple or non-simple configurations.
As typical examples we analyze honeycomb and staggered square lattices, though more complex lattices can be considered within the framework we present. In doing so, we are able to study a wide range of lattice dynamics including periodic but non-synchronized (out of phase) waveguide motion with phase offset across the two sublattices or counter rotating lattices etc. The Floquet bands are associated with intriguing edge mode dynamics which include: flat (stationary), non-unidirectional, oscillatory, and simultaneous topological and non-topological eigenstates. These modes complement the well-known traveling edge states present when the lattice waveguides rotate in-phase with each other.
For a particular rotation pattern we find common structure in the dispersion bands of both the honeycomb and staggered square lattices. This suggests that for a particular waveguide rotation there is an associated edge wave dynamic that is independent of the underlying lattice configuration. Additionally, our tight-binding model incorporates the geometry of the individual waveguides. When the waveguides are stretched in a preferred direction we find non-topologically protected modes whose edge modes reflect at lattice defects, rather than simply moving around the defect like a protected mode.
We present an asymptotic analysis which reveals that along the lattice boundary the edge modes can behave as one-dimensional solutions to the nonlinear Schrödinger (NLS) equation. Hence we refer to these traveling modes as edge solitons, where applicable. The NLS equation is a universal model for the envelope description of dispersive waves in weakly-nonlinear media, like the systems we consider in this paper. The nonlinear edge mode envelope can be seen as a balance between the lattice induced linear dispersion effects and sufficiently strong beam focusing or defocusing. Indeed, many properties associated with the NLS equation can be found in photonic topological insulators, such as modulational instability [21], band gap solitons [22] and soliton propagation in helically driven photonic graphene [23]. Moreover, we construct nonlinear modes that, to leading order, satisfy the focusing NLS equation and whose corresponding linear dispersion is topologically protected. These modes travel stably around defects. Hence we term these solutions topologically protected edge solitons.
The outline of the paper is as follows. In Sec. II we derive a tight-binding model that describes paraxial beams in general longitudinally driven honeycomb and staggered square lattices. The linear dispersion bands for one-dimensional edge states are computed via Floquet theory in Sec. III. Here we study several complex sublattice rotation patterns and their effect on the edge mode dynamics. In Sec. IV we consider nonlinear edge beams. Both linear and nonlinear edge modes are found to be modulated by a slowly-varying envelope function that satisfy the linear and nonlinear Schrödinger equations, respectively. In the nonlinear regime an asymptotic analysis shows that edge solitons exist for solutions with narrow spectral support around some central frequency. The edge wave behavior across lattice defects is explored in Sec. V. There we observe topologically protected modes to be scatter-free at lattice barriers, in contrast to non-topologically protected modes. We conclude in Sec. VI.

II. DERIVATION OF TIGHT-BINDING APPROXIMATION
Nonlinear quasi-monochromatic light beams propagating through photonic lattices are described by the paraxial wave or nonlinear Schrödinger (NLS) equation where ∇ 2 ≡ ∂ 2 x + ∂ 2 y and ψ(r, z) is a complex envelope function defined on the transverse plane r = (x, y) and propagation direction z. The function ∆n(r, z) is the lattice potential and it models the variations in the refractive index. The physical parameters are the wavenumber k 0 = 2πn 0 /λ, bulk index of refraction n 0 , beam wave-length λ, deviation from bulk index |∆n|, and Kerr coefficient n 2 .
We may think of a lattice with two generating sites per unit cell as the combination of two interpenetrating sublattices with potentials V 1 (r) and V 2 (r). These two sublattices have regions of high refractive index concentrated at the white and black lattice sites, respectively, as shown in Figs. 1 and 2. To model such a lattice configuration we employ the potential function that consists of well-localized dips (minima) located at the lattice sites. We approximate these sublattices by the sum of Gaussians where the positions of the white and black lattice sites are given, respectively, by the sets Asymmetry in the waveguide geometry can be explored by adjusting the width parameters σ x and σ y . When σ x = σ y we say the lattice isotropic, whereas when σ x = σ y the lattice is anisotropic. Detuning (difference in the refractive index) between the sublattices can approximated by including a coefficient 0 < q < 1 in front of one of the sublattice potentials i.e. V 2 (r) → qV 2 (r). Investigations of detuned lattices is beyond the scope of this paper and we do not consider it here. We point out that this approach can be readily extended to lattices with more than two lattice sites per unit cell by adding an additional sublattice V 3 (r) function to potential function (2) that is defined similar to Eq. (3). For well-separated lattice sites, and when r is near r j or r k , the potential function (2) is approximated by the first few terms of its Taylor series and thus takes a paraboloid form, namely ∆n(r) ≈ |∆n| The orbital functions used to construct the ansatz in Eq. (11) are normalized Gaussians centered at the lattice sites and satisfy a linear Schrödinger equation with potential function of the form in Eq. (4).
To break time-reversal symmetry the two sublattices V 1 (r) and V 2 (r) are rotated according to the functions h 1 (z), h 2 (z), respectively, i.e.
The only restrictions we place upon the functions h 1 (z), h 2 (z) is that they are smooth. This is quite general. Since most researchers have studied periodic longitudinal variation, below we will assume h 1 (z), h 2 (z) are periodic; but they need not be so. In particular, we take h i (z) = R i (cos (Ω i z + χ i ) , sin (Ω i z + χ i )) , i = 1, 2 (5) where R i is the helix radius, Ω i is the angular frequency of the oscillation, and χ i a phase shift.
Next we move into a coordinate frame co-moving with with the V 1 (r, z) sublattice by performing the change of variable which (after the prime notation is dropped) gives the equation (6) Introducing the phase transformation for the pseudo-field defined by A(z) = −k 0 h ′ 1 (z) yields (after dropping the tilde notation) We non-dimensionalize by where l is the distance between nearest-neighbor lattice sites, z * = 2k 0 l 2 is the characteristic propagation distance, I * is the peak input beam intensity, nonlinearity coefficient σ = 2γk 0 l 2 I * and lattice amplitude V 2 0 = 2k 2 0 l 2 |∆n|/n 0 . Dropping the ′ notation, these rescalings give the dimensionless equation with the dimensionless potential In this paper we focus on periodic rotation in z. The most general form of the dimensionless rotation functions is taken to be where η i ≥ 0 is the ratio of the helix radius to the distance between adjacent lattice sites, ǫ −1 i is the helix rotation frequency, and χ i is a phase shift. The sign of γ is taken to be positive corresponding to self-focusing Kerr nonlinear media [24] e.g. fused silica. To simplify Eq. (8) the phase transformation ψ(r, z) = φ(r, z)e −ir·A(z) is introduced and gives The dimensionless pseudo-field is given by where κ = k 0 ℓR 1 Ω 1 = η 1 /(2ǫ 1 ). For each of the examples considered below we simply drop the subscript and call these parameters η, ǫ.
To simplify the analysis of Eq. (9) a tight-binding approximation is applied. From a physical point of view this assumption is justified by the fact that many photonic experiments are performed in strong lattice regimes where V 2 0 ≫ 1. In the deep lattice limit the scalar field φ(r, z) is well-approximated by coupled evanescent modes centered at the lattice sites. Rigorous studies of these types of approximations were carried out in [25].
The translational symmetry of the lattice motivates the following ansatz [26] where the Gaussian functions φ 1,v and φ 2,v satisfy the equations for the local paraboloid potential (see Appendix A for more details). The two non-simple lattices configurations we consider in detail are honeycomb (see Fig. 1) and staggered square (see Fig. 2). The perfect interpenetrating square lattice, as shown in Fig. 2, is actually simple, but becomes nonsimple when there are necessarily different sublattices; e.g. by introducing a phase offset in one of the sublattices. Other lattice configurations can be considered with the method we present. The lattice sites for the honeycomb lattice are related via the characteristic vectors where |d| = 1 is the distance between nearest neighbors. The staggered square lattice is defined by the basis vectors and has a unit distance between neighboring lattice sites as well.
We substitute ansatz (11) into governing equation (9) and multiply the resulting equation by φ * i,p (r, z), i = 1, 2, where * indicates complex conjugation. The resulting equation is integrated over the entire domain. For a static isotropic lattice the interaction strengths are proportional to exp −(V 0 |v| 2 )/(4σ) , where |v| is the distance between the different lattices sites. For a typical set of honeycomb lattice parameters (e.g. V 2 0 = 50, σ = 1/2) the nearest-neighbor coefficients are on the order of O(10 −2 ), while next-nearest neighbor interactions are considerably smaller at O(10 −5 ). In the case of the staggered squared lattice the next-nearest neighbor interaction is on the order of O(10 −4 ). For this reason we only consider self and nearest neighbor interactions and neglect all others.

A. Honeycomb Tight-binding Approximation
For the honeycomb lattice displayed in Fig. 1 we derive a tight-binding approximation that takes into account the interactions of a Gaussian orbital with itself and its three nearest neighbors. In order to present the simplest picture possible, only the dominant terms are retained. The tight-binding approximation for this system is given by Honeycomb : where with θ j (z) = (k + A(z)) · v j . Without loss of generality, we take k = 0 [23]. The definitions of the coefficients L j (z), R j (z), j = 1, 2, 3 are given in Appendix B 1. The z-dependent coupling coefficient, c j (z), is a function that decays exponentially with V 0 . The linear terms that are induced from the lattice potential and rotation effects are included in L j (z). For a deep lattice (V 0 ≫ 1), the magnitude of c j is small and the magnitude of L j large, however for the values of V 0 considered here, their product is numerically O(1). When the two sublattices have the same rotation, i.e. ∆h 21 = 0, this system reduces to that considered in [23] by rescaling Z = c 0 L 0 z, where all c j and L j are now constant and R j is zero (see Appendix A).
Here we focus on the evolution of edge modes along the zig-zag edge of a semi-infinite strip domain, whose left side is shown in Fig. 1. As the beam evolves down the waveguide it is assumed to be well-confined inside the lattice region and negligibly small outside. The boundary conditions chosen to model edge modes on the left zig-zag boundary are: The boundary conditions for edge modes on the right zig-zag edge mirror these. We consider solutions of the form such that

B. Staggered square Tight-binding Approximation
Next we give a tight-binding approximation describing beam propagation in the staggered square lattice shown in Fig. 2. Only the dominant terms resulting from interactions of a Gaussian orbital with itself and its four nearest neighbors are considered. We point out that the indexing given for the square lattice in Fig. 2 is different than that of the honeycomb lattice. Here the white lattice sites (i.e. the 'b' sites) are all located at even points, and black sites (i.e. the 'a' sites) are located at odd positions. This is due to the fact that the underlying lattice is rectangular. The coupled system describing this is Again, for simplicity we take k = 0. We focus on the edge modes propagating along the edge of a semi-infinite strip domain. For localized modes traveling on the left side of the domain we take the boundary conditions a 2m+1,2n+1 = 0 for n < 0 , b 2m,2n = 0 for n < 0 , a 2m+1,2n+1 → 0 as n → ∞ , b 2m,2n → 0 as n → ∞ .
The boundary conditions on the right side mirror these. We take solutions of the form for real ω which yield the following coupled system where

III. LINEAR FLOQUET BANDS AND EDGE STATE DYNAMICS
In this section we consider a low amplitude linear limit (i.e. |a mn | 2 , |b mn | 2 ≈ 0) of the full nonlinear systems given in the previous section. The dispersion relation α(ω) is computed numerically via Floquet theory [27]. The Floquet multipliers γ = e −iαT +i2πτ , τ ∈ Z are obtained from the eigenvalues of the fundamental matrix solution at z = T , where T is the period of the lattice. A fourth-order Runge-Kutta method is used to integrate. For all band diagrams below, 40 lattice sites (in both a n and b n ) are used. The Floquet exponent is calculated up to an additive constant by To be specific we will focus on five different rotation patterns among the sublattices: • same rotation, same phase • different radii, same phase • π-phase offset rotation • counter rotation • different frequency, same phase We point out that only h 2 (z) is adjusted to take into account the non-synchronized sublattice motion above.
The pseudo-field, A(z), is defined in Eq. (10). The physical parameters chosen in the simulations below are presented in Table I and reflect the experimental setup used in [5]. For these values the dimensionless parameters are ǫ ≈ 0.75/π and V 2 0 ≈ 45. Moreover, one unit in the dimensionless z is equal to 6.5 mm in physical units.

A. Honeycomb Floquet bands
In this section we explore the linear band structure induced by various rotation patterns for the semi-infinite honeycomb lattice shown in Fig. 1. The Floquet exponents defined in (27) are determined for the linear one-dimensional system given in Eqs. (20)- (21). The first Brillouin zone for several Floquet bands is shown in Fig. 3. Each band structure consists of bulk/extended (solid regions) and edge/localized (curves) modes. The parameters used in the simulations are given in Fig. 3.
The first case we consider is that of isotropic waveguides rotating in phase with each other (28) [see Fig. 3(a)]. We find this case supports unidirectional edge states. Typical eigenfunctions from both curve branches are displayed in Fig. 4. The eigenmode shown in Fig. 4(a) corresponds to the left zig-zag edge and has strictly negative group velocity i.e. α ′ (ω) < 0. The situation is reversed for edge modes on the right zig-zag edge, shown in Fig. 4(b), that have positive group velocity (α ′ (ω) > 0). Taken together, this indicates that the edges modes on a finite domain propagate along the boundary of the lattice in a counter-clockwise fashion, the same direction of the rotating waveguides. Moreover, the outer boundary modes [|b 0 | in Fig. 4(a)] are found to have approximately 20 times larger magnitude than the inner boundary mode [|a 1 | in Fig. 4(a)].
An asymptotic theory, given in Appendix C, yields (to leading order) left zig-zag edge modes a mn (z; ω 0 ) = 0 , b mn (z; ω 0 ) = C(z)b s n (ω 0 )e imω0 , (33) evaluated at the wavemode ω 0 , where C(z) = C 0 exp (−iǫαz) for constant C 0 andα(ω 0 ) is found explicitly in terms of integrals; it is given in Eq. (C10). The theory is valid in the regime where the waveguides are rapidly oscillating (|ǫ| ≪ 1) and describes edge modes in the central gap near α = 0; we refer to this gap as the Floquet center. Where applicable, the asymptotic solutions are compared against those computed numerically.
Next we consider when the waveguides are rotating in phase with each other (28) but are anisotropic (σ x = σ y ). The band structure in Fig. 3(b) has lost its topologically protected edge character and we expect backscatter. This observation comes from the fact that the slope of edge curves are no longer sign-definite in one Brillouin zone. We point out that this structure only occurs when σ x > σ y i.e. when the waveguides are elliptical with the major axis in the x-direction (also n-direction). In terms of the results presented in [23], this regime corresponds to ρ < 1/2.
We now examine when the two sublattices have the same phase, but one of the waveguides has a radius that is three-quarters or half the size of the other (29) : Fig. 3(cd). In the first case (R a = 3/4), the crossing modes observed in Fig. 3(a) have deformed into the non-crossing curves seen in Fig 3(c), which still have a sign-definite slope. The fundamental Floquet exponent, τ = 0, as well as the first non-fundamental Floquet bands, τ = ±1, defined in Eq. (27) are shown in Fig. 3(d). The point at which these two branches meet is α = ±π/T . We refer to this point as the Floquet edge. We see that there is a family of edge modes that live in this Floquet edge gap; Now let us consider a π−phase offset between the two sublattices (30). The band structure for this arrangement is shown in Figs. 3(e-f). Included is the fundamental Floquet band (τ = 0) as well as its periodic extensions, τ = ±1. Presentation of the different branches of the Floquet bands in Fig. 5 illustrates how one can understand the meeting/merging of the the inter-band curves near the Floquet edge. At sufficiantly small helix radii the band structure consists of bulk bands separated by gaps and exhibits no edge modes. As the radius of the lattice helix, measured by η, increases the gap between the adjacent Floquet bands begins to close until the bands "kiss" each other at a single point in Fig. 3(e). This point marks a transition from which any additional driving introduces a new family of topologically protected edge modes at the Floquet edge [see Fig. 3(f)]. We point out that these edge states are distinct from those at the Floquet center. Another feature that we observe in this case, unlike the case in Fig. 3(a), is that the edge states do not remain as ǫ → 0 (This is shown in Fig. 11).
The next case examined is that of counter rotation (31). The bands corresponding to this case are shown in Fig. 3(g). We find two flat bands, separated by a small gap, that correspond to non-traveling states. The opposing motion of the two sublattices evidently cancel each other out yielding zero net movement; i.e. essentially a stationary state. The final scenario is when one sublattice has twice the frequency of the other sublattice (29). The band structure, displayed in Fig. 3(h), exhibits a set of topologically protected edge modes. We find that this configuration may possess a weak instability[, namely, |Im α| = O(10 −3 ); i.e. it would be a weak exponentially growing/decaying mode that would only manifest itself after a long distance.  Fig. 3(a), (c-d) Fig. 3(b), and (e-f) Fig. 3(f).

B. Linear edge mode dynamics in honeycomb lattice
In this section we explore the dynamics of the linear edge states found in Fig. 3. Specifically, we numerically integrate the full tight-binding system (17)-(18) corre-sponding to solutions of the form given in Eq. (33). In all simulations we omit nonlinearity: σ = 0. The solutions are initialized at a chosen frequency, ω 0 , by where µ = 0.3 denotes the slowly varying envelope and b n (ω 0 ) is the numerically computed mode whose Floquet bands are shown in Fig. 3. For consistency, we normalize all edge modes so the two-norm is one, namely < b n , b n >= 1 via the discrete inner product < f n , g n >= n f * n g n . We take periodic boundary conditions in the m-direction and zero boundary conditions in the n-direction. For all z-evolutions shown in this paper we take N = 100 sites in the n-direction to ensure the mode has sufficiently decayed. The system is integrated using a fourth-order Runge-Kutta method.
Several typical evolutions are displayed in Fig. 6. In the left column we show the solution magnitude at the left boundary of the domain i.e. n = 0. In the right column the maximum magnitude for both a mn (z) (in blue) and b mn (z) (in black) over the entire domain are given. In Figs. 6(a-b) the propagation of a topologically protected edge mode, whose dispersion curves are given in Fig. 3(a), is shown. The traveling mode moves in the negative direction at a constant speed. At a point of inflection, α ′′ (ω) = 0, a method of stationary phase calculation shows that the edge mode decays like ∼ z −1/3 . Next in Figs. 6(c-d) a stationary edge state that is not topologically protected [corresponding to the bands in Fig. 3(b)] is shown. In contrast to the previous case, this state is diffracting and losing amplitude at a rapid pace. Asymptotically, when α ′′ (ω) = 0 the mode decays like ∼ z −1/2 . The next mode we consider is located near the Floquet edge in Fig. 3(f) and corresponds to π-offset rotation. Unlike the previous two cases, whose energy is always in the b m0 lattice sites, this solution [shown in Figs. 6(ef)] is observed to oscillate back-and-forth between the b m0 and a m1 lattice sites. Energy is regularly transferred back-and-forth each cycle of the lattice waveguides (here the period is T = 1.5). A similar evolution pattern (not shown here) was found to occur for the edge state mode corresponding to different sublattice radii whose Floquet bands are shown in Fig. 3(d).
Floquet theory does tell us something about this coupled-mode dynamic. Typically, the 1D honeycomb Eqs. (20)- (21) are integrated over the period [0, T ] to find the fundamental matrix solution. The corresponding edge mode eigenfunction [similar to that in Fig. 4(a)] has considerably more energy in the outer b n sites than those of a n . If we instead calculate the Floquet multipliers over the interval [T /2, 3T /2], then the energy is primarily concentrated in the a n sites, rather than b n . From this observation we infer that these solutions are not of the form given in Eq. (33), and instead resemble a n (z) = A(z)r n and b n (z) = B(z)r n for |r| < 1, n = 0, 1, 2, . . . , where A(z) and B(z) are T -periodic envelopes. In other words, the modes are truly coupled (a n = 0) and can not be assumed to be scalar (a n ≈ 0).

C. Staggered square Floquet bands
In this section the Floquet bands for the staggered square lattice in Fig. 2 are computed for the rotation patterns described earlier in equations (28)-(32). The Floquet exponent (27) is computed from 1D staggered square system (25)- (26). Several band structures are shown in Fig. 7.
The first lattice configuration we examine is when the two sublattices rotate in phase with each other given in equation (28). In this case, the system given in Eqs. (25)-(26) is degenerate and reduces to a single equation. In terms of the lattice structure, the lattice sites form a simple lattice. There are no edge modes at all, as indicated in Fig. 7(a); there are only extended bulk modes. One way to create a non-simple lattice configuration is to make the radius of one sublattice smaller than the other. The dispersion bands for this scenario are shown in Figs. 7(c) and 7(d) where the radii of one sublattice is 75% and 60%, respectively, the size of the other. As the radius disparity grows, localized edge state curves manifest themselves in the gap at the Floquet edge. This latter band structure is distinguished from all other cases in that for certain positive ω values there exist two different Floquet exponents: one topologically-protected near the Floquet edge and the other non-protected near the Floquet center. In most of the topological systems we investigated, either one or the other mode types exist at a fixed ω, but not both.
The two rotation patterns we consider next are counter rotation (see equation (31) and π-phase offset rotation (see equation (30). First, we examine the counter rotation case whose Floquet bands are shown in Fig. 7(b). If the radius of the sublattices are driven hard enough, then, similar to the honeycomb case, the bulk band splits apart and inside the central gap an edge mode corresponding to a nearly flat band is found to exist. Next we consider a π-offset rotation. We point out that this scenario was explored in [19,20]. Qualitatively, our results agree with theirs, namely there is a transition point in the topological structure of the bands that is displayed in Fig. 7(e). Below this transition point, in a weak rotation regime, the band structure takes a "trivial" form similar to that seen in Fig. 7(a). If, on the other hand, the spiral radius is increased beyond this threshold a set of topologically protected edge modes emerge and give the "nontrivial" band structure shown in Fig. 7(f). This transition point in band structure resembles the honeycomb cases given in Figs. 3(e-f).
The final case we consider is that of different frequency between the two sublattices (see equation (32). For a large enough helix radius, edge states are found to occur near the Floquet edge. Similar to the honeycomb case in Fig. 3(h), these states appear to possess a weak instability [|Im α| = O(10 −3 )] that could become relevant over very long propagation distances.
At this point we summarize the commonalities in the linear band structures between the honeycomb and staggered square lattices. In the case of same rotation pat- terns we find that edge modes are present if the underlying lattice is fundamentally non-simple (like honeycomb, or staggered with different sublattice size or frequency). Taking sufficiently different radii among the sublattices creates edge states at the Floquet edges, for both lattice types. Stationary (or nearly stationary) edge modes are generated in the case of counter-rotating sublattices. A family of edge modes is found to exist between different branches of Floquet bands (27) for π-offset rotation when the waveguide parameters exceed a certain thresh-old. When the two sublattices have different frequencies edge modes can be found (with a possible weak instability).  Fig. 7(b), (c-d) Fig. 7(d), and (e-f) Fig. 7(f).
In this section we propagate the linear modes associated with the dispersion curves shown in Fig. 7. The governing equations (22)-(23) are integrated using initial conditions of the form where µ = 0.3 and the decaying mode b 2n is numerically computed from the 1D system (25)- (26). As was the case in honeycomb lattice, we take N = 100 lattice sites in the n-direction for both b 2m,2n and a 2m+1,2n+1 . For all cases considered here we display the edge magnitude as well as the maximum magnitude over the domain as a function of z.
The first nontrivial case we consider is that of counter-rotation whose corresponding dispersion curves are shown in Fig. 7(b). The Floquet bands are nearly flat and so very little translation of the localized mode is expected. The edge mode evolution is shown in Figs. 8(ab) and indeed the mode is well-localized, has a small positive velocity and most of the energy is in the 'b' mode. Next we look at the case of different radii among the sublattices corresponding to the bands displayed in Fig. 7(d). The edge wave propagation and magnitude evolution is shown in Figs. 8(c-d). Similar to the honeycomb dynamics above, the edge modes whose Floquet exponents reside in the Floquet edge gap are observed to pivot back-and-forth among the b 2m,0 and a 2m+1,1 lattice sites. The energy is observed to oscillate with the same frequency as the lattice, T = 1.5. There is an additional submode excited in Fig. 8(c) that is moving left-to-right and corresponds to the mode located near α(π/10) ≈ 4 in Fig. 7(d). Moreover, this right moving mode is not a protected solution i.e. it will scatter at defects. It was numerically verified that the left-moving mode interacted with a defect like the topologically protected mode shown in the top row of Fig. 12; i.e. this mode travelled through the defect. On the other hand, the right-traveling mode was found to scatter off a defect like the mode displayed in the bottom row of Fig. 12. For the cases we studied it was unusual to have a system that simultaneously supports both topological and non-topological edge states.
In Figs. 8(e-f) which corresponds to Fig. 7(f) the πoffset edge profile is found to oscillate back-and-forth between the a and b sites with the same period as that of the helix (T = 1.5). This dynamics pattern resembles the π-offset honeycomb evolution shown in Fig. 6(e-f). Between the honeycomb and square evolutions we can see that similar lattice rotation patterns often yield similar edge mode evolutions, when the underlying band structures are similar.

IV. EDGE SOLITONS
In this section we explore the nonlinear (σ = 0) equations given in Sec. II. Using an asymptotic analysis, valid in the small ǫ regime, we find true edge soliton solutions, i.e. edge modes modulated by a slowly-varying envelope function that satisfies the NLS equation. Direct numerical simulations are used to validate our asymptotic results. On the other hand, this asymptotic theory does not describe every possible edge mode found in Sec. III. As we have remarked, modes located and the Floquet edge can disappear as ǫ → 0, and therefore are outside the scope of our small ǫ analysis.
To leading order, we find 2D zig-zag edge mode solu-tions to Eq. (17)-(18) of the form where the envelope function, C(y, z), satisfies the NLStype equation such thatα (j) 0 = d jα dω j ω=ω0 and the function C(y, z) varies slowly in y, i.e. |∂ y C| ≪ 1. The variable y is slowly-varying in the m-direction. The details of the analysis are given in Appendix C. In the direction perpendicular to the zig-zag edge, the solution in (36) decays like the stationary mode b s n = r n for |r| < 1 [the value of r(ω 0 ) is defined in Eq. (C6)]. This solution is derived under a narrow-band approximation, which assumes that only the frequencies, ω, near ω 0 make substantial contributions to the solution. Then by setting C(y, Z) =C(y, Z)e −iα0Z we obtain (to leading order) the NLS equation which for α ′′ 0 > 0 (focusing) has the traveling soliton solutioñ with µ ≥ 0. At points of inflection [like at ω = π/2 in Fig. 3(a)], α ′′ 0 = 0 and therefore the leading order dispersive term is now third-order and we get the higherorder NLS equation which does not support solitons. Next, we numerically verify these asymptotic results by direct numerical simulations performed on the full tightbinding system (17)- (18). Here we restrict our attention to the nonlinear modes found using our asymptotic theory (corresponding to red curves in Fig. 3). Analysis of the remaining modes is outside the scope of this paper.
To initialize the simulations we take functions of the form in Eq. (36). For the figures shown here the linear decaying mode, b s n , were computed numerically, but we did check that the asymptotic solution, b s n = r n , gave similar results. We normalize the edge mode two-norm so that < b n , b n >= 1. The envelope is initialized by the localized function for an amplitude, A, chosen to balance the leading-order dispersion term (either α ′′ 0 µ 2 or α ′′′ 0 µ 3 ) with the cubic nonlinearity, A 2 α nl .
Two edge mode evolutions are displayed in Fig. 9. In the first case, a traveling mode whose corresponding linear edge mode [located at ω 0 = π/2 in Fig. 3(a)] is at a point of inflection (α ′′ 0 = 0). This implies that the envelope is governed by the higher-order NLS Eq. (40) which does not support pure soliton modes. A closer comparison of the full discrete solution on the domain boundary, b m0 , and the envelope C(y, z) is shown in Fig. 10, where the initial and final solution profiles are shown in Figs. 10(a) and 10(b), respectively. The envelope captures the translation and a small dispersive tail which becomes more pronounced as z grows.
The next case we consider is when the envelope function satisfies the classic NLS equation (38). We choose the same lattice parameters that correspond to the linear Floquet bands shown in Fig. 3(b). At the point ω 0 = 3π/4 the value of α ′′ 0 is positive indicating that the NLS equation (38) is focusing and thus admits soliton modes. The evolution of such a nonlinear edge mode is shown in Fig. 9(b). The mode is stationary because α ′ 0 = 0. A closer look at the initial and final edge soliton profiles is given in Figs. 10(c-d). Comparing these two figures we see that the final magnitude is nearly the same as that which was initially injected. In both this case and previous one there is a small gap between the envelope and discrete mode peaks that may be attributed to some energy being transferred from the b mn mode into the a mn mode. We also remark that in the latter case it so happens that |b 0 | < 1, so there is also an initial gap at z = 0.
The asymptotic results shown here only cover edge modes located near the Floquet center. Moreover, the curves at the Floquet edges in Figs. 3 and 7 are the result of finite ǫ and disappear as ǫ → 0. This is highlighted in Fig. 11. In this case there are two distinct families of edge modes: one that bifurcates from the stationary mode [located near α = 0 in Fig. 11(a)] and another [located near α = ±π/T in Fig. 11(a)] separate one that arises when the helix period is not necessarily small in comparison to the characteristic propagation distance scale (not oscillating too fast). As ǫ approaches zero the edge modes at the Floquet edge are observed to disappear, leaving , and envelope (red curve), C(y, z). The parameters in panels (a-b) are the same as those in Fig. 9(a), and similarly for panels (c-d) and Fig. 9(b).
only those modes predicted by the asymptotic theory in Appendix C. Indeed, to understand a mode that oscillates like the linear edge state in Fig. 6(e) requires a true coupled-mode theory, i.e. not setting a mn = 0, to account for this back-and-forth energy transfer. To understand this finite ǫ edge modes is beyond the scope of this paper.  Fig. 3(f).

V. LATTICE DEFECTS
Lattice defects and imperfections are common in any real photonic lattice. Here we consider a scenario where many lattice sites are absent along the boundary, thereby introducing a boundary notch, or wall. The absence of any lattice sites along the boundary means the wave field is effectively zero there. To implement these boundary effects into the tight-binding model above we set a mn = b mn = 0 at the lattice defect locations. To track the evolution dynamics of the edge mode as it confronts the defect two quantities are monitored: the maximum intensities, |a mn (z)| 2 and |b mn (z)| 2 , and the 'participation' ratio The first quantity measures the edge mode peaks, while the latter gives a measure of the pulse width since (42) is proportional to it. We focus our attention on three illuminating examples. In the first situation we evolve a linear (σ = 0) topologically protected edge mode whose corresponding Floquet exponent is given in Fig. 3(a) at ω 0 = π/2. The evolution of the edge mode when it comes into contact with the lattice defect is shown in Fig. 12. In most nontopological systems significant scattering would be expected. Here, however, the edge mode tracks around the boundary notch and exhibits no backscattering. The intensity and participation ratio evolutions in Fig. 13 shed some light on the edge mode-defect dynamics. Along the defect boundaries perpendicular to the n = 0 axis the boundaries have armchair configurations. It is here that the energy begins to evenly distribute between the modes [see Fig. 13(a)]. In addition, the outer edge mode, b m0 , is found to spread to many sites (become wider) [see Fig. 13(b)] along these armchair boundaries. The outgoing intensity is observed to have nearly the same magnitude as the incoming value. The second case we consider is that of a nontopologically protected mode. The corresponding Floquet exponent is shown in Fig. 3(b) at ω 0 = π/2. Recall that this mode resulted from elongating the waveguide geometries in the direction perpendicular to the zig-zag boundary. The edge mode dynamics shown in Fig. 12 contrast from those in the previous case. Upon making contact with the defect boundary the edge mode scatters back in the opposite direction. From Fig. 13(c) the mode intensity is found to double at the defect corner, meanwhile Fig. 13(d) shows that the pulse width is cut in half (becomes more localized). Upon reflecting backwards significant dispersion in the edge state is observed.
Finally, we consider the interaction of a nonlinear edge mode with a defect. In particular, we focus on states we refer to as topologically protected edge solitons. These modes are characterized by linear Floquet bands with sign-invariant group velocity [see Fig. 3(a)] whose envelopes, to leading order, satisfy the NLS equation (38). In general, we expect these states to combine the robustness of solitons with the unidirectionality of a protected state. Consider dispersion bands like those shown in Fig. 3(a). At the point ω = 5π/8 ≈ 1.963 we have α ′′ 0 > 0 and |α ′′′ 0 | ≪ 1, yet α ′ 0 < 0. In this case the traveling envelope satisfies the focusing NLS equation and therefore it is a true soliton. One such case is displayed in Fig. 14. This nonlinear edge mode is observed to track around the defect just like the linear protected mode in Fig. 12. The outgoing soliton intensity is found to be nearly the same as the incoming magnitude. The maximum intensity and the participation ratio for this case are shown in Fig. 15. These values closely resemble the linear ones in Figs. 13 (a-b). One difference is that some loss is in the magnitude is observed in nonlinear case. This may be attributed to the fact that we are not in the pure NLS (38), but instead of a small amount of higher-order dispersion, since α ′′′ 0 = 0. This extra dispersion results in a small amount of radiation being emitted from the mode Decreasing the value of µ should improve the asymptotic NLS approximation.

VI. CONCLUSIONS
We have introduced a direct method for deriving tightbinding approximations of beam propagation in general 2D longitudinally driven waveguide arrays. As prototypes we examined periodically driven honeycomb and staggered square lattices. The new set of governing equations we derived allows us to find and examine unexplored edge modes with rather complicated sublattice rotation patterns such as: counter rotation, different radii, different frequency, different structure and π-phase offset. We computed the linear Floquet bands for these rotation configurations and also examined their evolutions.
An asymptotic theory was developed; it showed that the nonlinear edge modes behave as linear edge states modulated by a slowly-varying envelope that satisfies the 1D nonlinear Schrödinger equation. Ideal nonlinear edge modes combine the scatter-free topological properties of the linear problem with the robust nature of solitons.
Finally, topologically protected modes, both linear and nonlinear, were shown to possess unidirectionality (absence of backscatter), even when encountering strong lattice defects. In the absence of topological protection the mode simply reflects off a barrier. The approach for deriving tight-binding equations and analyzing photonic topological insulators allows the exploration of many new Floquet and non-Floquet lattice systems.
Here we provide some additional comments/details on deriving the honeycomb lattice tight-binding approximation in Sec. II. The derivation for the staggered square lattice follows in a similar manner. The ansatz in Eq. (11) is defined in terms of the two Gaussian orbital functions φ 1,v (r, z) and φ 2,v (r, z) centered at the points v and v+d+∆h 12 (z), respectively. These evanescent functions are taken to satisfy a linear (σ = 0) and local (r → r j , r k ) set of weakly-rotating (|A z | ≪ 1) equations where In particular, the normalized orbital functions used are where (A6) Substituting ansatz (11) into Eq. (9) for the Gaussian orbitals that satisfy Eqs.
where subscripts x and y denote theî and vector components, respectively. We next multiply this equation by φ 1,p (r) and φ 2,p (r, z) and integrate to get the 'a' and 'b' mode equations given in Eqs. (17) and (18), respectively. For deep lattices (V 0 ≫ 1) we may consider only the self and nearest-neighbor interactions to good approximation.
Following integration of the orbital functions, several simplifications to the remaining equations can be made. All self-interaction phase terms can be removed by the simple transformation a p+v (z) =ã p+v (z) exp(iΦ p+v (z)) and b p+v (z) =b p+v (z) exp(iΘ p+v (z)) where Next it is observed that |da p /dz| = |db p /dz| = O exp[−V 0 |d| 2 /(4σ)] . Since the off-diagonal derivative terms are of the form exp[−V 0 |d| 2 /(4σ)]da p /dz they are exponentially smaller than all other terms and may be neglected. We only consider the remaining dominant terms. These terms are numerically identified to be those with a coefficients of V 0 or V 2 0 (since we are in the deep lattice limit) and all other terms are dropped. Finally, to simplify working with these equations we relabel our grid in terms of a two-dimensional square lattice indexed as in Fig. 1.

Appendix B: Tight-binding Approximation Coefficients
In this section the tight-binding approximation coefficients for the honeycomb lattice given in Eqs. (17)- (18), and the staggered square lattice (22)- (23). The subscripts x and y denote theî and vector components, respectively. The nonlinearity coefficient in both cases is given by

Appendix C: Asymptotic Analysis
Details for deriving asymptotic solution (36) to the honeycomb lattice system (17)- (18) are presented here. This analysis generalizes the calculation performed in [23] to cover the more general non-synchronized rotation patterns discussed in this paper. The periodic functions ∆h 21 , ϕ and A are all assumed to depend only on the fast variable ζ = z/ǫ, where |ǫ| ≪ 1 and weak nonlinearity of σ = ǫσ is assumed. To begin, we take functions of the form given in Eq. (19) and then expand a n = To eliminate secularities these equations are rewritten as i ∂a i ∂b where κ 0 (ζ, ω) = e id·A+iϕ c 0 (L 0 − iR 0 ) , The terms f ± are expanded in series of ǫ and at the next order we get where we have rescaled Z = ǫz. Removing secularities, as we did in Eqs. (C3)-(C4), yields i db where Eq. (C4) gives a and α nl (ω) = < (b s n ) 2 , (b s n ) 2 > < b s n , b s n >σ , using the inner product < f n , g n >= n f * n g n . We now convert the spectral NLS equation (C9) into its corresponding spatial version. To begin, we expand the Floquet exponent in a Taylor series expansion around the central frequency, ω = ω 0 , bỹ ′′ 0 + · · · (C13) Next we take the inverse Fourier transform of this equation. When a narrow band is taken in ω this corresponds to a wide spatial profile. In other words, the spatial mode is slowly-varying along the zig-zag boundary, or in other words, it is nonzero at many lattice sites. Hence we take a Fourier transform that is continuous in space and replace the term (ω − ω 0 ) with the derivative −i∂ y . Doing this yields the nonlinear Schrödinger-type equation where derivatives in y are slowly-varying, i.e. |∂ y | ≪ 1.