1 Introduction

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

## Abstract.

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 Schrödinger 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.

## 1. 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 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 has been obtained by measuring the period of the Bloch oscillations of the atoms in the vertical optical lattice; recalling that

 T=2πℏmgb, (1)

where is the mass of the Strontium atom, is the Planck constant and is the lattice period, then a precise value of the constant 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 ; in fact, the scattering length of atoms 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 was consistent with the one obtained by classical gravimeters; but it was affected by a relative uncertainty of order 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

 Vper(x)=V0sin2(kLx) (2)

where is the period, and .  The one-dimensional BEC is governed by the one-dimensional time-dependent Gross-Pitaevskii equation with a periodic potential and a Stark potential

 iℏ∂tψ=HBψ+fxψ+γ|ψ|2ψ, f=mg, (3)

where the wavefunction is normalized to one:

 ∥ψ(⋅,t)∥L2=∥ψ0(⋅)∥L2=1,

and where

 HB=−ℏ22m∂2xx+Vper(x)

is the Bloch operator with periodic potential .  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. , then the dominant term of the wavefunction exhibits a periodic behavior with Bloch period within an interval with amplitude , where is the width of the first band and where is the strength of the external homogeneous field (in the case of then takes only positive values, obviously).  Therefore, for times of the order of the Bloch period 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 (with a fixed number of wells) and that the Stark potential is replaced by a Stark-type potential due to an homogeneous external field which acts only in a bounded region containing the 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)

 {iϵ∂tψ=HNψ+fWN(x)ψ+γ|ψ|2ψ, HN=−ϵ2∂2xx+VNψ(x,0)=ψ0(x) (4)

where 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 ).  We assume that the wells have all the same shape and we denote by 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 localized on the th cell of the lattice; that is

 ψ(x,t)∼N∑ℓ=1cℓ(t)uℓ(x).

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 .  Eventually, may coincide with the Wannier function associated to the first band of the Bloch operator or with the semiclassical single well ground state eigenfunction .  By means of such an approach the unknown functions turn out to be the solutions of a system of time-dependent equations which dominant terms are given by (here we denote )

 iϵ˙cℓ=−λDcℓ−β(cℓ+1+cℓ−1)+γ∥u0∥4L4|cℓ|2cℓ+fbℓcℓ, ℓ=1,…,N (5)

where 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 is equal to the amplitude of the first band [25].  In (5) we’ll fix .  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 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 -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 -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 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 .  In particular, when is a symmetric wavefunction then the oscillating period is almost constant for small and it practically coincides with the Bloch period ; on the other hand when is an asymmetrical function the oscillating period actually depends on .  This fact is in contradiction with the Bloch Theorem (which holds true when ), which implies that the Bloch period 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 for times of the order of the Bloch period 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 -mode 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.

### Notation

Let be a quantity depending on the semiclassical parameter .  In the following

 g=~O(e−S0/ϵ)

means that for any and any there exists such that

 |g|≤Ce−(S0−ρ)/ϵ, ∀ϵ∈(0,ϵ⋆).

Hereafter, by we denote a generic positive constant independent of .

Let , then by we denote the set of first positive integer numbers.

By we denote the norm of the Banach space , by we denote the scalar product of the Hilbert space .

## 2. Derivation of the DNLS (5)

### 2.1. Assumptions and preliminary results

We consider the time-dependent non-linear Schrödinger equation (4) where is a multiple-well potential and is a bounded Stark-type potential.  In particular we assume that

###### Hypothesis 1.

Let be an even (i.e. ) smooth function with compact support with a non degenerate minimum value at :

 v(x)>vmin=v(0), ∀x∈R, x≠0.

The multiple-well potential is defined as

 VN(x)=N∑ℓ=1v(x−xℓ)

for some fixed , where and where is such that .

Hence, by construction the potential has exactly wells with not degenerate minima at , .

###### Remark 1.

We assume that is an even function just for argument’s sake.  As discussed in Remark 6 this assumption may be removed.  Furthermore, we assume that is a smooth function as usual; in fact, a lessere regularity (e.g. ) would be enough.

###### Hypothesis 2.

Let be the monotone not decreasing function defined as

 WN(x)=⎧⎪⎨⎪⎩−L if x<−Lx if x∈[−L,L]L if x>L

for some .

That is the Stark-type potential is linear in the region containing the wells and it is a constant function outside this region (see Fig. 1).  In the “limit” where goes to infinity the potential becomes a periodic potential with period and the external potential becomes the Stark potential .

###### Remark 2.

