Coherent matter waves emerging from Mott-insulators

We study the formation of (quasi-) coherent matter waves emerging from a Mott insulator for strongly interacting bosons on a one-dimensional lattice. It has been shown previously that a quasi-condensate emerges at momentum kcond = π/2a, where a is the lattice constant, in the limit of infinitely strong repulsion (hard-core bosons). Here, we show that this phenomenon persists for all values of the repulsive interaction that lead to a Mott insulator at a commensurate filling. The non-equilibrium dynamics of hard-core bosons is treated exactly by means of a Jordan–Wigner transformation, and the generic case is studied using a time-dependent density matrix renormalization group technique. Different methods for controlling the emerging matter wave are discussed.


Introduction
Cold atoms in optical lattices currently constitute one of the most flexible and hence most promising ways to investigate areas of many-particle physics for which no established knowledge prevails. One such area is that of strongly correlated systems, for which neither the ground-state nor excited states can be easily described and for which no generic analytic tools with predictive power are currently available. The situation is even more difficult for non-equilibrium processes, where, in general, a full knowledge of the eigenstates of the system is needed in order to properly describe its dynamics out of equilibrium.
Recent experimental work has succeeded in producing one of the hallmarks of strong correlations, namely a Mott insulator [1,2]. Remarkably, this was achieved with bosons, a situation not easily encountered in condensed matter physics. Even the limit of very strong interactions, where the bosons can be considered to be impenetrable particles (hard-core bosons or a Tonks-Girardeau gas [3]) could be reached experimentally [4,5]. Such a limit is attainable in elongated traps for large positive three-dimensional scattering lengths, at low densities, or with very strong transversal confinement [6]- [8].
On the other hand, the study of the non-equilibrium dynamics of quantum gases was instrumental for understanding their properties [9]. More recently, a controlled study of the evolution of one-dimensional Bose gases prepared initially out of equilibrium was performed [10], opening up the possibility of investigating the role of integrability in the relaxation dynamics of a many-body system [11]. In the case of hard-core bosons in one dimension, the evolution of a system with an arbitrary initial state can be described theoretically in an exact way [12]- [14] via an exact mapping onto free fermions, the Jordan-Wigner transformation [15].
We concentrate here on the out-of-equilibrium evolution of a one-dimensional Bose gas with strong interactions whose initial state is a Mott insulator. It was previously shown that, in the hard-core limit, an initial Fock state develops quasi-long-range correlations when the bosons are allowed to evolve freely on a lattice [12,14]. In particular, the momentum distribution function n k ≡ n k develops sharp peaks at momenta k = ±π/2a, where a is the lattice constant. An examination of the one-particle density matrix shows that, after the formation of the peaks in n k , it decays as 1/ √ x at long distances, as it does for hard-core bosons in equilibrium [16,17], demonstrating that, in fact, the peaks in n k signal the emergence of quasi-coherence at a finite wavevector. A detailed picture of the (quasi-) coherent part is obtained by examining the lowest 3 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT natural orbital (NO), i.e., the eigenvector of the one-particle density matrix corresponding to the largest eigenvalue. Once the peaks in n k form, the NO evolve at a constant velocity v NO = ±2at/h, where t is the nearest-neighbour hopping amplitude, without appreciable change in their form. These are the maximal group velocities on a lattice with a dispersion k = −2t cos ka. The process of formation of the quasi-condensate is also characterized by a power law. The population of the quasi-condensate increases in a universal way as ∼1. 38 √ tτ/h, as a function of the evolution time τ, independently of the initial number of particles in the Fock state. The time τ m at which the maximal occupation of the NO is reached depends linearly on the number of particles N b in the initial Fock state, and is given by The appearance of quasi-condensates at k = ±π/2a can be understood on the basis of total energy conservation. Given the dispersion relation of hard-core bosons on a lattice, since the initial Fock state has a flat momentum distribution function, its total energy is E T = 0. If all the particles were to condense into one state, it would be to the one with an energy k = E T /N. Taking into account the dispersion relation k = −2t cos ka, k = 0 corresponds to k = ±π/2a. Actually, since there is only quasi-condensation in the one-dimensional case, the argument above applies only in that the occupation of a given state is maximized. In addition, the minimum in the density of states at these quasi-momenta strengthens the quasi-condensation into a single momentum state.
Since hard-core bosons can be treated exactly [4], [12]- [15], they are extremely well-suited for a theoretical study of non-equilibrium dynamics because large systems (with hundreds to thousands of bosons) can be examined over long times. However, the experimental investigation is hampered by the quite stringent requirements for the realization of such systems. We therefore consider here the case of finite interactions, modelled by the one-dimensional Hubbard model where b † i and b i are bosonic creation and annihilation operators, respectively, and n i = b † i b i is the density operator. The hard-core limit corresponds to U → ∞. The value at which the Mott insulator appears has been estimated as U c /t ∼ 3.5 in one dimension for a commensurate density n = 1 [18]. Hence, all the cases considered here correspond to U > U c .
As in the hard-core case, we start with bosons in a Mott-insulating state spread over several lattice sites and monitor the free expansion on a lattice. For the time evolution of the system, it is, in principle, necessary to treat the whole Hilbert space of the system, restricting exact treatments to extremely small systems. Instead of fully diagonalizing the Hamiltonian matrix, efficient iterative eigensolvers such as the Lanczos or the Jacobi-Davidson procedure are commonly used [19,20], enabling one to treat somewhat larger systems.
Unfortunately, the Lanczos method is limited by the exponential growth of the Hilbert space as a function of the number of degrees of freedom. Therefore, a more efficient way of representing the relevant subspaces for the time evolution is needed. Recent progress in this direction was achieved by extending the density matrix renormalization group (DMRG) method [21,22] to treat the time evolution of correlated systems [23]- [28], leading to the so-called t-DMRG.
Here, we apply the t-DMRG to study the expansion of soft-core bosons out of a Mott insulator and compare it to the hard-core case. It will be shown that the essential features obtained with hard-core bosons are preserved, while new control possibilities are opened by tuning the strength of the interaction U. We will also show that further control of the momentum of the emerging quasi-condensates can be achieved by introducing a superlattice, which is obtained by superimposing an extra periodic potential on to an already existing lattice potential. Such systems have been recently realized experimentally with ultracold gases trapped on optical lattices ( [29,30] and references cited therein), and have been studied theoretically by various mean field approaches [31], quantum Monte Carlo simulations [32], and exact diagonalization [33].
In section 2, we discuss the theoretical treatment of the time evolution of a general quantum system both within the framework of the Lanczos method and of the t-DMRG. The results are shown in section 3, where the evolution for parameters in the range 6 U/t 40 are considered. As in the hard-core case, the momentum distribution function n k displays maxima at finite wavevectors. However, those wavevectors are displaced to lower values of k as the strength of the interaction is reduced. In section 4, we study how the introduction of a superlattice allows further control of the emerging quasi-condensates. For these systems we restrict our analysis to the hard-core regime. Finally, a concluding discussion is given in section 5.

