Nonlinear Schrodinger equations with a multiple-well potential and a Stark-type perturbation

A Bose-Einstein condensate (BEC) confined in a one-dimensional lattice under the effect of an external homogeneous field is described by the Gross-Pitaevskii equation. Here we prove that such an equation can be reduced, in the semiclassical limit and in the case of a lattice with a finite number of wells, to a finite-dimensional discrete nonlinear Schrodinger equation. Then, by means of numerical experiments we show that the BEC's center of mass exhibits an oscillating behavior with modulated amplitude; in particular, we show that the oscillating period actually depends on the shape of the initial wavefunction of the condensate as well as on the strength of the nonlinear term. This fact opens a question concerning the validity of a method proposed for the determination of the gravitational constant by means of the measurement of the oscillating period.


Introduction
Laser-cooled atoms have drawn a lot of attention as for potential applications to interferometry and high-precision measurements, from the determination of gravitational constants to geophysical applications [13,16,17,22], see also [10,29] for a recent review. The idea of using ultracold atoms moving in an accelerated optical lattice [4,5,21,23,27] has opened the field to multiple applications. In particular, by means of the method proposed by Cladé et al [9], a value for the constant g has been measured using ultracold strontium atoms confined in a vertical optical lattice [12]; such a result has been improved by using a larger number of atoms and reducing the initial temperature of the sample [20]. Determination of g has been obtained by measuring the period T of the Bloch oscillations of the atoms in the vertical optical lattice; recalling that where m is the mass of the Strontium atom, is the Planck constant and b is the lattice period, then a precise value of the constant g has been obtained by means of the experimental measurements of the oscillating period. Since Bloch oscillations with period (1) have been predicted by the Bloch Theorem [8] only for a one-body particle in a periodic field and under the effect of a Stark potential then it has been chosen, in the experiments above, a particular Strontium's isotope 88 Sr; in fact, the scattering length a s of atoms 88 Sr is very small and thus it has been assumed by [12,20] that the effects of the atomic binary interactions are negligible. The obtained value for the constant g was consistent with the one obtained by classical Date: July 16, 2015. gravimeters; but it was affected by a relative uncertainty of order 6 × 10 −6 because of a larger scattering in repeated measurements, mainly due to the initial position instability of the trap. Such a technique is also proposed to measure surface forces [28], too. The critical point of this experimental procedure concerns the validity of the Bloch Theorem and the estimate of the effect of the atomic binary interactions on the oscillating period of the BEC. In order to discuss this point here we are inspired by a realistic model of a one-dimensional cloud of cold atoms in a periodical optical lattice under the effect of the gravitational force. The periodic potential has the shape where b = 1 2 λ L is the period, and λ L = 2π k L . The one-dimensional BEC is governed by the one-dimensional time-dependent Gross-Pitaevskii equation with a periodic potential and a Stark potential where the wavefunction ψ(·, t) ∈ L 2 (R, dx) is normalized to one: ψ(·, t) L 2 = ψ 0 (·) L 2 = 1 , and where is the Bloch operator with periodic potential V per (x). By γ we denote the effective one-dimensional nonlinearity strength.
It is a well known fact (see §6.1 by [8]) that when the wavefunction ψ is prepared on the first band of the Bloch operator and if the nonlinear term is absent, i.e. γ = 0, then the dominant term of the wavefunction ψ exhibits a periodic behavior with Bloch period T within an interval with amplitude B1 |f | , where B 1 is the width of the first band and where f ∈ R is the strength of the external homogeneous field (in the case of f = mg then f takes only positive values, obviously). Therefore, for times of the order of the Bloch period T we may assume that the motion of the BEC occurs in a finite interval. Hence, we can restrict ourselves to the analysis of equation (3) in a suitable finite interval and then we may assume to consider a multiple-well potential V N (x) (with a fixed number N of wells) and that the Stark potential x is replaced by a Stark-type potential W N (x) due to an homogeneous external field which acts only in a bounded region containing the N wells (see Fig.  1). That is, instead of (3) we consider, as a model for a BEC in an optical lattice under an external homogeneous field, the time-dependent non-linear Schrödinger equation (NLS) where > 0 plays the role of the semiclassical parameter (we prefer to denote here the small semiclassical parameter by instead of the usual notation because in a subsequent section we'll discuss a real physical model where will assume its fixed physical value; with such a notation it turns out that the Bloch period is given by T = 2π |f |b ). We assume that the N wells have all the same shape and we denote by b > 0 the distance between the adjacent absolute minima points. The study of the dynamics of the wavefunction ψ, solution of (4), is then achieved by means of a discrete nonlinear Schrödinger equation (DNLS). The idea is basically simple and it consists in assuming that the wavefunction ψ may be written as a superposition of vectors u (x) localized on the −th cell of the lattice; that is Such an approach has been successfully used in the cases of semiclassical NLS with multiple-well potentials [24] or with periodic potentials (see [14,18,19]), without the external field with potential W N . Eventually, u (x) may coincide with the Wannier function u W (x) associated to the first band of the Bloch operator H B or with the semiclassical single well ground state eigenfunction u sc (x). By means of such an approach the unknown functions c (t) turn out to be the solutions of a system of time-dependent equations which dominant terms are given by (here we where λ D is the ground state of a single cell potential and where β is the hopping matrix element between neighboring sites. In fact, the parameter β is expected to be such that 4β is equal to the amplitude B 1 of the first band [25]. In (5) we'll fix c 0 ≡ c N +1 ≡ 0. Equation (5) represents a discrete nonlinear Schrödinger equation (DNLS). Our approach is both semiclassical and perturbative. It is semiclassical in the sense that it holds true in the semiclassical regime of small enough; and it is perturbative in the sense that the external field f and the nonlinearity power strength γ must be small when goes to zero (see Hyp. 3 for details). Under these conditions we prove the validity of the N -mode approximation (5) with a rigorous estimate of the remainder term for times of the order of the Bloch period. Then, we numerically solve the N -mode approximation (5), and we compute the oscillating period taking into account the nonlinear interaction. In fact, the behavior of the wavefunction is not simply periodic in time; it turns out that the center of mass x t = ψ, xψ shows an oscillating motion with modulated amplitude. The oscillating period turns out to be depending on the nonlinearity parameter strength γ and we see that it also depends on the distribution of the initial wavefunction ψ 0 . In particular, when ψ 0 is a symmetric wavefunction then the oscillating period is almost constant for small γ and it practically coincides with the Bloch period T ; on the other hand when ψ 0 is an asymmetrical function the oscillating period actually depends on γ. This fact is in contradiction with the Bloch Theorem (which holds true when γ = 0), which implies that the Bloch period T does not depend on the shape of the initial wavefunction, and it may explain the relatively large uncertainty observed by [20] in their experiments, as discussed in the Conclusions.
The paper is organized as follows. In Section 2 we derive the DNLS (5) from the NLS (4) in the semiclassical limit → 0 for times of the order of the Bloch period T with a rigorous estimate of the remainder term. In particular: in §2.1 we introduce the assumptions and we recall some preparatory results; in §2.2 we derive the DNSL by making use of some ideas previously given by [24] and adapted to the case of multiple-well potential with an external Stark-type perturbation. In Section 3 we consider a realistic experiment and we compute the wavefunction dynamics by making use of the DNLS. In particular: in §3.1 we discuss the validity of the Nmode approximation for different values of the parameters; in §3.2 we numerically compute the wavefunction for times of the order of the Bloch period. In Appendix we write the Wannier functions in terms of the Mathieu functions.
Hereafter, by C we denote a generic positive constant independent of . Let N ∈ N, then by N N := {1, 2, . . . , N } we denote the set of first N positive integer numbers.
By · L p we denote the norm of the Banach space L p (R), by ·, · we denote the scalar product of the Hilbert space L 2 (R).

Derivation of the DNLS (5)
2.1. Assumptions and preliminary results. We consider the time-dependent non-linear Schrödinger equation (4) where V N is a multiple-well potential and W N (x) is a bounded Stark-type potential. In particular we assume that ) smooth function with compact support with a non degenerate minimum value at x = 0: v(x) > v min = v(0), ∀x ∈ R , x = 0.
The multiple-well potential is defined as Hence, by construction the potential V N (x) has exactly N wells with not degenerate minima at x = x , ∈ N N . Remark 1. We assume that v(x) is an even function just for argument's sake. As discussed in Remark 6 this assumption may be removed. Furthermore, we assume that v is a smooth function as usual; in fact, a lessere regularity (e.g. C 2 ) would be enough.
That is the Stark-type potential W N is linear in the region containing the wells and it is a constant function outside this region (see Fig. 1). In the "limit" where N goes to infinity the potential V N becomes a periodic potential with period b and the external potential W N becomes the Stark potential x.

