Superfluid dynamics of a Bose–Einstein condensate in a periodic potential

We investigate the superfluid properties of a Bose–Einstein condensate (BEC) trapped in a one-dimensional periodic potential. We study, both analytically (in the tight binding limit) and numerically, the Bloch chemical potential, the Bloch energy and the Bogoliubov dispersion relation, and we introduce two different, density dependent, effective masses and group velocities. The Bogoliubov spectrum predicts the existence of sound waves, and the arising of energetic and dynamical instabilities at critical values of the BEC quasi-momentum which dramatically affect its coherence properties. We investigate the dependence of the dipole and Bloch oscillation frequencies in terms of an effective mass averaged over the density of the condensate. We illustrate our results with several animations obtained solving numerically the time-dependent Gross–Pitaevskii equation.


Introduction
The study of the superfluid properties of Bose-Einstein condensates (BECs) trapped in periodic potentials attracts a fast growing interest. The main reason is that the control parameters of such systems are widely tunable in realistic experiments, allowing for the investigation of different and fundamental issues of quantum mechanics, ranging from quantum phase transitions [1] and atom optics [2,3] to the dynamics of Bloch and Josephson oscillations [4]- [6]. Several efforts are also focused on the realization of new technological devices, such as BEC interferometers working at the Heisenberg limit [3] and quantum information processors [2].
The dynamics of BECs in lattices is highly non-trivial, essentially because of the competition/interplay between the discrete translational invariance of the periodic potential and the nonlinearity arising from the interatomic interactions. For deep enough optical potentials, interactions induce a quantum transition from the superfluid to a Mott-insulator phase [1,7,8].
In this work we study the system in a region of parameters such that its ground state stands deeply in the superfluid phase, with the dynamics governed by the Gross-Pitaevskii equation (GPE). Because of the discrete translational invariance, the excitation spectrum of the system exhibits a band structure which has several analogies with the electron Bloch bands in metals [9]- [11]. On the other hand, the coexistence of Bloch band and nonlinearity allows, for instance, solitonic structures [12]- [14] and dynamical instabilities [15]- [17] which do not have an analogue either in metals or in Galilean invariant systems.
Exact, time-dependent solutions of the GPE with an external periodic potential, equation (1), can be written as Bloch states, namely as plane waves modulated by functions having the same periodicity of the lattice. The dynamics of small amplitude perturbations on top of these states satisfies two coupled, linear Bogoliubov equations, which can be solved numerically. However, when the interwell barriers of the periodic potential are high enough, the system can be described in a nonlinear tight binding approximation and several important properties of the system can be 112.3 retrieved analytically [18]. Indeed, in the nonlinear tight binding approximation the continuous GPE can be replaced with a discrete nonlinear (DNL) equation (equation (5)), where the relevant observables of the system are the number of particles N l (t) trapped in well l and the relative phases φ l,l+1 (t) = φ l+1 (t) − φ l (t). In this paper, we rewrite the results derived in [18] in a more convenient form, namely in terms of two effective masses and group velocities. Furthermore, we compare our analytical expressions with full numerical solutions, and we extend our analysis to investigate the behaviour of the system at low optical potential depths, where the nonlinear tight binding approximation breaks down. We show that the phenomena predicted by the DNL equation (5) can be generalized to the case of shallow potentials, bringing new insights into the dynamics of the system.

Discrete nonlinear dynamics
In the 'classical' (mean field) approximation, the BEC dynamics at T = 0 is governed by the GPE [19] ih where g = 4πh 2 a/m, with m the atomic mass and a the s-wave scattering length: a > 0 (a < 0) corresponds to an effective interatomic repulsion (attraction). In the following we consider only a BEC with repulsive interatomic interactions. The external potential V ext includes the optical periodic potential V P , which is typically superimposed on a harmonic (or linear) potential V M . The periodic potential is where d is the lattice spacing and π/d is the wavevector of the lasers in the lattice direction. The lattice spacing determines the Bragg momentum corresponding to the boundary of the first Brillouin zone. The energy barrier between adjacent sites is expressed in units of the recoil energy E R = q 2 B /2m. From (2) we see that the minima of the laser potential are located at the positions x l = ld (l is an integer). Around these points, y y 2 + ω 2 z z 2 ]/2 includes the transverse confining field, and the 'driving' field V D = mω 2 x x 2 /2 gives the effective force acting on the centre of mass of the condensate wavepacket.
In order to understand the basic physics of the system, we first consider the case of deep optical lattices, where analytic solutions can be obtained in the tight binding approximation. Then we study the behaviour of the system beyond the tight binding limit, solving numerically the Gross-Pitaevskii and Bogoliubov equations with arbitrarily shallow periodic potentials.