Time evolution of many-body quantum systems
The general solution of the Schrödinger equation for a many-body system can only be given in a formal way, that for a time-independent Hamiltonian is where H is the Hamiltonian determining the evolution over a time τ from some initial time τ = 0. For sufficiently small systems, full diagonalization of H is possible, and a knowledge of all the eigenvalues and eigenstates allows for an exact determination of the evolved state. However, for a general many-body system, the size of the Hilbert space grows exponentially with the number of degrees of freedom, restricting the system size essentially to tens of atoms. More efficient ways of determining the time evolution are discussed in the following subsections.

Lanczos method
Here, we focus on the Lanczos procedure, which can be generalized in a straightforward way so that the time evolution of the system can be computed without calculating all eigenstates. The basic idea is to expand the time-evolution operator, and to focus on the set of states {|ψ(τ) , H|ψ(τ) , . . . , H m |ψ(τ) }, which spans the so-called Krylov subspace.
The key idea of the Lanczos method is to obtain a basis by orthogonalizing the vectors of the Krylov subspace only with respect to the previous two elements of the set, leading to the recursion relation

5
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT for the vectors |v i of the Lanczos basis. The projection of the Hamiltonian on to this basis set leads to a tridiagonal matrix T m = V T m HV m , where V m is a rectangular matrix containing the Lanczos vectors as column vectors. In this way, the Hamiltonian T m can be efficiently diagonalized.
Using this approach, an approximation for the time evolution of a given state |ψ(τ) at time τ over a small time interval τ can be given. For a time-independent Hamiltonian it reads Remarkably, an exact bound can be given for the error in the approximation [34]: valid for m ρ τ/2, where ρ is the width of the spectrum of H and m is the dimension of the Krylov subspace. For ρ τ 1 and m sufficiently large (of the order of 10), the formula above shows almost exponential convergence.

