Time-dependent currents of 1D bosons in an optical lattice

We analyse the time-dependence of currents in a 1D Bose gas in an optical lattice. For a 1D system, the stability of currents induced by accelerating the lattice exhibits a broad crossover as a function of the magnitude of the acceleration, and the strength of the inter-particle interactions. This differs markedly from mean-field results in higher dimensions. Using the infinite Time Evolving Block Decimation algorithm, we characterise this crossover by making quantitative predictions for the time-dependent behaviour of the currents and their decay rate. We also compute the time-dependence of quasi-condensate fractions which can be measured directly in experiments. We compare our results to calculations based on phase-slip methods, finding agreement with the scaling as the particle density increases, but with significant deviations near unit filling.


Introduction
Exciting progress in experiments with cold atoms in optical lattices [1] has not only paved the way for study of the quantum phases associated with strongly interacting many-body systems [2,3], but also for study of non-equilibrium dynamics in such systems. This is particularly true for transport properties, where the long coherence times associated with the experiments make it possible to gain new insight into phenomena such as spincharge separation [4,5,6], and currents in the presence of impurities and junctions [7,8,9,10] by studying them in a new environment. Particularly in the case of 1D systems, this connection to non-equilibrium dynamics is further strengthened by the recent development of numerical methods based on matrix product states [11], including the time evolving block decimation (TEBD) algorithm [12,13] and the related timedependent density matrix renormalisation group (t-DMRG) methods [14,15]. These make possible the quantitative prediction of time-dependent dynamics for size scales that are typical in experiments, as well as the identification of parameter regimes in which specific phenomena can be observed.
A key characteristic for the transport of bosons in a lattice is the stability of currents in the presence of interactions which can arise in regimes where the system is superfluid. However, this stability depends not only on the current and the interactions [16,17,18,19], but also on the dimensionality of the system. This was demonstrated in recent experiments by observing damping of the centre of mass motion for bosons oscillating in an harmonic trap [20]. For sufficiently small initial displacements, the motion was found to be stable up to a critical lattice depth in higher dimensions, whilst for bosons confined to move along one dimension, decay of the oscillatory motion was observed for arbitrarily small depth. This was explained by the increased role of quantum fluctuations in the 1D system [21,22,23,24,25,26,27]. A similar situation occurs in the case where a homogeneous current is created, e.g., by starting in the lowest energy state, and then accelerating the lattice [28,29,30,31,32] (as shown in figure 1). A mean-field stability diagram as a function of the magnitude of the acceleration and the strength of inter-particle interactions was computed by Altman et al. [33,34], predicting a sharp transition between stable and unstable regions. For bosons confined in 3D, such transitions were observed recently [31]. However, the same experiments indicated strong deviations from this behaviour with a 1D gas, for which a crossover was observed between the two regimes. The stability of the resulting current will depend on the value of this quasi-momentum and the inter-particle interactions.
Here we perform a detailed numerical analysis of the decay of currents for a Bose gas moving uniformly in 1D, making quantitative predictions for the time-dependence of currents that can be measured in an experiment. This analysis is made possible via the use of the infinite TEBD (iTEBD) algorithm [35], which until now has primarily been applied to study the ground states of homogeneous systems. Here it is used to investigate the time-dependence of homogeneous current flow, allowing predictions to be made over longer timescales without boundary effects. We begin from the ground state of the Bose-Hubbard model [2], which describes the lattice gas, and consider an acceleration applied in order to produce a finite current. As a function of the interaction strength in the gas and the magnitude of the acceleration, we observe a broad crossover between regions of stable and unstable current, which we characterise via the variation in the rate of decay of the current for different parameters. We find that the crossover region occurs at significantly smaller values of interaction strength and initial acceleration than are predicted by mean-field methods. We also compare our decay rates to those obtained for systems with large filling factors via phase-slip calculations by Polkovnikov et al. [34].
The rest of this paper is organised as follows: In section 2 we discuss the Bose-Hubbard model for the system we study, and the methods we use to characterise the currents in the system. In section 3 we describe the use of the iTEBD algorithm and the procedure we use in our calculations. In section 4 we present a quantitative analysis of the decay of currents within the system, as well as the dependence of the decay on the initial acceleration and the inter-particle interactions. The scaling properties of the decay rates are investigated in section 5, and compared to calculations based on phase-slip methods, and in section 6 we give a summary and outlook for this work. Appendix A provides detailed information about translationally invariant block sizes within the iTEBD algorithm.

