9 fore, apart from an unimportant phase change, shape- invariant beams remain invariant for an arbitrary trans- verse translation on the lattice. Let us tentatively search for a solution to Eq.(72) of the form where J$_{n}$ is the Bessel function of first kind of order n , and α = α ( z ), σ = σ ( z ) are unknown functions which de- pend on propagation distance z , but not on lattice site n . Note that, as ∑ $_{n}$n | J$_{n}$ ( α ) | 2 = 0 and [ ∑ n | J$_{n}$ ( α ) | $^{2}$n$^{2}$] = α $^{2}$/ 2, the parameter α is related to the beam size w$_{x}$ [Eq.(6)] by the simple relation w$_{x}$ = $_{α/}$√$_{2, whereas}$ σ defines a transverse tilt of the beam ’phase front’. Sub- stitution of Eq.(74) into Eq.(72) and taking into account the identities of Bessel functions J$_{n}$$_{+1}$ ( α ) + J$_{n}$$_{−}$$_{1}$ ( α ) = (2 n/α ) J$_{n}$ ( α ) and J$_{n}$$_{−}$$_{1}$ ( α ) − J$_{n}$$_{+1}$ ( α ) = 2 J ′ $_{n}$( α ), one ob- tains that Eq.(74) is indeed a solution to Eq.(72) pro- vided that α and σ satisfy the coupled equations c$_{n}$ ( z ) = J$_{n}$ ( α ) exp( − iσn ) (74) ˙ α = − 2Δ sin σ (75) ˙ σ = − 2Δ α cos σ − f. (76) q$_{1}$$_{x}$ = i ∑ n,m (2 n − 1) c ∗ $_{n,m}$(0) c$_{n}$$_{−}$$_{1}$$_{,m}$ (0) (68) q$_{2}$$_{x}$ = n,m c ∗ $_{n,m}$(0) c$_{n}$$_{−}$$_{2}$$_{,m}$ (0) (69) q$_{1}$$_{y}$ = i ∑ n,m (2 m − 1) c ∗ $_{n,m}$(0) c$_{n,m}$$_{−}$$_{1}$ (0) (70) q$_{2}$$_{y}$ = n,m c ∗ $_{n,m}$(0) c$_{n,m}$$_{−}$$_{2}$ (0) . (71) IV. SHAPE-INVARIANT DISCRETE BEAMS where we have set 〈 y $^{2}$( z ) 〉 = 〈 y $^{2}$(0) 〉 + 2 | Ω$_{y}$ | 2 + 2Re { q$_{1}$$_{y}$ Ω$_{y}$ − q$_{2}$$_{y}$ Ω 2 y } (67) The existence of shape-invariant beams, i.e. families of field distributions that propagate without changing their functional shape, is well-known for paraxial propagation in Gaussian optics or in continuous lensguide media (see, for instance, [27]). Here we address the related problem of investigating the existence of shape-invariant discrete beams, i.e. field distributions that do not change their functional shape when propagating along waveguide ar- rays with arbitrarily curved optical axis. This is a rather challenging problem because no general method capable of constructing shape-invariant beams seems to be avail- able. However, for the simple waveguide array models considered in the previous section, a family of shape- invariant discrete beams can be introduced in a simple manner. Owing to their functional form, such beams are referred to as discrete Bessel beams. Owing to the functional form of c$_{n}$ , we will refer such shape-invariant beams to as discrete Bessel beams. Let us define a complex- q parameter for the discrete Bessel beam (74) according to q ( z ) = α ( z ) exp[ iσ ( z )] (77) dq dz = − 2 i Δ − if ( z ) q. (78) so that the modulus of the complex q parameter gives the beam spot size at propagation distance z , whereas its phase corresponds to the phase front gradient. From Eqs.(75) and (76) one readily obtains for the complex q parameter the following simple evolution equation Let us consider a one-dimensional waveguide array with an arbitrarily curved optical axis. In the tight- binding and nearest neighboring coupling approxima- tions, light propagation is described by the following set of coupled-mode equations A. Discrete Bessel beams in one-dimensional arrays The general solution to Eq.(78), for a given initial value q (0) at the z = 0 input plane, is given by q ( z ) = exp[ − iφ ( z )] { q (0) − 2 i ∫ z 0 dξ Δ exp[ iφ ( ξ )] } (79) φ ( z ) = ∫ z 0 dξf ( ξ ) . (80) where i ˙ c$_{n}$ = − Δ( c$_{n}$$_{+1}$ + c$_{n}$$_{−}$$_{1}$ ) − nf ( z ) c$_{n}$ (72) g$_{n}$ ( z ) = c$_{n}$$_{−}$$_{n}$$_{0}$ ( z ) exp { in$_{0}$ ∫ 0 dξf ( ξ ) } (73) The propagation of a discrete Bessel beam along a curved waveguide array is thus reduced to the prop- agation of its complex q parameter, which plays an analogous role of the complex- q parameter for Gaussian beams in lensguide media. The propagation law of the q parameter admits of a simple geometrical interpretation in the complex q plane. According to Eq.(78), for an infinitesimal propagation distance δz the change of q ( z ) is given by the superposition of the two paths AB and where f ( z ) describes the rate of transverse index gradi- ent induced by waveguide bending [13], lumped waveg- uide tilting [3] or locally imposed phase changes among adjacent waveguides [26] as discussed previously. Let us fist observe that, if c$_{n}$ ( z ) is a solution to Eq.(72) corre- sponding to a given initial field distribution c$_{n}$ (0), then for an arbitrary integer n$_{0}$ is the solution to Eq.(72) corresponding to the trans- lated initial field distribution g$_{n}$ (0) = c$_{n}$$_{−}$$_{n}$$_{0}$ (0). There-