Remark 2.
We restrict ourselves to a multiple-well potential V N with a finite number of wells only for sake of simplicity; one could consider the case of a periodic potential by making use of the tools developed by [14]. On the other side, the assumption on W N is not merely for the sake of simplicity; actually, the Stark-type potential W N is a bounded operator while the Stark potential x is not a bounded operator and this fact is a source of several technical problems. In fact, in real experiments the BEC are trapped in a finite spatial region.
Hypothesis 3. We assume to be in the semiclassical limit, that is we look for the solution of (4) in the limit of that goes to zero. We assume also that the other two parameters γ and f are small for small. That is we assume that there exists > 0 such that for some s > 2, C > 0 and ρ ∈ (0, S 0 ) independent of ; furthermore, we assume also that for some positive constant C and for any ∈ (0, ).
The self-adjoint extension of the linear Schrödinger operator formally defined on L 2 (R) as has an almost degenerate ground state with dimension N . More precisely, let λ , ∈ N N , be the lowest eigenvalues of H N with associated normalized eigenvectors v . In particular we have that (see Lemma 2 [25]) is the Agmon distance between two wells and λ D is the ground state of the single well operator − 2 ∂ 2 xx + v, where the single well potential v has been introduced by Hyp. 1. The numerical pre-factor β is the hopping matrix element between neighboring wells, and it is such that 4β is asymptotic to the amplitude of the first band of the periodic Bloch operator 1 are, respectively, the bottom and the top of the first band. Such a numerical pre-factor is going to be exponentially small, i.e.
β =Õ(e −S0/ ) as → 0 + . Remark 3. Hyp. 3 means that, from a practical point of view, the parameter f cannot be arbitrarily small, but it has a lower bound of order β. On the other hand, the parameter γ may be arbitrarily small.
The associated normalized eigenvectors are given by [25] where α j, = α ,j = 2 N + 1 sin j π N + 1 and where u sc j (x) is the semiclassical single well ground state eigenfunction localized on the j-th cell; by construction and since v(x) is an even function then Now, let Π be the projection operator associated with the N eigenvalues λ , i.e.
Remark 4. Let σ(H N ) be the spectrum of H N ; then it is a well known semiclassical result that Remark 5. By [14] it has been proved that there exists a suitable orthonormal base u , ∈ N N , of the space F . The functions u are practically localized on the -th well. More precisely, they are such that The matrix with elements u , H N u j can be written as where T is the tridiagonal Toeplix matrix such that for some α > 0.
We finally assume that the initial state is prepared on the first N "ground states". That is It is well known that under the assumptions above the NLS (4) is locally well posed, and the conservation of the norm and of the energy [6,7] Furthermore the following a priori estimate follows, too. Lemma 1. There exists a positive constant C > 0 such that Proof. Indeed, from Theorem 2 by [24] and its remarks it follows that for some C > 0 and small enough, where Hence, the global well-posedness of the NLS follows [6,7].