Moving 1D bosons in an optical lattice
In this section we introduce the system of bosons moving in 1D in an optical lattice. We first introduce the Hamiltonian describing the system in 2.1, before discussing the current and general expectations for its stability during time evolution in 2.2. In order to make stronger connections to quantities that can be measured in an experiment, we introduce the quasi-condensate fraction in 2.3, the decay of which is strongly related to the decay of the current.

Bose-Hubbard model
The dynamics of bosons in the lowest band of an optical lattice is described by the Bose-Hubbard model [1,2,3] with Hamiltonian ( ≡ 1): Here,b i (b † i ) are the bosonic annihilation (creation) operators at lattice site i,n i =b † ib i , U is the on-site interaction energy shift, and J is the tunnelling amplitude between neighbouring sites, with ij denoting a sum over neighbouring lattice sites. This model is valid in the regime where Un, J ≪ ω T , where ω T is the energy separation between the lowest two Bloch bands, and we denote the filling factor of N particles on M lattice sites asn ≡ N/M.
For integern, a phase transition is observed in the ground state of this model as a function of u ≡ U/J between superfluid and insulating behaviour [36]. In one dimension withn = 1, the critical value of u is u c ≈ 3.37 [37]. The superfluid (SF) ground state is characterised by quasi off-diagonal long range order in the single particle density matrix (SPDM) b † ibj . For u > u c , the system is in a Mott insulator (MI) phase, with exponentially decaying off-diagonal elements of the SPDM b † ibj . In what follows below, we will assume that the system is prepared in the ground state of the Bose-Hubbard model for a particular choice ofn and u, and then accelerated to produce a finite initial current.

Currents in the Bose-Hubbard model
By accelerating the lattice or applying a linear gradient potential for a short period of time, it is possible to create a current of atoms moving with respect to the optical lattice. This can be quantified via the current operator which appears in the continuity equation Note that the average of the current expectation values M m ĵ m /M is proportional to the mean group velocity calculated via the quasi-momentum distribution of the particles.
In particular, we can define bosonic operators for the quasi-momentum modes,â q , which are related tob i asb m ≡ q e −iqamâ q / √ M , where a is the lattice spacing, and q = 2πr/L for r ∈ (−M/2, M/2] an integer, with L = aM the lattice length. We see that where n q ≡ â † qâq is the quasi-momentum distribution function at quasi-momentum q, for which the corresponding group velocity per lattice constant is given by 2J sin(qa). In this study, the current will always be translationally invariant, and so we omit the site label and write ĵ ≡ ĵ k , for any k.
If a finite current is induced to the system by accelerating the lattice so that the quasi-momentum distribution is shifted by an amount ka (see figure 1), then we expect the resulting behaviour to depend on the superfluidity of the gas. Even in the presence of interactions, a superfluid flow without dissipation can exist. However, if a current is generated in a non-superfluid initial state, or if the flow becomes dynamically unstable, then the current will decay. Whether the current is stable or not is thus a function of both the initial state of the system and the initial mean quasi-momentum, ka [34].
As has been predicted theoretically in the weakly interacting Gross-Pitaevskii regime [18,19] and demonstrated in experiments [17,28,29], the current will be unstable if the initial superfluid ground state (with all atoms near ka = 0) is accelerated to the inverted part of the lowest Bloch band, i.e., when ka > π/2 (this corresponds to a classical instability). We thus expect stable currents only when ka < π/2. If we begin with the ground state of a system with integer fillingn, then we also do not expect stable superfluid currents above the critical interaction u c for the SF-MI phase transition. A stability diagram, interpolating between these two regimes was investigated by Altman et al. [33] using a Gutzwiller mean-field technique, and they found a sharp transition between stable and unstable regions. This sharp transition was observed in an experiment by Mun et al. [31] for bosons allowed to move in the lattice in three dimensions. However, the same experiment observed behaviour more characteristic of a crossover between stable and unstable currents in 1D.
In this sense, in a 1D system it is not possible to characterise clear regions where the currents decay and do not decay, but rather the crossover must be characterised by how rapidly the current decays for different initial mean quasi-momenta and inter-particle interaction strengths. In section 4 we will investigate this crossover quantitatively in terms of such decay rates using the iTEBD algorithm.