We restrict ourselves to a multiple-well potential 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 is not merely for the sake of simplicity; actually, the Stark-type potential is a bounded operator while the Stark potential 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 are small for small.  That is we assume that there exists such that

 Ce−(S0−ρ)/ϵ≤|f|≤Cϵs, ∀ϵ∈(0,ϵ⋆),

for some , and independent of ; furthermore, we assume also that

 |γ|ϵ−1/2|f|≤C (6)

for some positive constant and for any .

The self-adjoint extension of the linear Schrödinger operator formally defined on as

 HN=−ϵ2∂2xx+VN

has an almost degenerate ground state with dimension .  More precisely, let , , be the lowest eigenvalues of with associated normalized eigenvectors .  In particular we have that (see Lemma 2 [25])

 λℓ=λD−2βcos(ℓπN+1)+O(ϵ∞)e−S0/ϵ, ℓ∈NN,

where

 S0=∫x1x0√VN(x)−vmindx>0

is the Agmon distance between two wells and is the ground state of the single well operator , where the single well potential has been introduced by Hyp. 1.  The numerical pre-factor is the hopping matrix element between neighboring wells, and it is such that is asymptotic to the amplitude of the first band of the periodic Bloch operator ; i.e. where and are, respectively, the bottom and the top of the first band.  Such a numerical pre-factor is going to be exponentially small, i.e.

 β=~O(e−S0/ϵ)  as  ϵ→0+.
###### Remark 3.

Hyp. 3 means that, from a practical point of view, the parameter 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]

 vℓ=N∑j=1αℓ,juscj+O(ϵ∞)e−S0/ϵ

where

 αj,ℓ=αℓ,j=√2N+1sin(jℓπN+1)

and where is the semiclassical single well ground state eigenfunction localized on the -th cell; by construction and since is an even function then

 uscj(x)=usc0(x−xj)  and  usc0(x)=usc0(−x). (7)

Now, let be the projection operator associated with the eigenvalues , i.e.

 Π=N∑ℓ=1⟨vℓ,⋅⟩vℓ

and let

 Πc=1−Π.

Let be the -dimensional space spanned by the eigenvectors , .

###### Remark 4.

Let be the spectrum of ; then it is a well known semiclassical result that

 C−1ϵ≤\rm dist({λℓ}Nℓ=1,σ(HN)∖{λℓ}Nℓ=1)≤Cϵ

for some positive constant .  Hence, since is a self-adjoint operator then

 ∥∥[HN−λD]−1Πc∥∥L(L2→L2)≤Cϵ−1

for some .

###### Remark 5.

By [14] it has been proved that there exists a suitable orthonormal base , , of the space .  The functions are practically localized on the -th well.  More precisely, they are such that

• for any and any ;

• for any ;

• , , and for any ;

• The matrix with elements can be written as

 (⟨uℓ,HNuj⟩)=λD% \sc 1N−βT+DN

where is the tridiagonal Toeplix matrix such that

 Tj,ℓ={0 if  |j−ℓ|≠11 if  |j−ℓ|=1

and where the remainder term is a bounded linear operator from to with bound

 Missing or unrecognized delimiter for \right

for some .

We finally assume that the initial state is prepared on the first “ground states”.  That is

###### Hypothesis 4.

.

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]

 E(ψ)=⟨ψ,HNψ⟩+12γ∥ψ∥4L4+f⟨ψ,WNψ⟩

easily follow:

 ∥ψ(⋅,t)∥L2=∥ψ0(⋅)∥L2  and%  E(ψ(⋅,t))=E(ψ0(⋅)).

Furthermore the following a priori estimate follows, too.

###### Lemma 1.

There exists a positive constant such that

 ∥ψ∥H1≤Cϵ−1/2  and  ∥ψ∥pLp≤Cϵ−p−24,∀p∈[2,+∞].
###### Proof.

Indeed, from Theorem 2 by [24] and its remarks it follows that

 ∥∇ψ∥L2≤C√Λ  and  ∥ψ∥Lp≤CΛp−24p

for some and small enough, where

 Λ=E(ψ0)−Vminϵ2

and where .  In particular, since for some , because is fixed, and since then ; therefore

 ∥∇ψ∥L2≤Cϵ−1/2  and  ∥ψ∥Lp≤Cϵ−p−24p.

Hence, the global well-posedness of the NLS follows [6, 7].

### 2.2. N-mode approximation

Let be the normalized solution of the NLS equation written in the formula

 ψ=ψ1+ψc, ψ1=Πψ=N∑ℓ=1cℓuℓ  and  ψc=Πcψ, (8)