N-mode approximation.
Let ψ be the normalized solution of the NLS equation written in the formula for some complex-valued functions c (t). By substituting (8) into the NLS (4) then it takes the formula We are going now to get an a priori estimate of the remainder term ψ c . First of all we rescale the time t → τ = β t and we redefine the wavefunction up to a gauge factor ψ( Hence, (9) becomes (where denotes the derivative with respect to τ ) Theorem 1. Let Hyp.1-4 be satisfied; then it follows that the remainder ψ c can be estimated for times of order of the Bloch period. That is for any fixed M ∈ N it follows that Proof. From the first equation of (10) and recalling that then a priori estimate follows because u L 2 = 1 and ψ L 2 = 1, and from Remark 5 iv. and Lemma 1.
Concerning ψ c it satisfies to the following integral equation where we set and where A and B are defined as By means of standard arguments [24] and making use of the fact that the operator W N is bounded then it follows that Lemma 2. Let Then the functions A and B are such that Proof. Indeed, Similarly, the estimate of the function B follows recalling that ψ c L ∞ ≤ C −1/4 from Lemma 1. Finally, the estimate concerning ∂A ∂τ immediately follows from (11). Hence, the estimates of the integrals I and II follow; in particular, integral II can be simply estimated as On the other hand, before to get the estimate of integral I we perform an integration by parts in order to gain a pre-factor β: From this fact and recalling that (Remark 4) Therefore, we have that From the Gronwall's Lemma it follows that In particular we observe that proving the Theorem since Γ ≤ C|f | from Hyp. 3 and since τ B = 2π b β |f | . We are going now to estimate the solutions c of the first equation of (10) which can be written as where the term u , (H N − λ D )ψ 1 can be represented by property iv. of Remark 5. Concerning the term u , B the following estimate uniformly holds with respect to the index We have that Lemma 3. The following estimates uniformly hold with respect to the indexes , j , m and k: (7) and Remark 4. Therefore where u 0 is normalized and I0 |u 0 (x)| 2 xdx =Õ e −S0/ because u 0 (x) = u 0 (−x)+ O e −S0/ . From this fact and since u L 2 (R\I ) =Õ(e −S0/ ) (see Lemma 4 iii. and Lemma 5 by [14]) then the asymptotic behavior i. follows. The other two asymptotic behaviors ii. and iii. similarly follow from property ii. by Remark 5; indeed and, where we assume that r = |j − |, proving so the estimates ii. and iii..

From this Lemma and from the previous computation it follows that the first equation of (10) becomes a DNLS of the form
where and where the remainder termsÕ e −S0/ are uniform with respect to the index .

Remark 6.
In fact, if v(x) is not an even function then by means of standard semiclassical arguments it follows that property Lemma 3 i. becomes for some constant c independent of the index . In such a case we must add the term f cc to the right hand side of the DNLS above and, by means of a gauge choice c → c e −i f c β τ , we can remove this term obtaining again equation (12). Now, we are able to prove that Proof. The proof is a simply consequence of equation (12) and from the fact that τ B = 2π b β f and Hyp.3.

Numerical analysis of a real model
We consider the experiment where a cloud of ultracold Strontium atoms 88 Sr are trapped in a one-dimensional optical lattice with potential (2). Realistic data for the experiment are [20]: -Lattice period: b = λ L /2 = 266 nm, λ L = 532 nm; -Lattice potential depth: where d ⊥ is the oscillator length of the transverse confinement; here a s denotes the scattering length of the Strontium 88 isotope: a s = −a 0 ÷13a 0 , where a 0 is the Bohr radius; N is the number of atoms of the condensate; in typical experiments d ⊥ ≈ 180 · 10 −6 m and N = 10 5 ÷ 10 6 ; -Acceleration constant g = 9.807 m/s 2 . The confined BEC is governed by Eq. (4) and here we make use of the N -mode approximation (13), that is the wavefunction ψ has the form ψ ∼ c u where c are the solutions of (13) and where u are functions localized on the -th lattice site. In order to justify the validity of such an approximation we'll check if the model is in the semiclassical regime, that is if the first band is almost flat and if semiclassical approximation u sc agrees or not with the Wannier function u W . Such a qualitative criterion has been also adopted by other authors [2,3,11] and we'll see that our results agree with the results contained in these papers. In particular, in [11] has been computed the hopping matrix elements u , H N u j too, where it has been numerically verified that these coefficients are negligible when |j − | > 1 for Λ 0 ≥ 10; thus, for such values of Λ 0 it is admitted that the N -mode approximation, consisting to describe (4) in terms of a nearest-neighbor model (13), works.
3.1. Validity of the semiclassical approximation. The semiclassical approximation u sc 0 (x) of the wavefunction has dominant behavior in the semiclassical limit, where µ = We may remark that the effective semiclassical parameter in adimensional units is given by and then the semiclassical approximation may be written as We'll see that for Λ 0 "large enough" (i.e. Λ 0 ≥ 10) then the first band is almost flat and the semiclassical function u sc 0 well approximates the Wannier function u W 0 , as we expect to observe in the semiclassical limit Λ 0 → ∞ (see, e.g., [30]).
Remark 7. By the scaling Equation (15) is equivalent, up to a change of scale of the time, to the equation where we set and where = Λ −1/2 0 plays the role of the semiclassical parameter.
We compute now the band functions and the Wannier functions for different values of Λ 0 . The semiclassical wavefunction u sc 0 is computed by (14), while the Wannier function u W 0 (x) may be computed by means of the Mathieu functions (see Appendix).
3.1.1. Model Λ 0 = 3. The first bands of the Bloch operator Hence, the values of the width of the first two bands are given by Furthermore it follows that the first gap has amplitude g 1 = E b 2 − E t 1 = 0.77 · E R of the order of the first band amplitude, while the width of the other gaps are very small (see Fig. 2, left hand side panel). If we compare the first Wannier function u W 0 (x) and the semiclassical approximation u sc 0 (x) it turns out that (see also .32 · E R and E t 1 = 4.58 · E R ; n=2) E b 2 = 7.02 · E R and E t 2 = 8.87 · E R ; n=3) E b 3 = 9.54 · E R and E t 3 = 14.07 · E R . Hence, the values of the width of the first two bands are given by Furthermore it also follows that the first gap has amplitude g 1 = E b 2 −E t 1 = 2.44·E R is much larger than the amplitude of the first band and that the width of the other gaps are very small (see Fig. 3, left hand side panel). If we compare the first Wannier function u W 0 (x) and the semiclassical approximation u sc 0 (x) it turns out that (see also Fig. 3, right hand side panel) 055 . Hence, we may conclude that for Λ 0 = 10 the N -mode approximation properly works.