The quasi-condensate fraction
In section 4 we show that the decay of the current results from two-particle scattering processes that broaden the quasi-momentum distribution. In 2D or 3D decay of the current is thus also directly linked to a decrease in the condensate fraction in the experiment, which was used as the key experimental observable in reference [31]. In 1D, there is no condensate for an infinite system, however as the system is finite (of the order of 100 occupied lattice sites), a peak that would correspond to the condensate is observed in time of flight measurements of the momentum distribution. In order to make connection to this observable, we will calculate not only the time dependence of the current ĵ below, but also a quasi-condensate fraction defined over a finite portion of the system. This fraction, C R , can be defined as the largest eigenvalue of the reduced SPDM considering only a range of R sites in the system i.e. b † ib j , with i, j ∈ [l, l + R] for some arbitrary site l. We typically use R = 100 below, reflecting the typical occupied number of lattice sites in experiments.

Numerical Calculations of the current decay
In the remainder of this manuscript, we analyse the decay of a current, which is created by suddenly imparting a quasi-momentum ka at time t = 0 to bosons loaded into the lattice in the ground state configuration (which has zero mean quasi-momentum). This resembles the experimental situation of reference [31], where a current is obtained by means of a moving optical lattice which is accelerated to a final quasi-momentum ka on a timescale long enough to ensure adiabaticity with the respect to inter-band transitions.
Here we assume that the acceleration is fast compared with the tunnelling timescale 1/J, resulting in a simple translation of the quasi-momentum distribution, (see section 4.1 and figure 2). This corresponds to an application of the operator to the initial ground state.
To compute both the initial ground state of the system, and the time dependence of the current once the initial mean quasi-momentum ka is imparted on the system, we make use of the iTEBD algorithm [35]. This algorithm is an extension of the TEBD algorithm [12,13], and makes possible the near-exact integration of the Schrödinger equation for 1D lattice Hamiltonians in an infinite, translationally invariant system. These methods have not only been applied to coherent dynamics in 1D, but also generalised to dissipative systems [38,39,40], and extensions to a 2D state ansatz have been recently considered [11].
Until now, iTEBD has been primarily used to compute the ground state of translationally invariant systems, however here this algorithm becomes crucial in computing time-dependent dynamics ‡. By computing the decay of the current in an effectively infinite system, we avoid difficulties arising due to either open boundary conditions, which can prevent propagation of moving bosons, and periodic boundary conditions, which restrict the length over which correlation functions can be computed. In our calculations we apply the iTEBD algorithm as described in reference [35], but with the modification that the length of a single block in the translationally invariant system is increased from two to 2π/ka. Although a matrix product state representation can reproduce phase relationships with any period, even with a single repeated matrix, we find it convenient to increase the block size so as to represent one full cycle of the phase when we apply the operatorK(ka) to the state (see Appendix A for more details).
Throughout our calculations, we performed convergence tests in the matrix product state bond-dimension χ, as well as the time step for integration, and the maximum number of particles allowed per site in our description of the Hilbert space, d. For a more detailed analysis of how errors enter this method, see [42].