Adaptive time evolution with the DMRG
The basic idea of the DMRG method is to represent one or more pure states of a finite system approximately by dividing the system in two and retaining only the m most highly weighted eigenstates of the reduced density matrix of the partial system. In combination with the numerical renormalization group (NRG) approach developed by Wilson [35] and the superblock algorithms developed by White and Noack [36], this leads to a very powerful and efficient tool for the investigation of one-dimensional strongly correlated quantum systems on a lattice. Here, we only give a rough sketch of the method and refer to recent reviews [20,37,38] for a detailed description.
As depicted in figure 1, the key steps are to increase the number of degrees of freedom of the partial system by adding sites, then to decrease the number of degrees freedom by retaining states below a cutoff. In this way, the method carries out a renormalization group procedure closely related to Wilson's NRG.
In the first step of the algorithm, a site is added to one of the subsystems, with its Hamiltonian exactly represented. The Hamiltonian of the subsystem is usually represented in an efficient reduced basis built up from the m most important eigenstates of its reduced density matrix. Note that the basis is incomplete due to the truncation. A measure of the error ε introduced by the cutoff is the discarded weight, which measures the total weight of the discarded states where λ j is the jth eigenvalue of the reduced density matrix. After convergence is reached, one typically obtains values ε < 10 −6 with m < 1000. Although one is only retaining a tiny fraction of the total Hilbert space of the system (which already for small systems can reach a dimension of several million), the desired observables can thus be obtained with high accuracy.
In the second step of the iteration, the states one is interested in are obtained. These states are called 'target states'. In the original ground-state algorithm, these are the ground state and the few The lattice is shown in the usual 'superblock' configuration, where the left part of the lattice is the subsystem which is used to compute the basis of density matrix eigenstates. At the 'dividing' bond, two 'exact' sites are added; the 'sweep' proceeds from left to right. The flowchart at the right-hand side shows the relevant steps of the DMRG procedure as described in the text.
lowest lying excited states of the system, which are obtained by diagonalizing the Hamiltonian of the total system, e.g., by carrying out the Lanczos diagonalization algorithm described in the introduction. However, states other than the eigenstates of the Hamiltonian may be obtained in this step. This flexibility is crucial for the time-evolution algorithms, in which the time-evolved state is obtained by other means than solving an eigenvalue problem.
In the third step, the new effective basis is obtained by diagonalizing the reduced density matrix of the extended subsystem, given by where the sum goes over all target states. In step four, only the m eigenstates with the largest eigenvalues are kept. The operators needed to represent the Hamiltonian of the subsystem, to form the pieces of the Hamiltonian connecting subsystems, and to calculate observables are transformed into this new reduced basis. This effective Hamiltonian of the subsystem is now the starting point for step one of the next iteration. In this way, every step improves the accuracy of the obtained eigenstates and energies by improving the reduced basis used for the representation of the target states. DMRG iteration schemes are usually divided into two classes. In the so-called infinitesystem algorithm, the system grows at each step. This can be used to build up the system up to a desired lattice size. In the finite-system algorithm, the size of the lattice is fixed, and the 'dividing bond', i.e., the position at which the system is cut in two parts, is moved from the right end of the lattice to the left end and back (other variations are possible). This is called a 'sweep'. In order to obtain the ground state (and optionally the lowest lying excited states) with a high accuracy, the sweeps are iterated until convergence is reached. The calculation can be significantly speeded up if the diagonalization in step 2 of the DMRG procedure is started with a good initial guess for the wavefunction. Such an initial guess can be constructed using the so-called 'wavefunction transformation', which approximately transforms the wavefunction obtained from the previous finite-system step into the basis of the current superblock configuration. As we will see next, the wavefunction transformation also plays a key role in the adaptive t-DMRG schemes.
The main difficulty in calculating the time evolution using the DMRG is that the restricted basis determined at the beginning of the time evolution is not able, in general, to represent the state well at later times [24] because it covers a subspace of the total Hilbert space which is not appropriate to properly represent the state at the next time step. Since both the Hamiltonian and the wavefunction |ψ(τ) at time τ are represented in an incomplete basis, the result for the next time step |ψ(τ + τ) will have additional errors because the reduced basis is not an optimum representation for this state. In order to minimize these errors, it is necessary to form a density matrix whose m most important eigenvectors are 'optimal' for the representation of the state |ψ(τ) , as well as for |ψ(τ + τ) in the reduced Hilbert space. The most straightforward approach is to mix all time steps |ψ(τ i ) into the density matrix [24,28]. However, this can be extremely costly computationally. A more efficient way is to adapt the density matrix at each time step.
An approach for adaptive time evolution based on the Trotter-Suzuki [39] decomposition of the time-evolution operator was developed in [25]- [27]. The idea is to split up the timeevolution operator in local time-evolution operators U l acting only on the bond l. For lattice Hamiltonians containing only terms connecting nearest-neighbour sites, this is easily obtained using the Trotter-Suzuki decomposition, which in second order is given by Here, H even and H odd is the part of the Hamiltonian containing terms on even and odd bonds, respectively. Since each bond term H l within H even or H odd commutes, e −i τH can then be factorized into terms acting on individual bonds. As depicted in figure 1, in the DMRG procedure usually two sites are treated exactly, i.e., the entire Hilbert space of the two sites is included. The Trotter variant of the t-DMRG exploits this feature by applying U l = e −i τH l at the bond given by the two 'exact' sites. In this way, the time-evolution operator has no further approximations other than the error introduced by the Trotter decomposition. In particular, the error introduced by the cutoff is avoided. The wavefunction of the lattice is then updated by performing one complete sweep over the lattice and applying U l at the 'dividing bond'. In this way, only one wavefunction must be retained and it is possible to work with the density matrix for a pure state. The flowchart is sketched in figure 2. However, the method is restricted to systems with local or nearest-neighbour terms in the Hamiltonian. A more general basis adaption scheme aims at adapting the density matrix basis by approximating the density matrix for a time interval [40], steps using the Lanczos method described in subsection 2.1. This can be done easily because the Hamiltonian of the system is usually constructed anyway in the DMRG scheme. Within the restricted basis, the Lanczos iteration, equation (4), can then be performed, leading to the desired intermediate time steps computed using equation (6). Using this approach, we find that, if the time step is small enough, it is sufficient to retain only the target state |ψ(τ + τ) , so that one can work with a pure-state density matrix like in the Trotter approach. For larger time steps, it is important to mix at least the states |ψ(τ) and |ψ(τ + τ) into the density matrix. We find that it is sufficient to perform only one half-sweep in order to adapt the restricted basis, which is the minimum requirement for updating the basis on the complete lattice. The flowchart of this approach is sketched in figure 2. With this approach, it is, in principle, possible to treat more general Hamiltonians, as long as they can be treated accurately using the DMRG. However, due to the fact that, in general, one cannot work with a pure state density matrix, the dimension m of the restricted basis needed to obtain a certain discarded weight during the time evolution is larger compared to the Trotter approach described above, making this variant slower. The errors in both adaptive schemes for intermediate and long times are comparable. For the problem at hand, we therefore choose the Trotter variant (in second order) of the adaptive time-dependent DMRG.
In this study, we control the error during the time evolution by fixing the discarded weight and varying the number of basis states kept. By keeping a maximum of m = 800 density matrix eigenstates, we obtain discarded weights smaller than 5 × 10 −7 during the time evolution.
The initial state has zero density over a wide region in the system. This may lead to difficulties for the DMRG. In order to control this, we compare the time evolution of an initial Fock state of hard-core bosons obtained using the t-DMRG with the exact results from the Jordan-Wigner transformation. The maximum error as a function of time is plotted in figure 3. As shown in this figure, the maximum deviation from the exact results is smaller than 0.01 for all times considered.

