Wannier functions and discrete NLS equations for nematicons

Abstract: We derive nonlocal discrete nonlinear Schrödinger (DNLS) equations for laser beam propagation in optical waveguide arrays that use a nematic liquid crystal substrate. We start with an NLS-elliptic model for the problem and propose a simplified version that incorporates periodicity in one of the directions transverse to the propagation of the beam. We use Wannier basis functions for an associated Schrödinger operator with periodic potential to derive discrete equations for Wannier modes and propose some possible simplified systems for interactions of modes within the first energy band of the periodic Schrödinger operator. In particular, we present the simplest generalization of a model proposed by Fratalocchi and Assanto by including a linear nonlocal term, and see evidence for parameter regimes where nonlinearity is more pronounced.


Introduction
We study the derivation of nonlocal discrete nonlinear Schrödinger (DNLS) models for a system describing the propagation of a laser light in a nematic liquid crystal substrate with a periodic structure in the direction normal to the optical axis. This is the setup of laser light propagation in an effectively planar waveguide array in a nonlocal medium. Our work is motivated by the experimental work of [1][2][3] and some of the discrete NLS systems proposed to model the experiments. We are particularly interested in a DNLS model proposed by Fratalocchi and Assanto [4] that has a Hartree-type cubic nonlinearity. Recent theoretical and numerical studies have shown interesting departures from the dynamics of the cubic DNLS, specifically new non-monotonic amplitude profiles of breather solutions, additional internal modes in the linearization around breathers [5,6], and possibly enhanced mobility of traveling localized solutions [7]. The goal of our study is to derive and generalize the Fratalocchi-Assanto model.
The derivation of DNLS equations can be made systematic by the use of Wannier functions [8], and this formalism has been also used to show some rigorous justifications of the DNLS equations in a special limit referred to as the tight binding approximation [9,10].
The present work adapts this formalism to a nonlocal problem. A first step is to simplify the original equations describing the interaction of the laser field with the liquid crystal orientation field, deriving an NLS equation coupled to a linear elliptic equation with periodic coefficients, see Section 2. Wannier functions are constructed and computed using a periodic Schrödinger operator that appears naturally in the elliptic equation. We see that we can solve the elliptic equation using the Wannier basis and formally reduce the system to an evolution equation for the Wannier coefficients. The coefficients describing linear and nonlinear interactions between the Wannier modes involve integrals of Wannier functions and related functions and must be evaluated numerically. We examine a subset of these interactions with a first goal of obtaining the simplest extension of the Fratalocchi-Assanto model. We arrive at a system combining a nonlinear nonlocality of Hartree type similar to the one in the Fratalocchi-Assanto model with a linear nonlocality that is studied mathematically by several authors in systems with power law nonlinearity [11][12][13]. Some related nonlocal equations with periodic coefficients were studied in [14,15].
The computed coefficients suggest that nonlocal effects are more pronounced when the periodicity in the transverse direction has a smaller amplitude, and we also see parameter regimes where the nonlinear part dominates. For larger amplitude periodicity the standard cubic NLS seems to be an adequate model. The model we examine in more detail here is restricted to interactions among Wannier modes of the first energy band of the associated Schrödinger operator. Preliminary calculations suggest that other interactions within the first band and with other bands may also be important, see also [1] for experimental evidence. Extensions will be pursued in further work.
The paper is organized as follows. In Section 2 we derive a simplification of the physical model, and identify a periodic Schrödinger operator used to define the Wannier basis. In Section 3 we use the Wannier basis to arrive at a general evolution equation for the Wannier coefficients of the laser field amplitudes. In Section 4 we present numerical results on the coefficients for a small subset of the possible indices, considering only the lowest energy band. In Appendix A we present additional material on the numerical computation of the Wannier functions.

