7 verse index gradient entering in Eq.(34) may be also re- alized by applying a thermal gradient, or may describe lumped phase gradients [26] or an abrupt tilt of waveg- uide axis direction [3], in which cases E ( z ) shows a delta- like behavior. After introduction of a continuous function ψ ( r , z ) such that ψ ( r$_{n,m}$ , z ) = c$_{n,m}$ ( z ), one can read- ily check that the discrete diffraction equations (34) are equivalent to the following continuous Hamiltonian prob- lem tributions (such as the Gaussian beams) that propagate maintaining unchanged their functional shape (shape- invariant beams), and that field propagation may be sim- ply described by means of algebraic equations ruling out the evolution of some complex-valued beam parameters (such as the complex q -parameter for Gaussian beams; see, for instance, [27]). A natural question is whether one can similarly introduce shape-invariant discrete beams, i.e. field distributions that do not change their functional shape when propagating along curved waveguide arrays. As the problem of discrete diffraction in waveguide arrays with curved axis or transverselly-imposed index gradi- ents is analogous to the problem of one-dimensional or two-dimensional Bloch oscillations of electrons in peri- odic potentials with an applied electric field or of cold atoms in optical lattices, some results are already avail- able in the literature. In particular, in recent works [28– 30] an algebraic approach has been developed, capable of providing rather general results for wave packet center of mass evolution and wave packet spreading in certain lat- tice models. In this approach, after the introduction of a dynamical Lie algebra, an explicit form of the evolution operator is first derived, and then the expectation val- ues of operators are calculated in the Heisenberg picture. However, the question of existence of shape-invariant dis- crete beams and of their propagation in curved waveg- uide arrays does not seem to have been addressed yet. In this section, we present a generalization of Eqs.(12) and (17) describing the evolution of beam center of mass and beam width in curved waveguide arrays using the method of moments. Though similar results have been previ- ously published in Refs.[28–30] using an algebraic oper- ator approach, they are here re-derived for the sake of completeness using the method of moments, which does not require the explicit calculation of the evolution op- erator and the formulation of the problem in terms of a Lie algebra. In the subsequent section a family of shape- invariant discrete beams will be introduced, proving that their propagation in a generally-curved waveguide array is simply described by the evolution of a complex- q beam parameter, which plays an analogous role of e.g. the complex- q parameter of Gaussian beams propagating in paraxial continuous optical systems. Let us consider monochromatic light propagation in a two-dimensional waveguide array with a curved axis de- scribed by the parametric equations x = x$_{0}$ ( z ) and y = y$_{0}$ ( z ); the coupled mode equations describing light transfer among coupled waveguides in the single-band and tight-binding approximations are an extension of Eq.(1) to include fictitious transverse index gradients in- duced by waveguide curvature and read explicitly Application of the moment method to the one- dimensional Hamiltonian model (37) yields a set of closed coupled equations for the expectation values of operators x , θ , and of x $^{2}$, ρ and σ , where A. One-dimensional array where Δ$_{x}$ (Δ$_{y}$) is the coupling rate between adjacent horizontal (vertical) waveguides of the lattice [31]. H ( r , p ) = − 2Δ$_{x}$ cos( p$_{x}$ ) − 2Δ$_{y}$ cos( p$_{y}$ ) − E $_{x}$( t ) x − E $_{y}$( z ) y, (38) where Δ is the coupling rate between adjacent waveg- uides, and p = p$_{x}$ = − i∂$_{x}$ . The second model, which has been for instance considered in the experiment of Ref.[32], is described by the Hamiltonian where H$_{0}$ is the Hamiltonian of the homogeneous ar- ray defined by Eq.(3). The laws governing the evolution of beam center of mass and beam variance can be ob- tained by extending the method of moments described in Sec.II.B for the free diffraction problem. In general, the cascade of equations that one obtains by applying the Eherenfest equation (9) to 〈 r 〉 , 〈 x $^{2}$〉 , 〈 y $^{2}$〉 - and to the commutators found throughout the calculations- turns out to be unlimited for a general form of H$_{0}$ , and a closed set of equations are found solely for special forms of H$_{0}$ . Such a special circumstance is encountered in case of a one-dimensional waveguide array in the nearest- neighboring approximation, and in case of a rectangular- lattice waveguide array neglecting diagonal interactions. The first model corresponds to the Hamiltonian H ( x, p ) = − 2Δ cos p − E $_{x}$( z ) x, (37) H = H$_{0}$ ( p ) − E ( z ) · r , (36) i∂$_{z}$ψ ( r , z ) = H ( r , p ) ψ ( r , z ) (35) ( p = − i ∇ $_{r}$) with Hamiltonian where E ( t ) = E $_{x}$( t ) u$_{x}$ + E $_{y}$( t ) u$_{y}$ , E $_{x}$( z ) = − ( n$_{s}$/λ )¨ x$_{0}$ ( z ), E $_{y}$( z ) = − ( n$_{s}$/λ )¨ y$_{0}$ ( z ), n$_{s}$ is the refractive index of the waveguide substrate, and λ = λ/ (2 π ) is the reduced wavelength of light. It should be noticed that the trans- i ˙ c$_{n,m}$ = − ∑ l,r Δ$_{n}$$_{−}$$_{l,m}$$_{−}$$_{r}$ c$_{l,r}$ − E ( z ) · r$_{n,m}$ (34) θ = exp( ip ) (39) ρ = 1 2 { 1 − exp( − 2 ip ) } (40) σ = i { x exp( − ip ) + exp( − ip ) x } . (41)