Analysis of Current decay
In this section we characterise the current and its decay as a function of the initial mean quasi-momentum ka and the strength of interactions u, as well as the evolution of the characteristics of the state as a function of time. We begin in section 4.1 by analysing the ‡ Time-dependent dynamics with the iTEBD algorithm have also been recently employed to study quantum quenches [41].
initial state of the system at a time t = 0 + , after the instantaneous acceleration of the lattice results in a mean initial quasi-momentum ka. In section 4.2 we then describe the dynamical changes in the quasi-momentum distribution during time evolution. Finally, we discuss the dependence of the decay on ka and u and provide a shaded plot of a stability crossover diagram in section 4.3.

Initial current
In this subsection we characterise the system configuration at a time t = 0 + immediately following the initial instantaneous acceleration (corresponding to an imposition of the initial mean quasi-momentum ka). We consider both the quasi-momentum distribution n q and the corresponding current ĵ .   Figure 2 shows n q before and after the instantaneous application of a momentum shift of ka = 0.4π at t = 0 for a situation where the lattice filling isn = 1. At t = 0 − (that is, before the shift is applied) the system is in a SF or MI configuration, with u = 1 and 4, respectively. At time t = 0 − and for u = 1, n q is strongly peaked around q = 0, as expected for a superfluid ground state, while it broadens as interactions are increased, as shown for u = 4. The shift ka is applied instantaneously, and the shape of n q is unaltered at time t = 0 + . This latter distribution is the initial condition for the subsequent analysis of the system dynamics. Figure 3 shows the current ĵ calculated at t = 0 + for a few values of u. In the limit of vanishing interactions, the current ĵ at t = 0 + is ĵ = 2Jn sin (ka), corresponding  to all particles occupying the same quasi-momentum state. In a MI state with integer n, this current will be zero in the limit of infinitely strong interactions, as the quasimomentum distribution n q becomes flat. In figure 3 the case with u = 0 is plotted as a solid black line. The figure shows that for finite values of the interaction strength u from figure 2, ĵ decreases from the value 2Jn sin (ka) as expected. However, we note that the current is comparatively large even for the initial condition u = 4, corresponding to a MI state at t = 0 − . This is related to the fact that the quasi-momentum distribution n q , whilst broader than in the superfluid case, is still finite in width (see figure 2), and far from the flat distribution corresponding to the limit u → ∞. As u is increased, the distribution continues to broaden gradually, but there is no strong qualitative change in the initial quasi-momentum distribution or the initial current values on entering the Mott Insulator regime.

Time evolution of n q
Due to the lack of Galilean invariance in a lattice, a state with a finite current does not correspond to an eigenstate of the system. We expect that in the presence of interactions any initial current ĵ will eventually decay to a current-free state. The dynamics of the current decay corresponds to the rearrangement of the quasi-momentum distribution by two-atom scattering processes. The quasi-momentum of two particles is always conserved in these scattering processes, but because of Umklapp processes that connect the two edges of the Brillouin zone, the mean quasi-momentum is significantly changed, and will reach zero in a final steady state.
In figure 4 we show results for the time evolution of n q for a system with the initial conditions ka = 0.4π and u =n = 1. We observe for increasing time that the width of the momentum distribution increases and its mean value shifts towards zero. We have checked that the off-diagonal elements of the SPDM decay exponentially in the final state. These results are consistent with those reported in reference [24], where the change in the redistribution of quasi-momentum in the final state is linked to the appearance of a finite temperature in the system [43,44,45,46].  Time evolution of the normalised ground state quasi-momentum distribution n q in the SF regime for u = 1 after a momentum ka = 0.4π is imparted. The particles are redistributed into a final state with a broader distribution at smaller mean momentum on the timescale 10/J (see text). The numerical parameters are χ = 100, d = 6.
We find the behaviour of n q and SPDM described above to be typical of all parameter values we have investigated, with the primary difference for different parameters being in the timescale on which these processes occur. This includes the Mott Insulator regime, where the main qualitative difference is in the timescale of the decay of the current, not in the form in which the quasi-momentum distribution is redistributed. The dependence of this timescale on the different system parameters will be investigated below.

