Discrete Nonlinear Schrödinger Systems for Periodic Media with Nonlocal Nonlinearity: The Case of Nematic Liquid Crystals

: We study properties of an inﬁnite system of discrete nonlinear Schrödinger equations that is equivalent to a coupled Schrödinger-elliptic differential equation with periodic coefﬁcients. The differential equation was derived as a model for laser beam propagation in optical waveguide arrays in a nematic liquid crystal substrate and can be relevant to related systems with nonlocal nonlinearities. The inﬁnite system is obtained by expanding the relevant physical quantities in a Wannier function basis associated to a periodic Schrödinger operator appearing in the problem. We show that the model can describe stable beams, and we estimate the optical power at different length scales. The main result of the paper is the Hamiltonian structure of the inﬁnite system, assuming that the Wannier functions are real. We also give an explicit construction of real Wannier functions, and examine translation invariance properties of the linear part of the system in the Wannier basis.


Introduction
We study properties of an infinite system of discrete nonlinear Schrödinger (DNLS) equations that is equivalent to a coupled Schrödinger-elliptic system of partial differential equations with periodic coefficients. The system was derived in [1] as a model for the propagation of laser light in nematic liquid crystal substrates with a periodic structure in one of the directions normal to the optical axis. The model was originally motivated by experimental studies of such waveguide systems [2][3][4] and leads to extensions of a nonlocal DNLS equation of Fratalocchi and Assanto [5,6].
The Fratalocchi-Assanto equation has a nonlocal nonlinearity that leads to new effects when compared to the cubic power DNLS model studied commonly in photonics and atomic physics [7]. These effects include non-monotonic amplitude profiles of static (breather) solutions, additional internal modes in the linearization around breathers [8,9], and enhanced mobility of traveling localized solutions [10]. On the other hand, the mathematical justification of the Fratalocchi-Assanto model, in particular the question of how well it approximates the partial differential equations with periodic coefficients used to describe the underlying physics, is less studied. The present paper is a step in studying this problem.
Schrödinger-elliptic systems of differential equations with a similar nonlocal structure in the nonlinear term arise in a variety of contexts. Examples from physics include Bose-Einstein condensates [11], thermal media [12], and matter-wave microwave systems [13]. The recent review [14] includes further examples describing laser beams in liquid crystals [5,[15][16][17]. A related area of application of such models concerns thermo-optical interactions induced by beams in liquid crystals [18][19][20]. The combination of nonlocal nonlinearity and spatial periodicity or more general inhomogeneity, and the analysis of relevant equations is therefore a problem of wider interest.
The Schrödinger-elliptic differential equation we study describes the coupling of the laser field amplitude to the nematic crystal director angle. The derivation uses approximations of the coupled Oseen-Frank-Maxwell system for a linearly polarized beam [1]. The periodicity of the medium in the direction transverse to the propagation of the laser beam leads to an elliptic (Poisson-like) equation with periodic coefficients. Our approach is to expand the laser field and director angle in a Wannier function basis. The system is subsequently written as an infinite system of coupled DNLS equations for the Wannier mode amplitudes. The Wannier functions we use are defined in terms of a periodic Schrödinger operator appearing in the elliptic equation [1]. Note that the Wannier functions are integer translates of an infinite set of localized functions with an increasing degree of oscillation. Thus Wannier mode amplitudes give information on both the location and the spatial scale of images. Wannier and the related Bloch functions are a standard tool in the analysis of periodic Schrödinger operators [21][22][23][24], and related linear problems in theoretical physics, e.g., solid state physics [25].
Wannier functions are increasingly used in the study of nonlinear waves in inhomogeneous media. The use of Wannier functions for deriving discrete Schrödinger equations for nonlinear wave systems with periodic coefficients was first proposed in [26] for the periodic Gross-Pitaevski equation (NLS with periodic potential). The Wannier expansion has been used to justify the approximation of this equation by the DNLS equation in the the tight binding approximation limit for the potential term in [23,27]. Related systems where the theory applies are described in [11,12]. In the present problem the Wannier basis leads to a heuristic derivation of the model of [6] and also allows us to derive more general DNLS-type equations and systems that include additional inter-and intra-band Wannier mode interactions [1]. However, the Wannier approach does not immediately justify truncation to the lowest band because the linear part does not have the band gaps assumed in [23,26,27]. Thus the question of justifying the derivation of finite systems of (possibly a few) DNLS equations from the infinite system requires some additional analysis, and also motivates a better understanding of the structure of the infinite system.
A first result of the paper is an outline of the global existence theory, that is the boundedness of a suitable norm of the solutions. This type of result is a mathematical way to describe the absence of catastrophic self-focusing (beam collapse) and the possibility of stable localized beams [28]. The result also implies an estimate for the energy (optical power) at different length scales and provides a heuristic justification of truncation to a finite number of DNLS systems, corresponding to Wannier modes of the first bands.
The main result of the paper is a proof that the infinite system resulting from the Wannier basis expansion is a Hamiltonian system. This fact implies the Hamiltonian structure of the finite band truncations and can useful in analyzing discrete soliton structures, using for instance methods from [8,9]. The proof assumes that the Wannier functions are real, and we subsequently give examples of an explicit construction of real Wannier functions in terms of explicit Bloch functions.
We also examine some features of the linear part of the problem, in particular we show that it is diagonalized by the trigonometric functions. This observation implies that the dispersion relation and the coupling between the modes can be computed with relative ease, and that the linear part of the problem is homogenized in the Wannier basis, i.e., is effectively a translation invariant [29,30]. This latter property is an additional motivation for further developing Bloch-Wannier analysis in nonlinear wave equations.
The paper is organized as follows. In Section 2 we outline the global existence theory for the coupled Schrödinger-elliptic system and show that the system in the Wannier basis is Hamiltonian. In Section 3 we discuss the construction of real Wannier functions. We also discuss translation invariance properties of the linear part of the system. In Section 4 we discuss some questions for further work.

Hamiltonian Structure of Periodic Nematicon Equations
We consider the system of equations ("nematicon equations") with α, β, g 2 positive constants, and V b−periodic and positive. The complex amplitude u describes the electric field amplitude of a linearly polarized laser beam through a nematic liquid crystal sample, while ψ describes the director angle deviation of the liquid crystal due to the laser beam. The geometry of the problem is indicated in Figures 1 and 2, see also [3,6]. In Figure 1 we show a vertical direction x, and the laser beam propagation axis z. The y−axis is perpendicular to the plane of the figure. The laser beam electric field is polarized along the x−axis, while the the angle ψ is on the x, z plane. The device (medium) is periodic along the y−axis. The periodicity can be imposed by an external electric field that is also along the x−axis, see Figure 2. We also simplify the problem mathematically by ignoring the dependence of u and ψ in x. Boundary effects in the directions transverse to the beam are also ignored. Equations (1) and (2) were derived in [1] from Maxwell's equations coupled to the Oseen-Frank equations for the director field [16,31]. Schrödinger operator −∂ 2 y + V. Similar equations with constant coefficients have been studied widely in the context of optical solitons in liquid crystals ("nematicons") and other nonlocal media [5,[14][15][16][17]. In model (1), (2) the transverse periodicity of the medium is captured by the b−periodic function V, and our study involves the analysis of the periodic Schrödinger operator −∂ 2 + V. More detailed models [1] involve more complicated operators with periodic coefficients in the second equation. An example is the operator considered in [3]. The simplification used here captures the fact that the periodicity of the medium appears in the nonlinear term of the beam Equation (1). Z d Figure 2. The periodicity of the device (or medium) in the horizontal direction (y−axis) is due to a vertical electric field applied externally. The parallel stripes represent capacitors that apply a voltage that is uniform along the beam propagation direction (z−axis), and periodic along the y−axis, see [3,6]. The dependence on the vertical direction is not included in model (1), (2). Horizontal boundaries are also ignored.
The local and global existence theory for the initial value problem of system (1), (2) follows from standard arguments and similar to the one in 2-D in [31][32][33]. This theory implies that the solution avoids catastrophic nonlinear collapse in finite length, see [28]. This is an important feature of nonlinear beam propagation in nematic liquid crystals and related nonlocal media, and is a prerequisite for the existence of stable nonlinearly focused beams [17], see [33] for mathematical aspects.
The main ingredient of the global existence theory is the conservation of the Hamiltonian of the system (1), (2) and of the (optical) power P = R |u| 2 , the squared L 2 −norm of u. We use the notation We can use these two conserved quantities to show the boundedness of some simpler quantities. By the Cauchy-Schwarz inequality we have and using the boundedness of G in L 2 we obtain We use that for all y ∈ R using the Cauchy-Schwarz inequality. We then have and by (5) we bound the quartic part of the Hamiltonian as Then with P the power, a constant. The conservation of H and (10) imply that ||u y || L 2 must remain bounded for all z ∈ R. Also, R V|u| 2 ≤ ||V|| L ∞ P, thus we have a bound for all z ∈ R, with C 0 depending on H and P at z = 0. We now consider an equivalent discrete system using expansions in Wannier functions. We also examine some consequences of the Hamiltonian structure of system (1), (2) and of the bound (11).
We start by defining the Wannier functions associated to the Schrödinger operator −∂ 2 y + V(y), with V b-periodic, see [21,22]. Bounded solutions φ n,k (Bloch functions) and eigenvalues E n,k of the periodic Schrödinger equation satisfy where φ n,k (y) = v n,k (y)e iky , with v n,k (y + b) = v n,k (y), (13) for all y ∈ R, n ∈ N, k ∈ R. By V > 0 we have E n,k > 0, furthermore for all n ∈ N, k, y ∈ R. Then we can consider k in any interval of length 2π/b. The index n is referred to as band index (or number). For any fixed k in an interval of length 2π/b, E n,k is the n − th largest eigenvalue of (12) with boundary conditions φ n,k (y + b) = e ikb φ n,k (y), implied by (13). Also, by (12), (13), the b−periodic functions v n,k satisfy thus for any k fixed, n labels the eigenvalues E n,k in an increasing order. This equation can be also be used to compute the v n,k , E n,k numerically for each k ∈ [0, 2π/b). For n ∈ N, y ∈ R, we consider the Fourier coefficients of φ n,· (y), and we also have the inversion formula The set of functions w m n : R → C, n ∈ N, m ∈ Z defined by (16) are known as Wannier functions [21,22,25]. Note that the Bloch functions are not unique. One of the basic results is that we can define the Bloch functions so that the Wannier functions form an orthonormal basis for L 2 (R; C) [23,24]. We discuss the construction of Bloch and Wannier functions in the next section.
Another property of the Wannier functions is that, by (16), (13), Thus, fixing n, the function w m n is a translation of the function w 0 n by mb. We use expansions of ψ and u in Wannier functions w m n as By the orthonormality of the Wannier basis, the coefficients c n,m , u n,m are obtained from the physical quantities ψ, u by The Wannier functions and the integrals must be evaluated numerically (or approximately).
Note that the definition w 0 n and the regularity of the φ n,k in k can also lead to strong localization of the w 0 n in y, see [23,24,34] and the discussion of the next section. The decay of w m n is more pronounced for larger oscillation V and for the first n. Numerical examples are shown in [1]. For rapidly decaying Wannier functions, the decay of the coefficients of c n,m , u n,m in m reflect the decay of the spatial profile of ψ, u respectively.
We can also use the orthonormality of the Bloch and Wannier functions to derive a bound on the optical power of each energy band. Let and (12) we have that H 2 of (11) satisfies Let ε n = min k∈[0,2π/b) E n,k . We have ε n > 0, ∀n ∈ N. Then By the orthornormality of the Bloch and Wannier functions for the n−th band and (20), (22) Combining with (25), (11) we have for all n ∈ N. For large n we have ε n ∼ n 2 , more precisely, there exist c, C > 0 such that cn 2 ≤ ε n ≤ Cn 2 , ∀n ∈ N. We discuss this estimate in the next section. Therefore for all n ∈ N. This bound gives us the optical power in the higher band components, e.g., we can estimate number of modes needed to have a given high percentage of the power in the lowest band modes. This is a heuristic justification of using a finite system where n ∈ {1, . . . , N}, i.e., a truncation to the Wannier modes of a finite, possibly large, set of bands. Note that (29) does not give us however an estimate for the difference between solutions of the full and truncated systems. This question will be examined in future work.
Equations (1) and (2) in the Wannier basis, see [1], are where with System (30) was obtained in [1], and we describe the steps in the Appendix A. To show that it is a Hamiltonian system we compare (30) to Hamilton's equations with the Hamiltonian H of (3) expressed in the Wannier basis.
By (11) and by (33) We show that Hamilton's equations for (36) coincide with (30), provided the Wannier basis functions are real.
We first see that the symmetry of G implies that the coefficients G m,m n,n of (33) satisfy G m,m n,n = (G m ,m n ,n ) * , for all n, n ∈ N, m, m ∈ Z. We use the double index notation j = (n, m), i.e., (19), (20) are written as with summation over j = (n, m) ∈ N × Z. Then G m,m n,n = G j,j with j = (n, m), j = (n , m ). Let g j,j = R (Gw j )w * j . We will show that G j,j = g j,j and that the symmetry of G implies g j,j = g * j ,j . We write (2) as ψ = Gv, v = |u| 2 , and by (38), v = ∑ j v j w j we have Then or By (33) we have g j,j = G j,j , ∀j, j ∈ N × Z. Symmetry of the real bounded operator G with respect to standard L 2 inner product implies We now examine Hamilton's equations. We write (36) as Hamilton's equation is We have thus recovered the linear part of (30).
For the nonlinear part we have Let f = w j 1 w * j 2 , g = w * j 3 w l . We omit the dependence of f , g on the indices for simplicity. Also let f = ∑ k f k w k , g = ∑ k g k w k . By (37) and so that by symmetry of G, By (47), (49) we then have Λ j 3 ,l,j 1 ,j 2 = Λ * j 1 ,j 2 ,j 3 ,l .
Clearly, the above hold for any double index j 1 , j 2 , j 3 , l ∈ N × Z. If the Wannier functions are real, the coefficients G k 1 ,k 2 and Λ * j 1 ,j 2 ,j 3 ,j 4 are real. By (50), (46) yields the nonlinear part of (30). This concludes the argument.
We remark that the Hamiltonian structure of (30) easily implies the Hamiltonian structure of finite band truncations of the (30). The same applies to truncations where we consider a finite set of sites m. It suffices to restrict the summations in (36) to a finite range of n, m, also setting modes outside the desired index range to zero.
Also the Hamiltonian of (36) is invariant under the global phase change u n,m → e iφ u n,m , for arbitrary real φ and all n, m. This fact justifies the terminology coupled DNLS for (30).
As seen in [1], the Wannier expansion leads to a natural extension of the Fratalocchi-Assanto model [6]. The coupled mode approach of [6] can be also extended to describe more degrees of freedom per site [7]. Generally, mode expansions have additional structure when they arise from the solution of some spectral problem. This is the case for Bloch and Wannier functions. This additional structure however requires substantial computational effort, e.g. we need to compute Bloch and Wannier functions and evaluate Wannier overlap integrals. We discuss some of the relevant issues in the next section. We emphasize however that the general structural features of the equations, e.g., Hamiltonian structure, symmetries, form of mode interaction terms, are key. Heuristic simplifications that preserve these features can yield useful models. It is also seems important to be able to justify truncations to a small number of bands. We have at the moment only a partial justification for such truncations, relying on the rate of decay of the power in the higher bands (29).

Real Wannier Functions and Dispersive Properties
The Hamiltonion structure of the infinite system for the Wannier coefficients (30) assumed real functions. In this section we describe the construction of real Wannier functions using explicit constructions of the Bloch functions. We also observe that the linear part of (30) will in general couple modes from different bands. This is a main difference between our system and the equations considered in [11,12,26,27]. We show that we can still however diagonalize the system using trigonometric functions. In that sense, the Wannier-Bloch analysis leads to a homogenized , i.e., effectively translation invariant, linear part, see [29,30].
To construct the Wannier functions we examine the Schrödinger equation with V nonegative and b-periodic as a second order ODE with a real parameter E ∈ R. We assume that V is also piecewise Lipschitz. The spectrum of −∂ 2 y + V is given by the set of real E for which all solutions of (51) are bounded, see e.g., [24,35]. The equations for the real and imaginary parts of Ψ decouple, and all complex valued solutions are of the form AΨ 1 (.; E) + BΨ 2 (.; E), A, B ∈ C, where Ψ 1 (.; E), Ψ 2 (.; E) are any two linearly independent real solutions.
Let Ψ(y; E) = AΨ 1 (y; E) + BΨ 2 (y; E) for A, B complex, then (54), (55) at y = 0 and (52), (53) imply the system Then λ must be an eigenvalue of the matrix M defined by If λ is an eigenvalue of M, and λ = e ikb , then k and E are related to where see [21]. The dependence of µ on the choice of Ψ 1 , Ψ 2 , i.e., the initial conditions (52), (53) is supressed from the notation.
A band is a maximal connected interval of R where µ(E) is monotone and satisfies |µ(E)| ≤ 1. By the above, there is an infinite number of bands B n , n ∈ N, and a natural way to enumerate them so that points n < n , E ∈ B n , E implies E n ≤ E n , with equality holding for n = n + 1 and E = maxB n , E = minB n . Every E satisfying (60) must belong to exactly one band. Also, for every n ∈ N, and k ∈ [0, π/b], there exists a unique E ∈ B n satisfying (60). We denote such E ∈ B n by E k,n . The large E behavior of µ(E) also implies that there exist c, C > 0 such that ε n = minB n satisfies cn 2 ≤ ε n ≤ Cn 2 , ∀n ∈ N.
The solutions of (60) can be then parametrized as E n,k , k ∈ [0, π/b], n ∈ N, and B n = [E n,π/b , E n,0 ]. Also, we let E(−k, n) = E k,n and extend E n,k to k ∈ R by 2π/b−periodicity. This notation is consistent with (60). For V non-negative all bands belong to R + . By the implicit function theorem, given n ∈ N, E n,k is real analytic for k ∈ (0, π/b), and is continuous in [0, π/b]. The even and 2π/b−periodic extension of E n,k to real k is continuous in R, and real analytic at all points outside the lattice Z π b , for all n ∈ N. Regularity of E n,k at points Z π b for given n follows under gap conditions for the edges of the band B n . Consider now E = E n,k as above a solution of (60) for some n ∈ N, k ∈ [0, π/b], and the corresponding real solutions Ψ 1 (·; E), Ψ 2 (·; E). Solving (56), (57) we have and we obtain with λ = e ikb , E = E n,k , n ∈ N, k ∈ (0, π/b). The expression can be extended to the endpoints k = 0, π/b, under conditions we discuss below. Also, the complex coefficient A is free, e.g., it can be also chosen to normalize Ψ(y; E). In general it may depend on n, k, and we write A = A n,k . Denote Ψ(·; E) = Ψ(·; E n,k ) by Ψ n,k , n ∈ N, k ∈ [0, π/b]. Clearly Ψ * n,k is also a solution of (51). We then let The functions φ n,k are extended by 2π/b−periodicity for k ∈ R \ Z π b and are Bloch functions, see (12), (13), (14).
We check that the corresponding Wannier functions are real. By (18), it suffices to show that the w 0 n are real. By (16) and by (63) ∀y ∈ R. It is assumed that the last integral is well defined. We now give a condition that makes the above construction well defined, leading to w 0 n ∈ L 2 (R, R), for all n ∈ N. In particular, assume that the limits lim k→0 + dE n,k dk , lim k→π/b − dE n,k dk exist (and are finite). Then the corresponding limits of the fraction Ψ 2 (b; E)/(λ − Ψ 1 (b; E)) also exist, and Ψ n,k given by the right hand side of (62), with A = 1, is well defined for all n ∈ N, k ∈ [0, π/b]. We further see that the Ψ n,· (·) are continuous in [0, π/b] × R, ∀n ∈ N. Defining φ n,k by (63), the integrals in (64), (65) are finite for all y ∈ R. It follows that φ n,· (y) ∈ L 2 ([−π/b, π/b]; C), for all n ∈ N, y ∈ R. Then by Percival and (18) we see that i.e., finite, for all n. Thus the Wannier functions constructed this way are square-summable. Normalized Wannier functions are obtained by choosing a suitable coefficient A n,k = A n for each n ∈ N.
Note that the condition we used is always satisfied if both E n,0 , E n,π/b belong to the boundary of the spectrum for some n. In such band gap situations, E n,· is real analytic in R, and we obtain an exponentially decaying Wannier function w 0 n [23,24]. The above suggest that several qualitative features of the Wannnier functions can be deduced by theoretical arguments. The main input is information on the energy band structure. This information is obtained by solving (59) numerically. The function µ(E) must be computed numerically from (60). The functions Ψ 1 (y; E), Ψ 2 (y; E) are computed by numerical integration of (51) in the interval y ∈ [0, b] for different values of E, using the initial conditions (52), (53) respectively. Explicit expressions for the Ψ 1 (b; E), Ψ 2 (b; E) are known for a piece-wise constant potential V with two steps, see [1,23,36], but are cumbersome in the general case. This calculation also yields E n,k , k ∈ [−π/b, π/b], for the lowest n, numerical plots can be found in several sources, see e.g., [1,26].
Wannier functions are obtained numerically from (62)-(64), see [1,23,26] for some examples indicating the decay of the Wannier functions (for small n) as V is varied. The evaluation of mode interaction coefficients (33), (34) also uses quadrature, see [1]. The main difficulty here is the large number of coefficients, and the combinatorial nature of their enumeration. At this stage we need some efficient cut-off criteria, and we typically opt for some heuristic truncation to a few mode interactions, e.g., with a few nearest neighbors, justified by the decay of Wannier functions. This part of the analysis is still not as developed.
In the case where we consider truncation to the first band modes, the main question is distinguishing between a power (on-site-only) and a nonlocal nonlinearity [6], see [1] for some results. As we already mentioned the two models have properties the can distinguish them [8][9][10]. The possibility of long range linear mode interactions, e.g., as in [37], was also considered. It would be desirable to have a similar study for a model with two or three bands, inter-band mode interactions could be a more important feature of the problem.
We now examine the linear part of the nematicon system (1), (2). Generally the Wannier modes of different bands interact, and we want to examine the effect of these interactions for finite band truncations of the general discrete system (A8).
In what follows we will consider expansions in real Wannier functions and use the Hamiltonian structure of the linear systems. The linear coefficients of (31), (35) Let D n 1 ,n 2 (m) = D 0,m n 1 ,n 2 , then by (66) we have the symmetries D m 1 ,m 2 n 1 ,n 2 = D n 1 ,n 2 (m 1 − m 2 ) = D n 2 ,n 1 (m 2 − m 1 ) = D m 2 ,m 1 n 2 ,n 1 , for all n 1 , n 2 ∈ N, m 1 , m 2 ∈ Z. Linear interaction coefficients for n 1 = n 2 = n depend on |m 1 − m 2 |. In general, the linear interactions between bands n 1 = n 2 do not vanish. The linear coefficients D m 1 ,m 2 n 1 ,n 2 are contrasted to those of the periodic or perturbed periodic Schrödinger equation withṼ a perturbation of the b−periodic potential V used to define the Wannier functions. ForṼ ≡ 0 the Hamiltonian is and we use expansion in the Bloch functions and (12), (13) and the definition of the Wannier functions to compute withÊ Also, E n,k real and even in k, impliesÊ n (−m) =Ê * n (m) =Ê n (m), for all n ∈ N, m ∈ Z. The interaction between different bands therefore vanishes.
The effect of the perturbationṼ is described adding to the Hamiltonian the part We have hṼ = The coefficients will in general couple modes from different bands. In the case wherẽ V is also b−periodic theṼ m 1 ,m 2 n 1 ,n 2 have the symmetries that are similar to the ones in (67), i.e., V m 1 ,m 2 n 1 ,n 2 =Ṽ m 2 ,m 1 n 2 ,n 1 =Ṽ 0,m 1 −m 2 n 1 ,n 2 and will also couple modes from different bands. In the case where V +Ṽ is periodic a new set of Wannier functions may be defined so that the new bands decouple.
We remark that the Hamiltonian of the linear part of is denoted by h 2 , see (43). Clearly, h 2 = h 2,V − h V , with the notation of (70), (73). Thus the coefficients D m 1 ,m 2 n 1 ,n 2 of (66) can be expressed in terms of theÊ n (m), V m 1 ,m 2 n 1 ,n 2 of (71), (74), as To find the dispersion relation we look for solutions u j,m = A j e −i(mbκ+ω κ z) , j ∈ {1, . . . , N}, m ∈ Z, κ ∈ R. Then (76) becomes Thus ω κ is an eigenvalue of the matrix T κ defined by For instance, truncation up to the second band yields By the symmetries (67) T κ is Hermitian, ∀κ ∈ R, so that the eigenvalues ω κ are real. Note that the solutions u j,m and T κ are 2π−periodic in κ so that we may consider only κ ∈ [0, 2π). Varying κ ∈ [0, 2π) for each of the N eigenvalues ω κ,j , j = 1, . . . , N of T κ will produce N intervals.
We finally note that substitution of (31), (67) into (78), and use of (17) leads to somewhat simpler expressions that involve the Bloch functions The linear part of the evolution equation is therefore diagonalized by plane wave (trigonometric) solutions, and is thus effectively translation invariant in the Wannier basis. This is expected as the periodicity of the medium is absorbed in the nonlinear term. The range of the ω κ,j , κ ∈ [0, 2π), j = 1, . . . , N, will produce N intervals that must be computed numerically. These intervals may overlap or have gaps, although in the limit N → ∞ we expect that their union is the positive real axis, i.e., the spectrum of −∂ 2 y . The computation of these intervals involves the computation of linear mode interaction coefficients and will be considered in future work.

Discussion
We have examined some properties of a coupled Schrödinger-elliptic system modeling optical waveguide arrays in a nematic liquid crystal substrate. The system is studied by first passing to an equivalent infinite system of discrete describing the interaction of Wannier mode amplitudes of the relevant physical quantities [1]. The Wannier basis functions are defined in terms of a periodic Schrödinger operator appearing in the system and must be computed numerically. The Wannier mode amplitudes are related to the observed quantities by quadrature formulas (21) that may be evaluated numerically or using approximations. The Wannier basis approach leads to the derivation of systems of discrete nonlinear Schrödinger (DNLS) equations by truncation to the mode amplitudes of the first bands of the periodic Schrödinger operator. We have shown that these systems are Hamiltonian. The proof uses the reality of Wannier functions, and we have also described an explicit construction of real Wannier functions. Finally we show that the linear part of the system of discrete equations is diagonalized by (trigonometric) plane waves.
Bloch-Wannier analysis is a classical subject in theoretical physics, with well-known applications in classical and quantum mechanics [24,25]. Recently it is increasingly applied to the study of nonlinear waves in inhomogeneous media [23,26], and in the theory of homogenization [29,30]. The paper considers a problem in nonlinear waves where the periodicity appears in the nonlinear interaction. The system includes a natural operator that allows the use of Bloch-Wannier analysis, but differs from the more commonly stud-ied systems [7] where the linear part of the beam propagation equation is the periodic Schrödinger operator. In that case the band gap structure of the periodic Schrödinger operator allows us to justify the truncation to a model for the lowest band modes through a non-resonance argument [23]. In the present case the the most likely justification of DNLS models will involve truncation to a finite (possibly small) system of DNLS equations. This fact is also suggested by experimental results of [3]. Note that [3] consider a more complicated periodic operator in the director field equation that also includes the second (vertical) direction of the plane perpendicular to the optical axis. The present paper uses a simplification of the director field equation and takes advantage of the more developed Bloch analysis of the periodic Schrödinger operator. The idea is that the nonlinear effect of the transverse periodicity is already present in the simpler version. In addition, the more tractable Wannier-Bloch analysis of the simpler problem has allowed us to derive multiband DNLS systems and to analyze their structure. This a first towards further analysis of such systems, e.g., along the lines of earlier studies of the Fratalocchi-Assanto model [8,9]. While the paper considers a problem arising from nonlinear optics in liquid crystals, the combination of nonlocal nonlinearity and periodicity we consider may appear in other areas where similar systems are studied [14,15].
The paper shows that the Wannier basis expansion is an effective tool for elucidating the structural features of simplified DNLS equations. We also saw that the Wannier-Bloch approach requires the numerical computation of several intermediate quantities.
For instance, the computation of Wannier functions uses numerically computed Bloch functions and numerical integration over Bloch functions, while the nonlinearity of DNLS systems, as well as the linear interaction between bands also involves the evaluation of Wannier overlap integrals. Possible simplifications may arise for some limits of V [1], and there are general ideas such a eliminating non-resonant interactions [23,26] that can be examined further in this problem. but possibly less practical. Our view is that a good understanding of the structure of the Wannier coupled mode systems may allow us to analyze their dynamical properties without computing everything. It is possible that heuristic simplifications lead to models that capture significant features of the dynamics, and that a better theoretical understanding can bridge the gap between the simplified and fuller models. These questions will be addressed in future work. We then multiply (A2) by e imk b and integrate both sides over k ∈ [−π/b, π/b]. Interchanging the order of integration in k, y, and using the definition of Wannier functions we obtain that ψ is given by (19)  φ * n,k (y)e imbk E n,k + g 2 dk dy.
Expanding φ n,k in the Wannier basis as in (17)  This is the system of (30), with the coefficients (31)- (34).