3.2.
Numerical analysis of the model for Λ 0 = 10. We have seen that for Λ 0 ≥ 10 the N -mode approximation is justified. For Λ 0 = 10 we have that Equation (  The Bloch period is given by Hence, the parameters f , γ and β are in a suitable range as discussed in Remark 3. Furthermore, the motion of the Bloch oscillator occurs in an interval with width Hence, the N -mode approximation with N = 40 properly works. We consider three different situations. In the first one we assume that the state is initially prepared on a single lattice site, that is ψ 0 is a Wannier type function. In the other two cases we assume that the initial wavefunction ψ 0 is a symmetric or asymmetrical wavefunction initially prepared on different lattice sites. 3.2.1. ψ 0 is initially prepared on a single lattice cell. We consider a numerical experiment where ψ 0 (x) = u 0 (x), that is c (0) = 0, for = N/2, and c N/2 (0) = 1 (where N = 40). In fact, in such a case we observe a breathing motion for the wavefunction; that is, the wavefunction, initially prepared in a Wannier state localized on a single site of the optical lattice, symmetrically spreads in space and it periodically returns to its initial shape (Fig. 4, top panel, obtained for η = 0.2). Then the expected value of the center of mass is practical constant x t ≈ 0 up to small fluctuations. Figure 4. In the top panel we plot the absolute value of the wavefunction ψ(x, t) initially prepared on a single Wannier state for η = 0.2, it turns out that it symmetrically spreads in space and periodically returns to its initial shape without motion of the center of mass. In the bottom panel we plot the absolute value of the wavefunction initially prepared on several lattice sites for η = 0.2; it turns out that the center of mass oscillates with no marked changes of the shape of the wavefunction. Here T denotes the Bloch period and b is the distance between two adjacent wells. Dark regions mean that |ψ(x, t)| is practically zero there, white regions mean that |ψ(x, t)| has its maximum value there.  Table 1. Initial values of the coefficients c := c (0) of the wavefunction. The initial wavefunction ψ0 has a symmetric shape and its width is of order of several lattice periods.

ψ 0 is a symmetric wavefunction initially prepared on different lattice cells.
We consider a numerical experiment where N = 40 and ψ 0 (x) = 40 =0 c u (x), where c have a symmetric Gaussian-type distribution around = N/2. That is the initial value of the coefficients c (t) are given in Table 1, the initial wavefunction ψ 0 is plotted in Fig. 5, left hand side panel. In such a case the center of mass x t oscillates in space and the wavefunction moves with no marked changes in shape (see Fig. 4, bottom panel). In particular, the function x t exhibits, for η = 0, an oscillating motion where the wavefunction amplitude is modulated (see Fig. 6) and where the oscillating (pseudo-)period (that is the time interval between two consecutive minima or maxima points) depends on η. In Fig. 7 we plot the mean value of the oscillating period of the motion of the center of mass after 14 oscillations for η in the range [−0.1, +0.2]; it turns out that the relative uncertainty with respect to the Bloch period is of order 2.4 · 10 −5 .
3.2.3. ψ 0 is an asymmetrical wavefunction initially prepared on different lattice cells. We consider a numerical experiment where N = 40 and ψ 0 (x) = 40 =0 c u (x), where c have an asymmetrical Gaussian-type distribution. That is the initial value of the coefficients c (t) are given in Table 2, the initial wavefunction is plotted in Fig. 5, right hand side panel. As in the symmetric case the center of mass x t oscillates in space and the wavefunction moves with no marked changes in shape. Even in such a case the function x t exhibits, for η = 0, an oscillating motion where the wavefunction amplitude is modulated. In contrast with the symmetric case the oscillating (pseudo-)period (that is the time interval between two consecutive minima or maxima points) actually depends on η; in Fig. 7 we plot the mean value of the oscillating period of the center of mass after 14 oscillations for η in the range [−0.1, +0.2] and it is not almost constant like in the previous case, in particular it turns out that the relative uncertainty with respect to the Bloch period is of order 4.6 · 10 −4 , which is 20 times the relative uncertainty observed in the symmetrical case.   Table 2. Initial values of the coefficients c := c (0) of the wavefunction. The initial wavefunction ψ0 has an asymmetrical shape and its width is of order of several lattice periods.

Conclusion
In this paper we have proved that in the semiclassical limit the N -mode approximation (13), corresponding to a discrete nonlinear Schrödinger equation with a finite number of modes, gives the solution of the Gross-Pitaevskii equation (4) for a BEC in a multiple-well lattice in a Stark-type external field. Furthermore, we have numerically solved the N -mode approximation considering a real model, where for some values of the physical parameters the validity of the N -mode approximation (13) seems to be justified. In particular, we have seen that a state initially prepared on several wells have an oscillating behavior with modulated amplitude, the oscillating (pseudo-)period is computed for different values of the nonlinear strength and it turns out that such a period is practically constant when the initial state Figure 6. Here we plot the motion of the center of mass of the wavefunction initially prepared on several lattice sites. The initial wavefunction is a symmetric function. Top panel corresponds to the case of η = 0.1, bottom panel corresponds to the case of η = 0.2. The center of mass rapidly oscillates with modulation of the amplitude. The width of the escillations is in a range lesser or equal to 3b is agreement with (16). is a symmetric one; on the other side, such a period actually depends on the nonlinear strength when the initial state is an asymmetrical one. This observation opens a question about the validity of the method proposed by Cladé et al [9] for the deterimantion of the gravitational constant g by means of the measurement of the oscillating period [12,20], where it has been assumed that the oscillating period coincides with the Bloch period T independently from the shape of the initial wavefunction and of the value of the nonlinearity strength parameter. Here we plot the mean value of the pseudo-period of the oscillating motion of the center of mass after 14 oscillations, as function of the effective nonlinearity parameter η. Broken line corresponds to the case of a symmetric wavefunction prepared on several lattice sites; it turns out that in such a case the oscillating period is almost constant. Full line corresponds to the case of an asymmetrical wavefunction prepared on several lattice sites; it turns out that it actually depends on η. Here T denotes the Bloch period, while t denotes the oscillating period.

Appendix A. Band Functions and Wannier functions
For a generic one-dimensional Bloch operator H B the spectrum is given by a sequence of infinitely many closed intervals named bands. These intervals are the image of functions named band functions. The band functions of H B are denoted by E n (k), where the quasimomentum k runs in the Brillouin zone − π b , + π b . The spectrum of the Bloch operator H B is given by the bands σ( if n is odd . In the case of potential (2) the band functions may be explicitly computed. In particular let us look for the Bloch functions of the equation Hence, the band functions E n (k) associated to the spectral problem (17) are the solutions of the equation µ(E) = cos(kb) where Let λ = e ikb , then the equation µ(E) = cos(kb) can be written as µ(E) = 1 2 (λ + λ −1 ). We observe that for k ∈ 0, π b then sin(kb) = 1 − µ 2 (E). The Bloch function is given by [15] ψ(x, E) = χ(x, E) and We recall that the Bloch function ψ n (x, k) = ψ(x, E n (k)), where E n is the band function associated toH B , is normalized to one: and furthermore it is such that ψ n (x, −k) = ψ n (x, k) .
Finally, the Wannier function on the zero-th cell associated to the n-th band is given by N (E n (k)) dk since the Mathieu functions are real valued when their arguments are real numbers. In particular, of the Isaac Newton Institute for Mathematical Sciences where part of this paper was written.