Stability diagram for the current
We now investigate the decay of the current as a function of the initial mean quasimomentum ka and the inter-particle interactions u for an initial ground state with n = 1. We focus on the time-dependent dynamics of ĵ and also of the condensate fraction C 100 introduced in section 2.3. We extract a stability diagram for the crossover between parameter regions where currents are stable over long time periods, and regimes where currents decay rapidly.
In the following we focus on systems with interaction strengths 0.5 ≤ u ≤ 4.0 and 0.05π ≤ ka ≤ 0.50π, and we compute the real time evolution of the boson current and the condensate fraction C 100 for times 0 ≤ tJ ≤ 10. This time interval is sufficiently short to allow currents to be probed in the experiment within typical decoherence times, e.g., due to incoherent light scattering [2], and also allows for accurate numerical results. Figure 5(a) shows the time evolution of ĵ for a constant initial momentum ka = 0.1π and increasing values of the on-site interaction u, and figure 5(b) shows the decay of ĵ for a fixed on-site interaction of u = 1.0 and various values of ka. On the considered timescale we find that a constant superfluid current (with no decay) only exists for small values of ka and u, for example ka = 0.1π with u = 1.0. For a slightly increased initial quasi-momentum or on-site interaction a small decay of about 10% on the timescale of the simulation becomes visible. Further increasing ka and/or u causes the decay rate to increase rapidly until for u 3.0 at ka = 0.1π and ka 0.3π at u = 0.1, the current has completely vanished at time tJ = 10. Figure 6 shows that the dependence of the time evolution of the condensate fraction C 100 on ka and u is qualitatively very similar to that of ĵ .
An important result seen in both figure 5 and figure 6 is that we do not find a sharp transition point for u and ka separating parameter regions where the current is stable or unstable. This lack of a sharp transition is expected for one-dimensional systems, and is in contrast to the results found in higher dimensions [33]. Figures 7 and 8 summarise the dependence of the current decay on ka and u, and amount to stability diagrams for the crossover between parameter regimes where the current is stable and unstable. In particular, figure 7 is a shaded plot of the loss of the condensate fraction ∆C R (τ ), which we define as We plot ∆C R (τ ) as a function of ka and u at a specific time τ J = 10, as may be measured directly in an experiment. We chose R = 100, but found that ∆C R (τ ) is independent of R, based on a comparison of calculations with R = 50, 100, and 200 sites. We see that stable currents exist on this timescale only for small values of u 1.5 and ka 0.15π, while no stable current is present for u 2.5 and ka 0.25π. In the intermediate region we observe a smooth crossover between those two regimes. In the classical limit of small on-site interactions, the stability/instability crossover tends to occur at large values of ka ≈ 0.5π, which corresponds to the dynamical instability of reference [19], while for ka ≈ 0 the instability sets in close to the value u c ≈ 3.37, which corresponds to the SF-MI transition at zero current [37]. As a reference, in figure 7 we plot the mean-field result of references [33,34] as a solid black line. The latter provides a reasonable indication of the position of the crossover region for small u 0.5 only. Outside of this region, the decay appears typically to be much faster than is expected from the mean-field estimates.
Note that ∆C R (τ ) discussed above can be directly observed via interference patterns in momentum distributions [1], which can be measured via time-of-flight measurements in experiments. We can estimate the corresponding experimental timescales by taking a typical lattice depth of 10 E R , where E R denotes the recoil energy of the atoms E R ≡ p 2 r /2m with the recoil momentum p r ≡ h/λ l . For these lattice depths, the tunnelling amplitude J/ for Rb atoms is of the order of 100 Hz. Thus, the timescale τ J = 10 corresponds to experimental timescales of the order of 100 ms, which is within typical coherence times. Figure 8 shows a time-independent stability diagram for the current, which we compute from decay rates Γ C for the condensate fraction C 100 . As was shown in figure 6, the decay of C 100 as a function of tJ can be separated into three regions. Initially, the system shows an initial decay behaviour on short timescales tJ < 1, followed by an intermediate region where the decay is found to be approximately linear in time. Finally a saturation of the decay process takes place while the system approaches the zero-current steady state. We therefore extract the decay rate Γ C by determining the slope of the decay in the approximately linear intermediate region. Note that in the short time region where we fit, a linear decay behaviour is also equivalent to an exponential decay C 100 (t) ∝ exp(−Γ C t). Figure 8, which is a shaded plot of Γ C , for Γ C /J ≤ 0.1, shows results which are qualitatively similar to those of figure 7. That is, stable currents are found to exist only for small values of ka and u, and the stability/instability crossover in general occurs at values of ka and u significantly smaller than predicted by mean-field theory (see solid black line in the figure).
In the next section we compute decay rates for various lattice fillingsn ≥ 1 and investigate the scaling with the interaction strength and filling factor. We find suprisingly that the scaling is very similar to those computed for a weakly interacting system withn ≫ 1 using beyond-mean-field (phase-slip) calculations.