112.4
As has been shown in [18], when the interwell barriers are much higher than the chemical potential, it is possible to write the condensate wavefunction as where the Wannier wavefunctions l are well localized in each well. The total number of atoms is N T = l N l ≡ l |ψ l | 2 . Replacing this ansatz in the GPE (1) and integrating out the spatial degrees of freedom, we find the DNL equation The 'local' chemical potential µ loc l is the sum of three contributions which depend on the atom number N l explicitly through |ψ l | 2 and implicitly through the shape of the l . The tunnelling rates K l,l±1 between the adjacent sites l and l ± 1 are where the onsite wavefunctions have been calculated with an average number of atoms per site, N 0 = |ψ 0 | 2 , namely l (N l ) ¯ l (N 0 ) (a discussion of the validity of this approximation is in [18]). On the same line, the coefficient χ is given by and the on-site energies V l , arising from any external potential superimposed on the optical lattice, are such that V l ∝ l 2 (V l ∝ l) when the driving field is harmonic (linear). The dependence of the local chemical potential on the number of atoms is affected by the effective dimensionality of the condensates trapped in each well of the lattice. This can be determined by comparing the interaction chemical potential µ int l = |ψ l | 2 g d r 4 l and the three frequencies,ω x , ω y , ω z obtained expanding the lattice potential around the minima of each well V L m[ω 2 x (x − x l ) 2 + ω 2 y y 2 + ω 2 z z 2 ]/2. For instance, whenω x , ω y , ω z µ int l , the spatial width of each trapped condensate does not depend (in the first approximation) on the number of particles N l in the same well, and the condensate wavefunction in each valley of the periodic potential is well approximated by a Gaussian. We consider this as a 0D (zero-dimensional) case: then, the nonlinear tight binding approximation (4) reduces to the usual tight binding approximation ( r , t) = l ψ l (t) l ( r) [12]. The 1D case arises when two frequencies are larger than the interaction chemical potential. In this case the system realizes an array of weakly coupled cigar-shaped condensates, with l factorized as Gaussians along the two tight directions and a Thomas-Fermi in the other direction. In the 2D case only one frequency is smaller than 112.5 the local interaction chemical potential: we have an array of pancake-like condensates, with l factorized as a Gaussian along the tight direction and a Thomas-Fermi in the two other directions. The 3D case is given by the condition µ int l ω x , ω y , ω z and the wavefunction in the lth well l is simply given by a three-dimensional Thomas-Fermi function. The crucial point is that the effective dimensionality of the condensates gives a different scaling of the local interaction chemical potential (6) with the number of atoms with α = 4/(2 + D), where D = 0, 1, 2, 3, |ψ l | 2 is the number of atoms in well l, and U α is a constant which does not depend on the number of atoms nor on the site index. When χ N 0 K and D = 0, the DNL equation (5) gives the discrete nonlinear Schrödinger equation [12].

Excitation spectra
In this section we derive the Bloch and the Bogoliubov excitation spectra of the system in the absence of any driving field (V l = 0). First we derive our results analytically in the tight binding approximation; then we solve the equations numerically for a wide set of parameters to extend our treatment beyond the tight binding regime.