Free expansion of soft-core bosons from a Mott insulator
We consider here the free expansion of N b interacting bosons described by the Hamiltonian (1) on a one-dimensional lattice with L sites, lattice constant a and open boundary conditions. In the cases considered here, we take N b = 20 and L = 60. Due to memory limitations, it is not possible to allow for all possible occupations of a given site. In general, the maximal number of bosons per site needed to have an accurate description of the system increases as U decreases. In all the cases treated here, a maximum of three bosons per site was sufficient. Even at the smallest interaction studied here (U/t = 6), no appreciable difference was observed when the cutoff was changed from 3 to 4 bosons per site. Since the system becomes more dilute in the course of the free expansion, the limitation in the number of bosons per site becomes even less important at later times. The time sequences shown are all limited to times shorter than the time it takes the matter wave to reach the boundary of the system. Figure 4 shows a comparison of the density ( figure 4(a)) and the momentum distribution function (figure 4(b)) for a system with U/t = 40 with hard-core bosons at four different times. At such high values of the interaction, it is expected that the particles behave as hard-core bosons. In fact, there is no noticeable difference in the density profiles at any time. However, the momentum distribution functions at τ = 0 show clear differences. While hard-core bosons are equally distributed over all momenta, the soft-core case has a maximum at k = 0, showing that even in a Mott insulator, the fluctuations of the number of particles at each side populate that state preferentially. Nevertheless, after the Mott insulator is allowed to expand, the difference between both systems becomes barely noticeable. Already at the second time shown in figure 4, where a Mott plateau still exists at the centre of the cloud ( figure 4(a)), the momentum distribution functions of the hard-and soft-core bosons are very close to each other. This is expected because the constraint of a hard-core should become less relevant when the system is diluted.  In order to see the establishment of coherence explicitly, we examine the spatial behaviour of the one-particle density matrix. We would expect a change from an exponential decay in the Mott-insulating state to a power law behaviour if (quasi-) coherence emerges. Figure 5 shows the spatial behaviour of the one-particle density matrix both at time τ = 0 ( figure 5(a) bosons are in a Mott-insulating state, and at time τ ∼ 7.5 ( figure 5(b)), when the peaks around k = ±π/2a are well established. Figure 5(a) shows the spatial dependence of the one-particle density matrix measured from the centre of the bosonic region in a semi-logarithmic plot for different values of U. As expected for a Mott insulator, an exponential decay with a correlation length that shortens as U is increased is observed. Figure 5(b) shows the decay of the one-particle density matrix on a log-log scale at a time long enough so that the peaks around k = ±π/2a are fully developed. The evaluation was made in the part of the system where the lowest NO is appreciable, i.e. in the region of the system where a well-developed quasi-condensate can be expected. Figure 6 shows the spatial dependence of the lowest NO at the same time as in figure 5(b). The correlations in figure 5(b) were measured from the site x j = 37a (a position where the NO is well developed) and with x i > x j . For comparison, we superimposed the oneparticle density matrix for U/t = 20, N b = 20, and L = 60 in equilibrium where, due to the lower density with respect to the initial state in figure 5(a), a quasi-condensate exists. Over the distances where the NO has an appreciable value, no difference with the corresponding quantity in equilibrium can be noticed. It can be clearly seen that the one-particle density matrix has developed a power-law decay (the same as the one in the system in equilibrium) at the later time, with a power that approaches the one of hard-core bosons. Unfortunately, the expansion of the cloud and the total system size are not as large in the soft-core (as seen in the departures from the power-law due to finite-size effects) as in the hard-core case, so that the exponent of the power-law decay cannot be as accurately determined. Nevertheless, it is clear that a change from an exponential to a power-law decay takes place, indicating that a quasi-coherent matter wave has developed. Finally, we discuss the behaviour of the expansion for smaller values of the interaction U. Although at U/t = 40, the momentum distribution function closely follows the shape of n k of hard-core bosons, a tiny asymmetry can be seen in figure 4(b) around the peaks at k = π/2a. Such an asymmetry indicates that the maximum of n k is not exactly at k = π/2a, but is shifted slightly. A more detailed analysis for 6 U/t 40 is presented in figure 7. Figure 7(a) shows n k around k = π/2a with the data points from t-DMRG denoted by symbols with spline interpolations between them. It is clearly seen that the maximum of n k is displaced to smaller momenta as U decreases. The spline interpolation allows for a better determination of the maxima in n k , since a denser set of k-points corresponds to having a much longer lattice in a physical realization. On the other hand, the actual set of k-points in the t-DMRG simulation corresponding to L = 60 is dense enough to allow for a smooth interpolation without introducing artefacts due to the spline procedure. Figure 7(b) shows the location of the maxima of n k as a function of U in units of 2a/π, giving a guide for a fine tuning of the wavelength of the matter wave via the interaction strength.
The results in this section show that the main feature found for the case of hard-core bosons [12,14], namely, the emergence of a quasi-coherent matter wave from a Mott insulator, persists under more general conditions. The wavelength of the matter wave is determined primarily by the underlying lattice on which the expansion takes place, as in the case of hard-core bosons. Additionally, finer tuning is possible by regulating the interaction strength of the bosons. This means that on an optical lattice, the wavelength of the corresponding laser beam would essentially determine the momentum of the matter wave and a fine tuning can be reached by varying its intensity. The next section discusses further control possibilities by introducing more complex structures in the optical lattice.

Expansion in a superlattice
So far we have studied how finite values of the on-site interactions between bosons modify the behaviour already known in the hard-core limit. The general finding has been that the physics is similar, and that an emergence of quasi-condensates can be obtained experimentally for a wide range of finite repulsive interactions. In this section, we analyse how the introduction of a superlattice potential allows a further degree of control over the system and can lead to a richer momentum distribution of the expanding cloud of bosons. We will restrict the analysis to the hard-core limit keeping in mind its relevance to the soft-core regime.
In the presence of a superlattice potential, the hard-core boson Hamiltonian becomes with the additional on-site constraints which exclude double or higher occupancy. The bosonic creation and annihilation operators at site i are denoted by b † i and b i , respectively, and the local density operator by n i = b † i b i . The brackets in equation (12) apply only to on-site anticommutation relations; for i = j, these operators commute as usual for bosons; [b i , b † j ] = 0. In equation (11), the hopping parameter is denoted by t and the last term represents the superlattice potential with strength A and sites per unit cell.
Using the Jordan-Wigner transformation [15] b one can map the HCB Hamiltonian on to that of non interacting spinless fermions, where f † i and f i are the creation and annihilation operators for spinless fermions at site i, and n f i = f † i f i is the local particle number operator. For periodic systems with N lattice sites, one needs to consider that so that when the number of particles in the system ( i n i = i n f i = N b ) is odd, the equivalent fermionic Hamiltonian satisfies periodic boundary conditions; whereas, if N b is even, antiperiodic boundary conditions are required.
The above mapping to non interacting fermions allows one to realize that the presence of an additional periodic potential opens gaps at the edges of the reduced Brillouin zones. This implies that new insulating phases appear at fractional fillings n i = i/ , with i = 1, . . . , − 1, in addition to the insulating phase at n = 1 which is present in the absence of the superlattice.
In the following, we address the question of what happens during the evolution of initially prepared insulating states when they are allowed to expand in a superlattice. For simplicity, we restrict our analysis to the case = 2. (The generalization to larger values of is straightforward.) To study these systems, we follow the exact approach already described in detail in [12]- [14].
where '+' stands for the upper band and '−' for the lower one. One can then calculate the group velocity at each momentum as which means that the wave vectors at which the maximum group velocity occurs satisfy For these values of k, the density of states attains its minimum values. Hence, under the appropriate initial conditions, we should expect the quasi-condensates to emerge at these values of k rather than at k = π/2, as in the absence of the superlattice (A = 0). Since for = 2 the only insulating phases occur at full filling and at half filling, we start studying the evolution of an insulating state initially prepared with one particle per lattice site, like in the previous sections and in [12,14]. The expansion of such a state with 101 HCB's is shown in figure 8. While the evolution of the density ( figure 8(a)) is very similar to that in the absence of the superlattice [12,14], the evolution of the momentum distribution function is completely different ( figure 8(b)). No peaks appear in n k , in contrast to the case analysed in [12,14] and in the previous sections. This is because in the superlattice the mean energy per particle ( = 0) for the fully filled insulating state lies in the band gap, and the states with the closest energy are the ones with momentum k = ±π/2, which have ν g ± = 0 and the maximum density of states, i.e., exactly the opposite of the case without the superlattice.
A scenario in the superlattice that is closer to the one in a Mott insulator with n = 1 (and finite U) in the absence of the superlattice, is the one in which the initial state is prepared with a mean density of 0.5. Such state has short-range (exponentially decaying) one-particle correlations like the ones seen in figure 5. In addition, its mean energy per particle lies within the lowest band, where the density of states is finite.
In figure 9, we show the evolution of the mean density per unit cell and momentum profiles during the expansion of a state prepared in a half-filled box with 100 HCB's and A = 2t. As can be seen in figure 9(b), the initial momentum distribution is not flat like the one in figure 8(b). Its maximum at k = 0 signals, the presence of short-range correlations like in the Mott insulator of figure 4. Remarkably, during the expansion of this state, sharp peaks appear in n k at ka = ±0.87 and ka = ±2.27, which are the momenta for which the group velocity is maximum and the density of states has a minimum, following equation (18).
The peaks in n k signal the emergence of quasi-condensates of HCB's at finite momentum. This can be seen by studying the NOs (φ η ), which can be considered to be effective single-particle states in interacting systems, and are defined as the eigenfunctions of the one-particle density matrix ρ ij [41], with occupations λ η . In dilute higher dimensional gases, where only the lowest NO (the highest occupied one) scales ∼N b , the occupation of this orbital can be regarded as the BEC order parameter, i.e., the condensate [42].  Figure 10. Evolution of (a) the NO occupations and (b) the modulus of the lowest NO wavefunction (averaged in the unit cell) for 100 HCB's initially prepared in an insulating state with mean density of 0.5 on 1000 lattice sites. The superlattice parameters are = 2 and A = 2t. The inset in (a) shows the initial exponential decay of one-particle correlations (averaged per unit cell), and its conversion to a power-law decay in the region where the quasi-condensate forms. The black line corresponds to 1/ |x i − x j |. The times are τ = 0 ( ), 50h/t ( ), 200h/t (•), and 400h/t ( ).
In figure 10(a), we show the values of λ η for the 200 highest occupied orbitals after the same expansion times as in figure 9. One can see that during the expansion the lowest NO becomes 'highly' populated with of the order of √ N b particles. This is because quasi-long range correlations develop in the system, as shown in the inset of figure 10(a). The power-law decay of the one-particle correlations is ∼1/ |x i − x j |, i.e., the same form that appears in the absence of the superlattice [12,14], and that has been proven to be universal in the ground state [16,17]. The wavefunction of the lowest NO during the expansion is depicted in figure 10(b). One can see that it exhibits exactly the same features observed in [12,14] when there was no additional periodic potential. After the Mott insulator melts, the shape of the lobes of the NO stops changing and they just move with the maximum velocity in the lattice given by equation (17) for k = k m .
One final remark on these emerging quasi-condensates is in order. In contrast to those emerging in a system without a superlattice potential which are mainly formed by particles with ka = ±π/2, we find that in the superlattice the quasi-condensates are mainly formed by HCB's with the four momenta k m . This is reflected by the four-peak structure of the Fourier transform of φ η at momenta k m , depicted in figure 11. Hence, the four peaks that appear in the momentum distribution in figure 9(b) reflect the formation of the quasi-condensates in figure 10(b), which have a richer structure in k-space than those that emerge in the absence of the superlattice.

Conclusions
In this study, we have presented non-trivial generalizations of a non equilibrium phenomenon originally found in systems of hard-core bosons that expand out of strongly correlated Mott-insulating states. First, we have shown that the emergence of coherent matter waves at finite momenta persists away from the hard-core limit, i.e., to finite values of the on-site repulsion U which range down to the critical values, as long as a Mott-insulating state is attainable for the initial state. The accurate treatment of this complicated many-body problem has been made possible by the advent of the t-DMRG. By comparing to exact results in the hard-core case, we have shown that our t-DMRG calculations, although approximate, are very reliable. Two new features appear in the soft-core case. (i) Although the time evolution of the density profiles is indistinguishable from those of hard-core bosons at large values of U, the initial momentum distribution functions are markedly different. Nevertheless, after the Mott region melts, the momentum distribution functions of both systems become nearly identical. (ii) As the strength of the interaction is reduced, a shift of the momentum of the coherent matter wave to values smaller than π/2a is observed. The appearance of a power law in the spatial decay of the one-particle density matrix demonstrates explicitly the (quasi-) coherent nature of the resulting matter wave. Furthermore, we have shown that a coherent matter wave can be also obtained in the presence of superlattice potentials given an adequate selection of the initial insulating state. In the latter case, the emerging quasi-condensates have a more complicated structure in momentum space, with a number of dominating momenta whose locations depend on the superlattice used.
The results of this study show that it is possible to engineer atom lasers with a high degree of control. The momentum of the coherent matter wave can first be regulated by setting the wavelength of the underlying optical lattice and further fine-tuned by regulating the depth of the potentials in the lattice, i.e., the intensity of the corresponding laser beam. A further finite shift of the momentum can be achieved by superimposing another laser beam with a commensurate