Scaling of the current decay rates
In this section, we study the dependence of the decay rate Γ C defined in section 4.3 and of an analogous time-independent decay rate Γ j for the current (defined below) on the lattice fillingn, the interaction strength u, and the quasi-momentum ka. We find scaling laws exhibiting exponential dependence on these quantities. Interestingly, these scaling laws are similar to those computed for current decay in a weakly interacting system with large filling factor n ≫ 1, which were computed via (beyond-mean-field) instanton calculations in reference [34]. This is especially true in the casen ≥ 3, or far from the SF-MI transition withn = 1. Naturally, significant deviations occur between the exact numerical results close to the SF-MI transition and predictions for a weakly interacting system. Note that here we do not attempt to verify these instanton predictions within their regime of validity, which has been done in a recent study by Danshita and Polkovnikov [47], using particle numbersn ∼ 1000. Instead, we are interested in exploring the regimes close ton = 1, which correspond to recent experiments with atoms confined in 3D optical lattices.
In reference [34], a decay rate Γ for the current in a system with large filling factor was calculated using instanton (phase-slip) techniques. This could be shown to be Γ ∝ e −S , with S the semi-classical action Equation (6) is strictly valid in 1D in the weakly interacting (Gross-Pitaevskii) limit, u/n ≪ 1, for un ≫ 1, i.e. in the limit of large filling factors, and close to the classical instability at ka ≈ π/2. This is very different from the parameter regimes that are typical when the 1D system is formed via a 3D optical lattice, as we study here in the present article. In principle, our filling factorsn ≤ 4, and acceleration values 0.2 ≤ ka/π ≤ 0.35, do not correpond to the region of validity of equation (6). However, as we will see below, we nonetheless obtain very similar scalings.  Figures 9(a) and 9(b) show our results for the decay rates Γ C and Γ j , respectively, where Γ j has been extracted from the time evolution of the boson current ĵ (see figure 5) in the same way as for Γ C as described in section 4.3. In panels (a) and (b) we plot − ln(Γ C /J) and − ln(Γ j /J) respectively as a function of n/u, for 1 ≤n ≤ 4 and a fixed ka = π/4. Both of these plots show an approximately linear dependence on n/u, for alln and n/u 0.6. In addition, for these values of n/u both panels show that by increasingn the slope of the lines tend to converge to a fixed value. This indicates that not only is the decay rate proportional to exp(−λ 1/u), but that the constant λ appears to scale as √n forn 2. This is the same scaling as predicted for the casen ≫ 1 via instanton methods. Averaging over the slopes of the results for the large filling factorsn = 3, 3.5 and 4, in order to obtain the prefactor in λ, we find values of λ = (2.1 ± 0.5) √n and λ = (2.7 ± 0.6) √n for panels (a) and (b) respectively. We note that the equality of the slopes obtained from the current decay and condensate fraction decay within the fitting error indicates the equivalence of these quantities as good dynamical observables, for the characterisation of the current decay both in the numerical simulations and in experiments, for large enoughn. For smalln ≈ 1 we observe slightly different behaviour (see below for more details).
The results above for largen surprisingly show good agreement with scalings from the instanton calculations for a very different parameter regime, although we note that the prefactor in λ is different, as the scaling law in equation (6), would imply λ = 3.9 √n . The most obvious deviation between the iTEBD results and equation (6) takes place for the Γ C results of panel (a) in the vicinity of the SF-MI transition withn = 1, which occurs at n/u ≈ 0.54. For n/u 0.6 the rate Γ C no longer increases as n/u decreases. We indicate this parameter region with an arrow in figure 9(a). This behaviour is not captured by the decay of the current Γ j in panel (b). This contrast between the the decay of the condensate fraction and the decay of the current could reflect a reduced role of the condensate fraction in the overall system dynamics when we have an initial MI state. Note that for n/u 0.6 the decay of ĵ occurs rapidly on the timescale of our calculation time steps, which makes computing its decay rate less accurate in this regime.