Coupled NLS-elliptic system for optical solitons in liquid crystals
Optical beams in a nematic liquid crystal are described by the system with γ, ν positive constants, see [16,17], and ∆ = ∂ 2 x + ∂ 2 y . The variable z is analogous to time and we are interested in solutions to (2.1), (2.2) with u, θ given at z = 0. Physically u represents the amplitude of the electric field of a laser beam that propagates along the z−direction. The field u is assumed to have only one component, along the (vertical) x−axis. θ is the macroscopic field of angles describing the liquid crystal molecular orientation at each point, i.e., a macroscopic local average orientation of the liquid crystal molecules, assumed to be parallel to the x, z−plane. θ 0 and E 0 are assumed to be known functions. θ 0 represents an orientation field produced by an additional electric field E 0 applied to the sample. The field E 0 is assumed to have only an x−component.
Equations (2.1), (2.2) follow from a Maxwell-Oseen-Frank system that couples the electric fields E 0 and u and the director angle θ, see [18], Appendix. Intuitively, (2.2) describes the tendency of the molecular and electric field orientations to align, as well as an elastic resistance of the material to localized changes in orientation, while (2.1) describes the focusing effect of aligned molecular and electric field orientations.
To enhance this basic effect, it has been understood that the system must be "prepared" by first producing a nontrivial "pretilt" angle field θ 0 in the sample. We are here considering an experimental setup where θ 0 is produced by a known external field E 0 , e.g., placing the sample between two capacitor plates. Other experimental ways to produce θ 0 are described in [19,20]. θ 0 and E 0 are thus related, see below. The main point, however, is that these two functions can be used to introduce variable coefficients in the system, and we will consider the case where θ 0 , E 0 are b−periodic in the horizontal direction y, and independent of the z coordinate.
To describe the relation between θ 0 and E 0 in the absence of the laser beam, we use (2.2) with θ = θ 0 , u = 0, see [18] for the boundary conditions. In the presence of a laser beam, the laser field interaction changes the angle field to θ = θ 0 + ψ, i.e., ψ is the change from the pretilt angle θ 0 observed in the absence of the beam. Letting θ = θ 0 + ψ and substituting (2.3) into (2.2) we obtain the system see [16]. The second equation has been examined in more detail in [18] for the case of E 0 , θ 0 constants, with θ 0 ∈ ( π 4 , π 2 ). This assumption is important for solving (2.5) uniquely, see [18]. Here we are looking for an explicit expression of ψ in terms of u and we simplify the system first assuming ψ small. Then system (2.4), (2.5) is, up to terms of O(ψ 2 ), see also [18]. Further analysis of the equations assumes that E 0 , θ 0 are known. To simplify the problem we will choose E 0 , θ 0 that are reasonably consistent with (2.3), i.e., avoiding solving (2.3) exactly. One way to do this is to assume E 0 , θ 0 of the form with 0 < ≤ δ small parameters. The applied electric field E 0 is thus assumed to consist of two parts, a contribution √ E 0,V that depends only on the vertical direction x, and a smaller contribution δE 0,H that will be assumed b−periodic in the horizontal direction y. E 0 models the electric field of a set of parallel metallic stripes that run along the z−axis and are placed at the two planes bounding the experimental domain from above and below, at x = ±L respectively, see [1]. The lower and upper stripes have opposite charges, producing an electric field that is b−periodic in the y (transverse) direction.
Matching orders in (2.3) up to O( 3/2 δ) leads to The first equation is similar to the one arising in the case of constant applied field, i.e., compare to (2.3) with E 0 constant, and is solved using boundary conditions for θ 0 at the vertical boundaries x = ±L, see e.g., [17]. In the second equation E 0,H ( √ x, y) is assumed b−periodic in y, and the solution θ 0,H ( √ x, y) can also be chosen b−periodic in y.
To use (2.8)-(2.11) in (2.6), (2.7), we note that since E 0 , θ 0 vary slowly in the vertical direction x, their value in the center of the sample can be approximated by constants. We therefore substitute γ |sin 2θ 0 | + δ 1/2 cos 2θ 0 θ 0 (y) |u| 2 , (2.12) up to terms of O(δ 2 ) on the left hand side, and terms of O(δ 2 ) on the right hand side. Thus the terms omitted on each side can be considered of O(δ 2 ). Note that E 0 (y), θ 0 (y) are b−periodic, and are related by the simplification of (2.11). Thus θ 0 can be determined from E 0 , moreover if E 0 is b−periodic, then θ 0 can be also chosen b−periodic.
We further simplify the analysis by considering u, ψ independent of the vertical variable x. Letting δ = , and considering terms of up to O( 3/2 ) in (2.12), system (2.6), (2.7) then becomes The first equation follows from (2.6) using the expansion (2.9) for θ 0 and the simplifications above, in particular β = | sin 2θ 0 | + O( 3/2 ), where the omitted O( 3/2 ) terms are b−periodic in y. The second equation is (2.12) up to terms of O( 2 ), so that g 2 is a positive constant, and α, V b−periodic functions, defined implicitly by comparing the second equation to (2.12).
Since we are primarily interested in the inhomogeneity of the second equation and the inversion of the operator −∂ 2 y + V(y) + g 2 , we will only consider the lowest order, constant parts of α, β of the respective α, β above, and consider the simplified system where assuming θ 0 ∈ ( π 4 , π 2 ), are positive constants, and We can assume that V has a negative minimum −V m that satisfies V m < g 2 , with g 2 as in (2.18). We can then suitably redefine g 2 , and V of (2.19) so that both are positive. These sign assumptions are valid also for the variable coefficients α(y), β(y) of the intermediate system above. The goal of the paper is to write further explicit simplifications of system (2.16), (2.17) using the Wannier basis defined by the Schrödinger operator −∂ 2 y + V(y), with V as in (2.19). An alternative approach to simplifying the system is to write the inverse of the operator −∂ 2 and keep the lowest order terms in the expansion Assuming |V| bounded, both A −1 , and V are bounded operators in L 2 (R, C), with norms g −2 (by Plancherel's theorem), and V M = max y∈R |V(y)| respectively. The g −2 V M < 1 implies that A + V has a bounded inverse given by the convergent series expansion (2.21). We can also use the explicit expression where f is a function and K * f denotes the convolution Then (2.16) becomes The lowest order approximation that includes the effect of the periodic potential V is The fact that the variable coefficient equation has additional terms seems inconvenient for computations, and in the next section we describe a different approach to solving (2.17).

The discrete nonlinear equation for nematicons
We consider the system of periodic nematicon equations derived in the previous section with α, β, g 2 positive constants, and V b−periodic and positive. We recall the definition of the Wannier functions associated to the operator −∂ 2 y +V(y), V b-periodic, see [10,21] for the theory. Bounded solutions φ n,k (Bloch functions) and eigenvalues E n,k ∈ R of the periodic Schrödinger equation satisfy for all n ∈ N, k, y ∈ R. By the above 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 (3.2) with boundary conditions φ n,k (y + b) = e ikb φ n,k (y), implied by (3.3). 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 (3.5) are referred to as Wannier functions [21,22], They form an orthonormal set of L 2 (R; C), also by (3.5), (3.3), Thus, fixing n, the function w m n is a translation of the function w 0 n by mb. The definition of w 0 n and the regularity of the φ n,k in k also leads to strong localization of the w 0 n in y, see [10,23]. We use expansions ψ and u in Wannier functions w m n as We first consider the second nematicon equation (3.1b). We multiply (3.1b) by φ * n ,k (y), integrate over y ∈ R, and use the Schrödinger equation We then multiply (3.11) 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 (3.8) with To expand the first nematicon equation (3.1a) in coefficients of Wannier functions, we substitute the series expansions (3.8), (3.9) into (3.1a), multiply (3.1a) by w m * n (y), and integrate over y ∈ R. We obtain du n ,m dz = where by (3.13), (3.15) The above infinite system of ordinary differential equations for the Wannier coefficients is a spatial discretization of the original system of partial differential equations, and a first step to simplifying truncations leading to discrete NLS equations and systems [8]. In the present work, the main goal is a systematic derivation of nonlocal versions of DNLS equations and systems, motivated by the Fratalocchi-Assanto equation [4]. Important features of the Wannier functions w m n include their spatial localization ( [10], see also Appendix), and the translation property (3.7). Truncations involving m ∈ Z, n ∈ {1, . . . n max }, are expected to preserve the translation symmetry of periodic systems, and we increase the spatial resolution by considering larger n max .

Comparison with the discrete model of Fratalocchi and Assanto
We now consider simplifications of (3.16), (3.17) leading to a DNLS model that is comparable to the equation   We now examine the symmetries of the linear and nonlinear coefficients of (3.21).
Considering the linear coefficients, by (3.15) we have (3.25) Letting m 2 = m 0 − m 1 and m 3 = m − m we then have The Hamiltonian structure of (3.21) is not surprising given that the continuous system (2.16), (2.17) we expanded in Wannier functions is also Hamiltonian. This follows from the fact that the operator (A + V) −1 in (2.24) is symmetric, see e.g., [24]. A more systematic derivation of the Hamiltonian structure of the discrete system, considering first the Hamiltonian structure of (3.16), will be presented elsewhere. Note that system (2.4), (2.5) is also Hamiltonian in the sense described in [25], see also [18] for a similar system.

Numerical evaluation of linear and nonlinear coefficients
To evaluate the linear and nonlinear coefficients of (3.21) we have used the Wannier functions of the periodic square-well potential of (A.2). Appendix A contains further details. We are interested in the dependence of the coefficients on the amplitude V 0 of the potential and consider values V 0 ∈ [0.1, 10]. In all cases, the potential has period b = 2π.
The linear coefficients D m,m 1,1 = D 1,1 (m − m ) are plotted in Figure 1 for different potential amplitudes V 0 . Their values for nearest neighbors are shown in Table 1. The nonlinear coefficients V 1 (m 1 , m ) of (3.26) are shown in Figure 2. They depend on the potential amplitude V 0 and g 2 . Their values for nearest neighbors are shown in Table 2.

Figures 1 and 2 show clearly that both linear and nonlinear coefficients, D m,m
1,1 and V 1 (m 1 , m ) respectively, decay more slowly as the amplitude V 0 is decreased. The nonlinear coefficients also decrease with g 2 , Figure 2, as expected from their definition, e.g., (3.26).
Also, by Figures 1 and 2 case V 0 = 0.1, the signs of the linear and nonlinear terms are compatible to the ones assumed in (3.18), and correspond to the focusing sign combination. The off-diagonal linear coefficients are also rather small (∼ 0.05) compared to the nonlinear coefficients (∼ 1), we are thus near the "anti-continuous" limit considered for the nonlocal problem in [5,6]. A similar observation applies to the larger potential amplitudes shown, V 0 = 1, 10. (The diagonal linear coefficient D m,m 1,1 is independent of m and can be eliminated by a phase change.) In those two cases the cubic power DNLS seems to be an adequate model. Figure 2 shows the decay of the nonlinear interaction is faster for larger V 0 .  super-exponential decay |V 1 | ∼ Ce −γ|∆m | ξ with γ ∼ 0.8 and ξ = 1.17, for V 0 = 1 and small g ∼ 0.1; exponential decay |V 1 | ∼ Ce −γ|∆m | with γ ∈ [1.93, 1.96] according to g, for moderate V 0 = 1; and subexponential decay |V 1 | ∼ Ce −γ|∆m | ξ with γ ∈ [11,14] and ξ ∈ [0.7, 0.9] according to g, for V 0 = 10, see main text for details. Table 2. Nearest neighbor nonlinear coefficients V 1 (∆m ) for g = 0.1, 1 and 10 for potentials V 0 = 0.1, 1, and 10.

Discussion
We derived a model NLS-elliptic system for an optical waveguide array in a nonlocal medium, and used Wannier functions defined by an associated Schrödinger operator with periodic potential to write an infinite system for the interaction of the Wannier modes of the laser beam amplitude. We also discuss a simplified system that considers a subset of the interactions of Wannier modes whithin the first energy band of the periodic Schrödinger operator. We obtain a generalization of a nonlocal DNLS model derived by Fratalocchi-Assanto [4]. We see that nonlocal effects are more pronounced for small amplitude oscillations of the system coefficients across the waveguide array.
The NLS-elliptic model we use is a variant to the one considered by [4], where the authors use alternative perturbative arguments to solve approximately the equations for the pretilt angle θ 0 and for an additional angle ψ produced by the laser beam. In the present work we focus on the equation for the additional angle ψ, and solve it by introducing expansions in Wannier functions. Another difference is in the treatment of the equation for the electric field, where [4] uses a coupled mode Ansatz, with mode couplings of a predetermined. The present paper considers a related but more general idea, expansion in a Wannier basis [8]. The computation of the Wannier functions is involved, but advantages include the fact that the Wannier mode interactions are determined by the continuous equation. Also the general discrete system that describes the interactions among all modes can be simplified by considering a finite (and small) number of Wannier modes per site, and a subset of the mode interactions. Calculations that will be reported elsewhere suggest that these simplifications can be also made at the level of the Hamiltonian, preserving the Hamiltonian structure of the continuous model.
The simplification leading to discrete equations with the structure of the Fratalocchi-Assanto model suggests nonlinear interactions that decay roughly in an exponential way, as in [4]. The construction of [4] considers only nearest-neighbor linear interactions. In cases of slower decay for the nonlinear coupling, i.e., more pronounced nonlocality, our results suggest that we could also consider nonlocal linear coupling terms. Even in these cases however, linear couplings decay faster that nonlinear couplings, justifying the simplification in [4], and also suggesting the relevance of the anti-continuous limit considered in [5,6].
In future work we plan to examine models describing a larger class of interactions within the first band, as well as interactions with the second band. Higher local modes may be especially relevant for higher nonlocality, and are also relevant to experimental observations [1][2][3]. It would be also of interest to construct discrete models with saturation effects [18], starting from (2.4), (2.5). The formalism described so far does not apply directly.

A. Appendix. Computation of Wannier functions and energies
A.1. Two linearly independent solutions ψ 1 (y) and ψ 2 (y) of the Schrödinger equation for the b-periodic square-well potential where V 0 > 0 is the potential amplitude and ϑ its width. Results for some spacial cases are in [10,26]. Solutions of (A.1) have eigenenergies that belong to intervals called energy bands. Ψ's solutions may be expressed in terms of the so-called Bloch wave functions u n,k (y), that have the same period as V, as Ψ n,k (y) = e iky u n,k (y).
k is the parameter that allows to explore values within the energy bands For the general solution of (A.1), we first find two linearly independent solutions ψ 1 (y) and ψ 2 (y) assuming the initial conditions For E n,k > V 0 , solutions are and for E n,k < V 0

A.2. Bloch functions
The Schrödinger equation (A.1) is a second-order linear differential equation. The general solution Ψ(y) = Ψ n,k (y) may be computed as a linear combination of the independent solutions ψ 1 (y) and ψ 2 (y). Such linearity also applies to the derivatives, leading to the system Ψ(y) = Aψ 1 (y) + Bψ 2 (y), (A.11) Ψ (y) = Aψ 1 (y) + Bψ 2 (y). (A.12) By assuming the initial conditions Ψ(y + 0 ) and Ψ (y + 0 ), using y + 0 = 0 + , we eliminate the constants A and B. The resulting system of equations is written in matrix form following the procedure of Bruno and Nacbar [23] Ψ(y) Ψ (y) = ψ 1 (y) ψ 2 (y) ψ 1 (y) ψ 2 (y) where the 2 × 2 matrix in the right-hand side is known as the transfer matrix T from y + 0 to y. Eigenfunctions Ψ(y) and Ψ (y) satisfy the Bloch theorem, i.e., Ψ n,k (y + b) = e ikb Ψ n,k (y). Evaluation of (A.13) at y = y + 0 + b, and use of the Bloch theorem leads to where (M i j ) (i, j = 1, 2) are elements of the primitive matrix M defined by . Ψ n,k of (A.23) depends explicitly on k through cos kb in the real part and sin kb in the imaginary part. Ψ n,k also depends implicitly on k though E n,k = E n (k) by Ψ(y + 0 ), M 12 , M 11 , ψ 1 (y) and ψ 2 (y). If E n (k) has even parity on k then all terms implicitly depending on k are also even by their definition. In such a case, Ψ n,k (y) obeys the inversion symmetry Ψ n,−k (y) = Ψ * n,k (y) and Wannier functions are real, according to (3.5)  We numerically solve (A.18) and (A. 19) to find E n,k as a function of k, for the square-well potential of period b = 2π and amplitude width ϑ = π. Solutions were computed using Newton's method. As potential amplitudes we chose V 0 = 0.1, 1, and 10. Figure 3 shows results of the lowest energy bands E n (k) in the range k ∈ [−1/2, 1/2]. Results show even parity on k.  We used the eigenenergies to evaluate the real and imaginary parts of the Wannier functions from (3.5) using the Bloch functions (A.23). Wannier functions are real, in agreement with the inversion symmetry of Bloch functions [23]. The integral in (3.5) is computed by Simpson's rule using a grid with 10 3 points. Figure 4 shows Wannier functions w 0 n (y) of the first five energy bands (for n = 1, . . . , 5) for potentials V 0 = 0.1, 1, and 10. Plots range from y ∈ [−15b, 15b] in lattice units. Wannier functions show oscillations that decay along the y-coordinate. The number of oscillations increases with the band number n, and decreases with the amplitude of the potential V 0 . Oscillations still persist after several periods from the origin, primarily for large energy bands.