for some complex-valued functions .  By substituting (8) into the NLS (4) then it takes the formula

 {iϵ∂tcℓ=⟨uℓ,HNψ1⟩+γ⟨uℓ,|ψ(⋅,τ)|2ψ(⋅,τ)⟩+f⟨uℓ,WNψ(⋅,τ)⟩iϵ∂tψc=HNψc+γΠc|ψ(⋅,τ)|2ψ(⋅,τ)+fΠcWNψ(⋅,τ) (9)

We are going now to get an a priori estimate of the remainder term .  First of all we rescale the time and we redefine the wavefunction up to a gauge factor .  The Bloch period becomes

 τB=βϵT=2πβ|f|b

Hence, (9) becomes (where denotes the derivative with respect to )

 {iβc′ℓ=⟨uℓ,(HN−λD)ψ1⟩+γ⟨uℓ,|ψ(⋅,t)|2ψ(⋅,t)⟩+f⟨uℓ,WNψ(⋅,t)⟩iβψ′c=(HN−λD)ψc+γΠc|ψ(⋅,t)|2ψ(⋅,t)+fΠcWNψ(⋅,t) (10)
###### Theorem 1.

Let Hyp.1-4 be satisfied; then it follows that the remainder can be estimated for times of order of the Bloch period.  That is for any fixed it follows that

 maxτ∈[0,MτB]∥ψc(⋅,τ)∥L2≤C|f|ϵ

for some positive constant .

###### Proof.

From the first equation of (10) and recalling that

 N∑ℓ=1|cℓ(τ)|2=∥ψ1∥2L2=1−∥ψc∥2L2≤1

then a priori estimate follows

 |c′ℓ| ≤ ⟨uℓ,(HN−λD)ψ1⟩β+|γ|β∥ψ∥2L∞+|f|β∥WN∥L∞ (11) ≤ C+|γ|ϵ−1/2β+|f|βL

because and , and from Remark 5 iv. and Lemma 1.

Concerning it satisfies to the following integral equation

 ψc=I+II

where we set

 I := −iβ∫τ0e−i(HN−λD)(τ−s)/βΠcAds II := −iβ∫τ0e−i(HN−λD)(τ−s)/βΠcBds

and where and are defined as

 A := γ|ψ1|2ψ1+fWNψ1 B := Missing or unrecognized delimiter for \right

such that

 A+B=γ|ψ|2ψ+fWNψ.

By means of standard arguments [24] and making use of the fact that the operator is bounded then it follows that

###### Lemma 2.

Let

 Γ=|γ|ϵ−1/2+|f|.

Then the functions and are such that

 ∥A∥L2≤CΓ, ∥B∥L2≤CΓ∥ψc∥L2  and  ∥∥∥∂A∂τ∥∥∥L2≤CΓ2β−1.
###### Proof.

Indeed,

 ∥A∥L2≤|γ|∥ψ1∥2L∞∥ψ1∥L2+|f|∥WN∥L∞∥ψ1∥L2≤C[|γ|ϵ−1/2+|f|]

since .  Similarly, the estimate of the function follows recalling that from Lemma 1.  Finally, the estimate concerning immediately follows from (11). ∎

Hence, the estimates of the integrals I and II follow; in particular, integral II can be simply estimated as

 ∥II∥L2≤CΓβ−1∫τ0∥ψc(⋅,s)∥L2ds

since

 Extra open brace or missing close brace

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 :

 I = [−ie−i(HN−λD)(τ−s)/β[HN−λD]−1ΠcA]τ0+ +i∫τ0e−i(HN−λD)(τ−s)/β[HN−λD]−1Πc∂A∂sds

From this fact and recalling that (Remark 4)

 ∥[HN−λD]−1Πc∥L(L2→L2)≤Cϵ−1

then

 ∥I∥L2≤Cϵ−1maxs∈[0,τ][∥A∥L2+τ∥∥∥∂A∂s∥∥∥L2]≤Cϵ−1Γ[1+Γβ−1τ].

Therefore, we have that

 ∥ψc∥L2≤CΓβ−1∫τ0∥ψc(⋅,s)∥L2ds+Cϵ−1Γ(1+Γβ−1τ).

From the Gronwall’s Lemma it follows that

 ∥ψ(⋅,τ)∥L2≤Cϵ−1Γ(1+Γβ−1τ)eCΓβ−1τ

In particular we observe that

 maxτ∈[0,MτB]∥ψ(⋅,τ)∥L2≤Cϵ−1Γ(1+β−1ΓMτB)eCΓβ−1MτB≤CΓϵ−1≤C|f|ϵ

proving the Theorem since from Hyp. 3 and since . ∎

We are going now to estimate the solutions of the first equation of (10) which can be written as

 iβc′ℓ=⟨uℓ,(HN−λD)ψ1⟩+⟨uℓ,A⟩+⟨uℓ,B⟩

where the term can be represented by property iv. of Remark 5.  Concerning the term the following estimate uniformly holds with respect to the index

 |⟨uℓ,B⟩|≤∥B∥L2≤CΓ∥ψc∥L2.

Furthermore

 ⟨uℓ,A⟩ = γN∑j,k,m=1¯cjckcm⟨uℓ,¯ujukum⟩+fN∑j=1cj⟨uℓ,WNuj⟩ = γ|cℓ|2cℓ∥uℓ∥4L4+fcℓ⟨uℓ,WNuℓ⟩+γraℓ+frbℓ

where

 raℓ=∑j,k,m∈NN : |j−ℓ|+|m−ℓ|+|k−ℓ|>0¯cjckcm⟨uℓ,¯ujukum⟩

and

 rbℓ=∑j,ℓ∈NN, j≠ℓcj⟨uℓ,WNuj⟩

are remainder terms.

We have that

###### Lemma 3.

The following estimates uniformly hold with respect to the indexes and :

• ;

• ;

• where

 r=max[|j−ℓ|,|m−ℓ|,|k−ℓ|,|j−m|,|j−k|,|k−m|].
###### Proof.

Indeed, let , then

 ⟨uℓ,WNuℓ⟩=∫Iℓ|uℓ(x)|2xdx+∫R∖Iℓ|uℓ(x)|2WN(x)dx

where from (7) and Remark 4.  Therefore

 ∫Iℓ|uℓ(x)|2xdx = ℓb∫I0|u0(x)|2dx+∫I0|u0(x)|2xdx+~O(e−S0/ϵ) = ℓb[∥u0∥2L2−∥u0∥2L2(R∖I0)]+~O(e−S0/ϵ)

where is normalized and because .  From this fact and since (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

 |⟨uℓ,WNuj⟩|≤∥WN∥L∞∥uℓuj∥L1=~O(e−S0|j−ℓ|/ϵ)

and, where we assume that ,

 |⟨uℓ,¯ujumuk⟩|≤∥um∥L∞∥u∥L∞∥uℓuj∥L1=~O(e−S0r/ϵ)

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

 iβc′ℓ=−βN∑m=1Tℓ,mcm+γ∥uℓ∥4L4|cℓ|2cℓ+fℓbcℓ+~O(e−S0/ϵ) (12)

where

 ∥uℓ∥4L4=∥u0∥4L4+~O(e−S0/ϵ)

and where the remainder terms are uniform with respect to the index .

###### Remark 6.

In fact, if is not an even function then by means of standard semiclassical arguments it follows that property Lemma 3 i. becomes

 ⟨uℓ,WNuℓ⟩=ℓb+c+~O(e−S0/ϵ)

for some constant independent of the index .  In such a case we must add the term to the right hand side of the DNLS above and, by means of a gauge choice , we can remove this term obtaining again equation (12).

Now, we are able to prove that

###### Theorem 2.

Let be the solutions of the DNLS

 iβd′ℓ=−βN∑m=1Tℓ,mdm+γ∥usc0∥4L4|dℓ|2dℓ+fℓbdℓ (13)

satisfying to the initial conditions , where and are the solutions of (10).  Then, for any fixed it follows that

 maxτ∈[0,MτB],ℓ=1,2,…,N|cℓ(τ)−dℓ(τ)|=~O(e−S0/ϵ)  as  ϵ→0.
###### Proof.

The proof is a simply consequence of equation (12) and from the fact that and Hyp.3. ∎

## 3. Numerical analysis of a real model

We consider the experiment where a cloud of ultracold Strontium atoms are trapped in a one-dimensional optical lattice with potential (2).  Realistic data for the experiment are [20]:

• Lattice period: , ;

• Lattice potential depth: where is the photon recoil energy and where is between and ;

• Mass of the strontium isotope: ;

• Effective one-dimensional nonlinearity strength: let be the effective nonlinearity strength for the three-dimensional Gross-Pitaevskii equation, then it is expected that the effective one-dimensional nonlinearity strength is of the order [26]

 γ≈γ3D2πd2⊥

where is the oscillator length of the transverse confinement; here denotes the scattering length of the Strontium isotope: , where is the Bohr radius; is the number of atoms of the condensate; in typical experiments and ;

• Acceleration constant .

The confined BEC is governed by Eq. (4) and here we make use of the -mode approximation (13), that is the wavefunction has the form where are the solutions of (13) and where 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 semi-classical approximation agrees or not with the Wannier function .  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 too, where it has been numerically verified that these coefficients are negligible when for ; thus, for such values of it is admitted that the -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 of the wavefunction has dominant behavior

 usc0(x)=(mμ)1/8(πℏ)1/4e−√mμx2/2ℏ (14)

in the semiclassical limit, where