Bloch energy, Bloch chemical potential, effective masses and group velocities
The Bloch states p (x) = e i px/h˜ p (x), where˜ p (x) is periodic with period d, are exact stationary solutions of the GPE (1). Both the energy per particle ε α ( p) (Bloch energy) and the chemical potential µ α ( p) of such solutions form a band structure, so that they can be labelled by the quasi-momentum p and the band index α.
The DNL (5) describes only the lowest band of the spectrum (in the following, we will consider only the lowest band α = 1, and we will omit, for simplicity, the band index α). Exact solutions of the DNL are the 'plane waves' ψ l = ψ 0 exp(i( pld − µt)/h), where p is the quasimomentum, and l is the site index (note that the ψ l are plane waves in the discrete l-space, but do not correspond to plane waves in real space). Within the DNL equation framework, the energy per particle ε( p) and chemical potential µ( p) corresponding to these solutions are where µ loc = µ loc l | ψ l =ψ 0 = ∂(N 0 ε loc )/∂ N 0 , with N 0 = |ψ 0 | 2 the number of atoms per well and ε loc = 2U α N α/2 0 /(α + 2). In the previous equations we have introduced the effective masses m ε and m µ , to emphasize the low momenta (long wavelength) quadratic behaviour of the Bloch energy spectrum and of the chemical potential [20]. It turns out that several dynamical properties of the system can be intuitively understood in terms of such effective masses. This approach is quite common, for instance, in the theory of metals, where m µ ≡ m ε . However in BEC, because of the nonlinearity of the GPE, the two relevant energies of the system, ε and µ, have the same cos(π p/q B ) dependence on the quasi-momentum p, but different curvatures. Therefore, where K and χ have been defined in equations (7) and (8). Sometimes it is convenient to extend the definition of the effective masses to the full Brillouin zone, introducing the quasi-momentum dependent masses where, following equations (13) and (14), m ε = m ε (0) and m µ = m µ (0). It is also useful to introduce, with the same line of reasoning, two different group velocities, defined as There is a simple, general relation between the two different group velocities (following from with, given equations (13) and (14), v µ > v ε . The analogue relation for the effective masses has been retrieved in [20]. Of course, the concept of effective mass, defined as the inverse of the curvature of the corresponding spectrum (as that of group velocity, defined as the first derivative) can be extended to shallow optical potentials, where the nonlinear tight binding approximation breaks down. In this case, the quasi-momentum dependence of ε and µ will not be simply described by a cosine function, but will still remain periodic in the quasi-momentum p. In particular, the value p where m ε ( p) changes sign (corresponding to ∂ 2 ε/∂ p 2 = 0) will be greater than q B /2 and will in general not coincide with the momentum where m µ ( p) changes sign (corresponding to ∂ 2 µ/∂ p 2 = 0). The presence of the two different effective masses (group velocities) raises an important problem: which effective mass (group velocity), and how, enters in the dynamical properties of the system? For instance, we anticipate that the current carried by a Bloch wave with quasimomentum p is ρ 0 v ε ( p), where ρ 0 is the average particle density; m µ , on the other hand, plays a crucial role in the Bogoliubov spectrum. To conclude this section, we remark that Bloch states are not the only stationary solutions of the GPE. Because of nonlinearity, indeed, periodic solitonic solutions can also appear for a weak enough periodic potential, introducing new branches in the excitation spectra [21,22]. 112.7

The Bogoliubov dispersion relation
In this section we study the Bogoliubov spectrum of elementary excitations. This describes the energy of small perturbations with quasi-momentum q on top of a macroscopically populated state with quasi-momentum p (stationary solution of equation (1)). To be explicit, let us consider first the case in which the radial degrees of freedom y, z are integrated out. The wavefunction in the x direction can be written as Because of the periodicity, the Bogoliubov amplitudes can be written as Bloch waves where q is the quasi-momentum of the excitation and {ũ,ṽ} pq (x) are periodic functions. The subscript { pq} indicates that both the amplitudes u,ṽ and the excitation frequencies ω pq depend on the quasi-momentum p of the carrying wave and on the quasi-momentum q of the excitation.
In terms of the periodic functions˜ ,ũ andṽ, the Bogoliubov equations take the form where n is the 3D-average density and (21) and (22) can be solved numerically in a very efficient way working with the Fourier components of˜ ,ũ andṽ.
In the tight binding regime, the Bogoliubov analysis corresponds to perturbing the large amplitude wave as ψ l = [ψ 0 +δψ l ] exp(i( pld −µt)/h), with δψ l = q U q exp(i(qld/h−ω pq t)). Retaining only first-order terms with respect to δψ, we get two coupled linear equations analogous to (21) whose eigenvalues can be calculated analytically [18]. The general solution (for any effective dimensionality of the system: D = 0, 1, 2, 3) is with the chemical potential given by µ( p, N 0 ) = µ loc − (q 2 B /π 2 m µ ) cos(π p/q B ) (see equation (12)), and µ loc l = U α |ψ l | 4/(2+D) . For D = 0 and in the limit χ = 0, we recover the well-known results for the DNL Schrödinger equation [16,24]. Equation (23) has been first written in [18] in terms of the parameter of the DNL equation (5), while, for small q and arbitrary p, has been derived in [17] for arbitrary values of s.
In figure 1, we compare the analytic results (dots) with the numerical solution of equations (21) and (22) (solid curve), for a system in the tight binding regime. In the numerical  (21) and (22) for s = 20 and gn = 0.5E R (green dots); analytic solution equation (23) in the tight binding approximation (solid blue curve); analytic solution of equation (23) where m µ is replaced with m ε (dashed red curve). The quasi-momentum of the carrying wave is calculations, the effective masses are obtained from the curvatures of the Bloch energy and chemical potential spectra, while the term N 0 ∂µ/∂ N 0 has been evaluated from the density dependence of the chemical potential. As was noted in [20], effects related to the difference between the two effective masses in the Bogoliubov spectrum of a condensate at rest ( p = 0) are usually negligible. In contrast, such difference becomes important when the condensate moves with a large quasi-momentum, as shown in figure 1(b).

Sound waves and instabilities
The small q (large wavelength) limit of the Bogoliubov dispersion relation becomes (we assume, for the moment, that 1/m ε (∂µ/∂ N 0 )N 0 cos(π p/q B ) > 0). The linear behaviour in q indicates that the system supports (low amplitude) sound waves, propagating on top of the large amplitude travelling wave p with velocity where the 'chemical potential group velocity' v µ has been defined in equation (18), and the 'relative sound velocity' c is defined as The two different velocities v s,± correspond, respectively, to a sound wave propagating in the same and in the opposite direction of the large amplitude travelling wave. As we have already noted, and discuss again in the next section, v µ is different from (it is larger than) the actual velocity of the large amplitude wave, see equation (19). Contrary to the case of a Galilean invariant system (s = 0), the sound velocity depends on the quasi-momentum p. Moreover, v s depends on the effective dimensionality of the condensates, since (from equations (10) and (12)) (∂µ/∂ N 0 )N 0 ∼ αU α N α/2 0 . In the limit α = 2, p → 0 and m ε , m µ → m we get the sound velocity in the uniform case.
The system is energetically unstable if there exist any ω pq < 0. In the limit s = 0, this corresponds to a group velocity larger than the sound velocity (Landau criterion for superfluidity).
When the system has a discrete translational invariance (s > 0) the condition for this instability is obtained from the Bogoliubov excitation spectrum equation (23). Then, we have that the system is not superfluid when ω pq < 0, corresponding to This result should be compared with the well known Landau criteria for a homogeneous system (s = 0), stating that the superfluid is energetically unstable when v 2 > c 2 , v ≡ ∂ε/∂p = ∂µ/∂p being the group velocity of the condensate, and c = √ 1/m(∂µ/∂ N 0 )N 0 the sound velocity. There is a further dynamical (modulational) instability mechanism associated with the appearance of an imaginary component in the Bogoliubov frequencies, which disappears in the absence of interatomic interactions, or in the translational invariant limit (if a > 0). The onset of this instability in the tight binding regime, coincides with the condition The dynamical instability drives an exponentially fast increase of the amplitude of the (initially small) fluctuations of the condensate. Since the initial phases and amplitudes of the fluctuation modes are essentially random, their growth induces a strong dephasing of the condensate, and dissipates its translational kinetic energy (which is transformed in incoherent collective and single particle excitations). The unstable modes q grow with a timescale given by the imaginary part of the excitation frequency We remark here on the different scaling of the energetic and dynamical instability with the interatomic interactions. With a decrease of the scattering length, the sound velocity decreases, and smaller and smaller group velocities can break down the superfluidity of the system (in the limit a = 0, the sound velocity c = 0: the non-interacting condensate is always energetically unstable for an arbitrary small group velocity). On the other hand, the dynamical modulational instability criterion does not depend on the scattering length. This apparent paradox is simply solved by noticing that the growth time of the unstable modes, equation (29), actually depends 112.10 on interactions, and diverges when the scattering length vanishes (τ → ∞ when a → 0). Therefore, a non-interacting condensate is always dynamically stable. There is a further point to note: if we consider a condensate moving with an increasing velocity, the system always becomes first energetically unstable, then it hits the dynamical instability. However, in real experiments the energetic instability can grow quite slowly (and at zero temperature only in the presence of impurities [15]), so that the dominant dephasing mechanism is given by the modulational instability. This aspect can also be highlighted with numerical experiments, studying, for instance, Bloch oscillations of a condensate with the interactions switched off. In this case, even though the system is energetically unstable, it remains coherent over many oscillations. If the interatomic interactions are switched on, however, the system dephases rather quickly, the dephasing occurring when the quasi-momentum of the condensate is in the dynamically unstable region of the Bloch spectrum. We have done such a numerical experiment, and the results are shown in figures 4 and 5. Of course, our prediction can be tested in real experiments, tuning the scattering length with a Feshbach resonance.
To summarize, the tight binding approximation predicts the onset of dynamical instability (complex excitation frequency) at p = q B /2. We point out that p = q B /2 also corresponds to the quasi-momentum where, in the tight binding regime, the effective masses m ε ( p) and m µ ( p) change sign. A system with a negative effective mass and positive scattering length can be, roughly speaking, seen as equivalent to a system with a negative scattering length and positive effective mass. It is well known that a BEC having a negative scattering length is dynamically unstable, and such parallelism might be thought to give a simple explanation of the instability. However, we will see that this coincidence between the onset of dynamical instabilities and the inversion of sign of the effective mass does not take place at lower optical potentials (see figure 3).
Let us concentrate now on the behaviour of the excitation frequencies for shallow optical potentials, where the tight binding expression derived in (23) is not applicable. For small optical potential depths, the Bogoliubov equations have to be solved numerically and the results show a more complicated behaviour. In figures 2(a)-(c), we show the numerical solutions of the 1D Bogoliubov equations for three different values of s (s = 1, 2 and 5); we plot the real and imaginary part of ω pq as a function of q and we vary p in time.
We point out a series of differences with respect to the tight binding regime: • the complex frequencies appear at the boundary of the first Brillouin zone (q = q B ) for a value of p > q B /2 (dots in figure 3) and they reach the centre of the zone (q = 0) for a higher value of p (orange region in figure 3); • the range of momenta p where the frequencies are complex for some q, but real around q = 0, decrease with increasing s; in the tight binding limit this range vanishes; • in the limit of our numerical accuracy, which is due to the discrete sampling of p and q, we found that the value of p where the effective mass m ε changes sign (squares in figure 3) corresponds to the value of p at which the frequencies with non-vanishing imaginary part reach q ≈ 0. In the tight binding approximation, this appears explicitly through the term cos(π p/q B )/m ε under the square root.
We would like to remark on two important results arising from our study of the excitation spectra. First, as shown in figure 3, we found the onset of the dynamical instability for values of the quasi-momentum where the effective mass m ε ( p) is still positive. Second, the range of momenta where the system has a positive effective mass and, at the same time, is dynamically unstable, increases with decreasing s (bearing in mind that the amplitude of the imaginary part vanishes for s → 0 or gn → 0, which implies that the growth in time of the instability diverges both for uniform interacting systems and for ideal gases in optical lattices). So, one can study the behaviour of the system at low s to distinguish between two possible dephasing mechanisms, one due to the sign of m ε , the other one due to the dynamical instability, as will be detailed in section 6. Various important aspects of the physics of energetic and dynamical instabilities of a BEC in a periodic potential have been studied in [12,15]- [17].

Newtonian dynamics
Using the results in [18], we can now rewrite the dynamics of a BEC wavepacket in terms of the energy effective mass. For the BEC wavepacket we use the following ansatz:  (21) and (22) for gn = 0.5E R as a function of s. Dots: value of the quasi-momentum p where the excitation frequencies ω pq (q = ±q B ) become complex; orange region: value of the quasi-momentum p where complex frequencies are found around q = 0; squares: value of the quasi-momentum p where the effective mass m ε changes sign. The dotted curves are a guide to the eye.
where ξ(t) and σ (t) are, respectively, the centre and the width (in lattice units) of the wavepacket, P(t) and δ(t) their associated momenta and K(σ ) a normalization factor such that l N l = N T (with |ψ l | 2 ≡ N l ). The function f is generic: for instance, we can choose f (X ) = e −X 2 or f (X ) = (1 − X 2 ) 1/α (with −1 ≤ X ≤ 1) to describe, respectively, the dynamics of a Gaussian or a Thomas-Fermi wavepacket. The equations of motion of the collective variables ξ(t), σ (t), p(t), δ(t) have been obtained in [12,18]. With V j = j 2 ( = md 2 ω 2 x /2), and neglecting the dynamics wavepacket width dynamics (σ (t) = 0), we find that the group velocitẏ ξ and the effective force acting on the centre of mass of the wavepacket are given bȳ where V d = (ξ 2 + σ 2 I 2 I 1 ) with I 1 = d X f 2 (X ) and I 2 = d X X 2 f 2 (X ). Since the effective masses depend on the local (on-site) density, we have to introduce an effective mass averaged over the local density of the condensate wavepacket with, according to equation (13) 112.14 (i) in the case of a homogeneous system (V D = 0, N l = constant), the tunnelling rate is given by (ii) the frequency of small amplitude oscillations of the wavepacket driven by a harmonic field V l ∝ l 2 is (iii) if the driving field is linear V l = mGdl, we have simple Bloch oscillations with This analysis does not take into account possible dephasing mechanisms such as those investigated in the previous section. In the collective coordinate approach, such dephasing mechanisms can be described by including the dynamics of the width of the wavepacket σ (t) and of the corresponding momentum δ(t) [12].

Numerical experiments on Bloch oscillations, dipole oscillations and free expansions in the lattice
In this section, we discuss some numerical simulations of the GPE in order to illustrate the phenomena described in the previous sections. We first consider Bloch oscillations (section 6.1): we create a condensate in a harmonic trap superimposed on the lattice and then switch off the harmonic trap and replace it with a linear potential. We expect the BEC to oscillate periodically in space (Bloch oscillations). The second numerical simulation consists in creating a condensate in a harmonic trap superimposed on the lattice, and suddenly displacing the centre of the harmonic trap (section 6.2). This experiment has already been studied theoretically [16] and performed experimentally in [6,23]. We discuss it again, generalizing the previous results to the case of shallow optical lattices.
The third numerical simulation consists in creating a condensate in a harmonic trap superimposed on the optical lattice and then switching off the harmonic trap in the lattice direction, letting the condensate expand in the periodic potential (section 6.3): for values of the mean-field energy large enough (with a fixed height s of the interwell energy barriers), the wavepacket is self-trapped [12] and spreading of the wavepacket does not occur.
In all cases, for an interacting BEC, we found some form of self-trapping and dynamical instabilities for some values of the depth of the periodic potential or the initial conditions of the BEC wavepacket. For instance, in the dipole oscillation experiment, the condensate may stop on one side of the harmonic potential and be unable to complete the oscillation. It is useful to look at the dynamical evolution of the relative phases of condensates trapped in neighbouring wells and the BEC evolution in momentum space. There is a clear correspondence between the distribution in momentum space and that in quasi-momentum space: the quasi-momentum distribution of the Bloch state p (x) = e i px/h˜ p (x) with quasi-momentum p is δ( p); its momentum distribution is | c δ( p + 2 q B )| 2 , where is integer and where c is the Fourier coefficient of the periodic function˜ p (x). An analogous relation is also valid when the condensate is not in a well-defined Bloch state, but in a superposition of Bloch states of the first band. In this case, the width of the peaks of the momentum distribution will be equal to the width of the quasi-momentum distribution. In the following we will work with the momentum distribution, which is simply obtained from the Fourier transform of the condensate wavefunction.

Bloch oscillations
Bloch oscillations can be explained in very simple terms. In the presence of a linear potential superimposed on the optical lattice, the behaviour of a particle with quasi-momentum p is described by the equation of motion p(t) = Ft, where F is the constant force due to the linear potential. Since the velocity is given by v g = ∂ε( p)/∂ p, when the effective mass is negative, the particle will respond to a positive (negative) force with a negative (positive) acceleration. Since the energy band ε( p) is periodic in p, this will result in periodic oscillations in coordinate and velocity space.
This simple explanation, even neglecting important effects like Landau-Zener tunnelling to a higher band, provides a useful model to interpret experiments with electrons [25], with cold atoms [26] and with BECs [4,5]. However, if interactions in the condensate play a major role, the scenario can change dramatically. First of all, the momentum distribution ( p, t) will not evolve just like ( p(t)), as approximately happens for non-interacting systems, but will also spread. Furthermore, in the presence of interactions, it might happen that the condensate gets dephased and, after a short while, the oscillations stop (see figure 4). For the situation described in there (s = 10, for which the tight binding approximation works well), the dephasing process begins when the centre of the momentum distribution reaches p = q B /2. This point corresponds both to the onset of the dynamical instabilities and the inversion of the effective mass. Since the momentum distribution has a certain width, one might think that the oscillations stop because the sign of the effective mass is not the same for the whole condensate. An alternative interpretation relies on the onset of dynamical instabilities.
In order to highlight the correct interpretation, we study the Bloch dynamics of a noninteracting condensate, which is always dynamically stable. The initial spatial width is chosen in order to get about the same momentum distribution as in the interacting case, in order to have similar effective mass effects. More specifically, since in the interacting case the width of the momentum distribution increases slowly, while in the non-interacting case it remains almost constant, we choose the initial conditions so that the two momentum distributions will be similar at the 'critical point', where p = q B /2.
The direct comparison is shown in figure 5. We observe, in the non-interacting case, regular, perfectly periodic Bloch oscillations, in spite of the finite width of the momentum distribution. This clearly shows that, in the interacting case, the onset of decoherence is due to dynamical instability.

Dipole oscillations
Dipole oscillations consist of motion of a condensate at the bottom of a harmonic trap. The average velocity is periodic in time and the momentum distribution, showing characteristic peaks due to the optical lattice, also oscillates periodically in time at the bottom of the band. During the time evolution, the phase differences between neighbouring condensates remain locked over the whole condensate.
For a given set of parameters corresponding to small displacements, small interactions or small optical potential depths, dipole oscillations remain periodic, with the condensate locked in phase. On the other hand, if we increase one of these quantities, we find that the oscillations get dephased during the time evolution, or even stop before the condensate reaches the bottom of the harmonic potential. For the sake of comparison we display in figures 6(a)-(c) the evolution of the density, of the phase difference and of the momentum distribution for the following sets of parameters: for fixed interactions (gn = 0.01E R ) and harmonic trap (hω x = 0.004E R ), we choose Looking at the phase difference between neighbouring condensates, we find that when the condensate oscillation is interrupted, the phases get scrambled. This corresponds to a randomized flux of atoms which are no longer able to flow coherently down to the potential. The evolution of the momentum distribution suggests that this phenomenon happens when the condensate reaches the instability region, given in the specific cases by p greater than q B /2.
To further explore this interpretation, we choose a shallow optical potential such that there is a broad range of p where the effective mass m ε ( p) is positive and at the same time the system is dynamically unstable (see figure 3). We increase the nonlinear interaction parameter to get a relevant imaginary part of the excitation frequencies, otherwise the timescale where instabilities manifest themselves is too long. In figure 6(d) (lower plot), we indicate with a red dotted line the quasi-momentum where the dynamical instabilities arise and with a green dotted line the quasimomentum where the effective mass changes sign. We actually observe the first signatures of decoherence when the momentum distribution is contained between the two lines, indicating that the decoherence is related to the dynamical instability point. Experimental evidence of dynamical instabilities is reported in [23].

Expansion in the lattice
After creating the condensate in the harmonic trap superimposed on the lattice, we switch off the harmonic trap and let the condensate, which is initially at rest, expand. During the expansion, the current of atoms is from the inside to the outside of the cloud and the phase differences increase, being positive for x < 0 and negative for x > 0. In [12] the occurrence of self-trapping was predicted in the tight binding: when interactions are larger than a critical value, the width of the wavepacket does not continue to increase with time (as for vanishing or small interactions) and the wavepacket remains localized around a few sites. A similar nonlinear self-trapping occurs in the two-site problem [27].
Increasing the interactions, the system enters into the self-trapped regime as shown in figure 7(b). If the interactions are strong enough, we see that after a first stage (whose duration depends on the strength of the interactions) the expansion stops and the condensate evolves as a random flow of atoms between the condensates localized at the bottom of the different potential wells, indicating the onset of a new dynamical instability. In this case, however, a x / d Bogoliubov-like stability analysis is much more problematic because of the non-trivial temporal evolution of the condensate wavepacket. A possible, though approximate, approach is to write down an effective Hamiltonian of the system in terms, for instance, of the collective coordinates introduced in section 5. Such a Hamiltonian would contain a limited number of degree of freedom, making the stability analysis a much easier task. This is the approach followed in [12] to study the dynamics of an expanding condensate in the discrete nonlinear Schrödinger equation framework (with m µ = m ε ). Within this approach one recovers, in a unified framework, the critical values of the parameters for the self-trapping conditions of a wavepacket of finite width initially at rest, and the onset of the modulational instability of a Bloch wave discussed in the previous sections. For instance, considering a Gaussian wavepacket with initial width σ 0 and quasi-momentum p 0 , the collective coordinates approach predicts the onset of self-trapping at a critical value of the parameter = U 2 N T /2K [12]. When cos( p 0 ) > 0, the critical value is c ≈ 2 √ πσ 0 cos( p 0 ); when cos( p 0 ) < 0, the critical value is c ≈ 2 √ π | cos( p 0 )|/σ 0 . When the width of the wavepacket is very large (σ 0 → ∞), c → ∞ if cos( p 0 ) > 0 (and the system is always dynamically stable), while c → 0 if cos( p 0 ) < 0 (and the system is always dynamically unstable), recovering the findings of section 4.
The study of the dynamical instabilities of a condensate trapped in a periodic potential is a rich problem, and deserves further investigation. As we have mentioned, this is connected to the general problem of the stability of a non-stationary state, which also includes for instance the propagation of sound waves in the nonlinear regime. First experimental results on self-trapping with weakly coupled BECs are reported in [28,29].

Conclusions
The Gross-Pitaevskii dynamics of a BEC trapped in a deep periodic potential can be studied in terms of a DNL equation. This mapping allows a clear and intuitive picture of the main dynamical properties of the system, which can be calculated analytically. We have calculated the effective masses of the system, connected to the Bloch energy and chemical potential spectra. We have calculated the Bogoliubov dispersion relation, and studied the sound velocity and the appearance of energetic and dynamical instabilities. We have generalized these concepts to the case of shallow optical lattice, which requires a numerical solution and provides insight into the understanding of the problem. Both in the tight binding limit and in the case of the shallow optical potential, we have investigated in detail the onset of dynamical instability, which seems to be the main mechanism of dephasing of the condensate in Bloch oscillation and dipole oscillation experiments.