Summary and Outlook
In summary, the iTEBD method allows us to make quantitative predictions for the time-dependence of currents for bosons moving in 1D in an optical lattice. In contrast to the behaviour in higher dimensions, we observe a broad crossover between regions of stable and unstable currents, with the typical decay rates for the currents being more rapid than mean-field predictions. We also find surprising agreement with the scaling of decay rates with interaction strength and filling fraction calculated for the case of large fillings and weak interactions via phase-slip methods. These results should be accessible in current experiments with filling factors near unit density. As the density increases, we find that the scaling of the decay rate of the current almost agrees with predictions from instanton calculations.
These results are strongly related to the decay of currents generated in a displaced harmonic trap, which has also been the subject of significant experimental [16,17,20] and theoretical [19,21,22] investigation. An understanding of these currents is also an important starting point for the investigation of more general transport properties of bosons in 1D. This could include behaviour of these currents in the presence of an impurity [7,8], or spatially varying interactions [9,10].

Acknowledgments
We thank D. Jaksch, S. Clark, E. Demler and P. Zoller for helpful discussions. This work was supported by the Austrian Science Foundation (FWF) through SFB F40 FOQUS and project I118 N16 (EuroQUAM DQS), the DARPA OLE program, STREP FP7-ICT-2007-C project NAME-QUAM, and the Austrian Ministry of Science BMWF as part of the UniInfrastrukturprogramm of the Forschungsplattform Scientific Computing and of the Centre for Quantum Physics at LFU Innsbruck.

Appendix A. Translationally invariant block size and the iTEBD Algorithm
The original iTEBD algorithm as introduced in [35] uses translationally invariant matrix product states to represent the quantum state of an infinite homogenous system. In order to conveniently compute time-evolution, the translationally invariant state is usually represented by a repeated block of two sites, as described in reference [35]. Note that despite the block size of two, excitations with any period can be represented in this form, because matrix product states effectively store phase relationships between neighbouring sites rather than phases corresponding to a single site. Indeed, the use of a block size of two is a matter of convenience in applying the algorithm, and if variational methods are applied [11], then a block size of one can be sufficient to represent a translationally invariant state. In the present work, we find it convenient to extend the block size beyond two to m sites, where for a particular simulation in which the momentum translation of the initial state is ka, we choose m = 2π/ka. This is convenient in the application of the operatorK(ka), but has no fundamental effects on the computation of time evolution of the state. In our case, time evolution is simulated by decomposing the time evolution operator into m twosite unitary operations, instead of two as in the original iTEBD algorithm. We expect that this, in case of a large enough bond dimension χ still results in an exact state representation for the infinite system during the whole time evolution.
To demonstrate the independence of the computation of time evolution on our choice of m, we plot example values for the time evolution of the condensate fraction C 100 , computed with increasing m in figure A1. We use an initial acceleration ka = π/4, and compare results for block sizes m = 8, 16, 24, for two different values of the interaction strength u = 1, 3. We see that the results are independent of the block size, and after testing for convergence in the parametes χ and d should represent exact time-evolved values.