Time evolution of Matrix Product States

In this work we develop several new simulation algorithms for 1D many-body quantum mechanical systems combining the Matrix Product State variational ansatz with Taylor, Pade and Arnoldi approximations to the evolution operator. By comparing all methods with previous techniques based on Trotter decompositions we demonstrate that the Arnoldi method is the best one, reaching extremely good accuracy with moderate resources. Finally we apply this algorithm to studying the formation of molecules in an optical lattices when crossing a Feschbach resonance with a cloud of two-species hard-core bosons.

In this work we develop several new simulation algorithms for 1D many-body quantum mechanical systems combining the Matrix Product State variational ansatz with Taylor, Padé and Arnoldi approximations to the evolution operator. By comparing with previous techniques based on Trotter decompositions we demonstrate that the Arnoldi method is the best one, reaching extremely good accuracy with moderate resources. Finally we apply this algorithm to studying how correlations are transferred from the atomic to the molecular cloud when crossing a Feschbach resonance with two-species hard-core bosons in a 1D optical lattice. The Density Matrix Renormalization Group (DMRG) method is a successful technique for simulating large lowdimensional quantum mechanical systems [1]. Developed for computing ground states of 1D Hamiltonians, it is equivalent to a variational ansatz known as Matrix Product States (MPS) [2]. This relation has been recently exploited to develop a much wider family of algorithms for simulating quantum systems, including time evolution [3,4], finite temperature [4,5] and excitation spectra [6]. Some of these algorithms have been translated back to the DMRG language [7] using optimizations developed in that field and introducing other techniques such as Runge-Kutta or Lanczos approximations of the evolution operator [8,9].
The MPS are a hierarchy of variational spaces, S D , [See Eq. (1)] sorted by the size of its matrices, D. MPS can efficiently represent many-body states of 1D systems, even when the Hilbert space is so big that the coefficients of a pure state on an arbitrary basis cannot be stored in any computer. While the accuracy of this representation has been proven for ground states [10], evolution of an arbitrary state changes the entanglement among its parties, and a MPS description with moderate resources (small D) might stop to be feasible.
We will take a pragmatic approach. First of all, MPS and DMRG algorithms can compute truncation errors so that the accuracy of simulations remains controlled. Second, we are interested in simulating physically small problems, such as the dynamics of atoms and molecules in optical lattices. For such problems small D are sufficient to get a qualitative and even quantitatively good description of the observables in our systems. As we will see below, the biggest problem is the accumulation of truncation errors and not the potential accuracy of a given space S D for representing our states.
The outline of the paper is as follows. We first present many recently developed simulation algorithms under a common formalism based on the optimal truncation operator. We then introduce two new algorithms: one of them is based on Padé expansions of the evolution operator while the other one uses linear combinations of MPS to increase the accuracy (similar to an "Arnoldi" method). A comparison of all algorithms using spin− 1 2 models shows that all methods are strongly limited by truncation and rounding errors and that the Arnoldi method performs best for a given accuracy.
In the last part of the paper we apply these algorithms to study a model of hard-core bosonic atoms going through a Feschbach resonance. Current experiments [11] with such systems have focused on the number and stability of the formed molecules. In this work we focus on the 1D many-body states and show that coherence is transferred from the atomic component to the molecular one, so that this procedure can be used to probe higher order correlations in the atomic cloud.
Summary of algorithms.-The space of MPS of size D is the set of states of the form Here i k = 1 . . . d labels the physical state of the k-th lattice site, the Hilbert space has a total of d N states and we sum over repeated indices. Since S D is not a vector space, we cannot solve a Schrödinger equation directly on it. Our goal is rather to approximate the evolution at short times by a formula of the kind |ψ(t + ∆t) ≃ P D [U n (∆t)|ψ(t) ] . The operator P D optimally projects an arbitrary vector onto the manifold S D , and U n (∆t) = exp(−iH∆t) + O(∆t n ) is itself an approximation to the evolution operator.
This formulation applies qualitatively to all recently developed MPS and DMRG algorithms. For instance, Refs. [3,4,7] use a Trotter decomposition of a nearest neighbor Hamiltonian H = m H m,m+1 while the Runge-Kutta algorithm in Ref. [8] is equivalent to a fourth order expansion of the exponential In practice there are differences between all cited implementations, such as truncating in between applications of the exponentials or replacing the previous Taylor expansion by a product of truncated first order monomials, k P D (1 + α k H). In this work we present two new methods. The first method uses a Padé approximation to the exponential computed with same order polynomials in the numerator and denominator [12]. The second method combines ideas from Refs. [6,9]. First of all, it is possible to prove that a linear combination of MPS, such as Nv k=1 c k |φ with initial condition |φ 0 := |ψ(0) . Finally, defining the matrices N ik := φ j |φ i and H ik := φ j |H|φ i we compute an Arnoldi estimate of the exponential This algorithm involves several types of errors. The error due to using only N v basis vectors is proportional the norm of the vector φ Nv as in Lanczos algorithms [13]. Truncation errors arising from P D can be computed and controlled during the numerical simulation. It is important to remark that while the errors in Eq. (4) may be compensated by adding more vectors, the biggest errors arise from the final truncation in Eq. (5).
Implementation.-All algorithms build on the same nonlinear operator, P D , which optimally approximates a state using the elements of S D . To truncate a linear combination of vectors we use the definition [4] If we rather want to approximate the action of an operator that can be decomposed as U = X −1 Y , we will apply a generalization of the correction vector method [1] This last formula has been used for the Padé (3) equation and for the Runge-Kutta methods, where X and Y are polynomials of the Hamiltonian. One may quickly devise a procedure for computing the minimum of Eq. (7) based on the definition of distance: (8) All scalar products in Eq. (8) are quadratic forms with respect to each of the matrices in the states |φ and |ψ . The distance is minimized by optimizing these quadratic forms site by site, or two sites at a time, sweeping over all lattice sites until convergence to a small value which will be the truncation error [4].
The Comparison.-We tested all algorithms by simulating the evolution under the family of spin-1 2 Hamiltonians with nearest-neighbor interactions of an initial state, |ψ(0) ∝ (|0 + |1 ) ⊗L , where |0 and |1 are the eigenstates of s z . By restricting ourselves to "small" problems (L ≤ 16), we could compare all algorithms with accurate solutions based on exact diagonalizations and the Lanczos algorithm [13].
In Fig. 1a we plot the errors of the different methods, ε := ψ D (T ) − U (T )ψ(0) 2 , for 8 spins, D = 16, J = 2B and ∆ = 0. Since S D contains all possible states, P D = 1 and we expect no truncation errors. Therefore, for medium to long time steps the Trotter, Taylor We have compared the best methods, that is the Trotter decomposition and the Arnoldi methods with 4 and 8 vectors, using a bigger lattice with 16 spins, and smaller matrices, D = 64, 80 and 128. As Fig. 1b shows, the XY model establishes entanglement along the lattice, the smaller matrices cannot account for this and the errors grow exponentially. Figure 2a shows how the errors in the Arnoldi method are correlated to the errors made when approximating the exact solution with a MPS of fixed size. One could think that by increasing the number of vectors we should be able to also increase the integration time-step and decrease the truncations even for fixed D, but as Fig. 2b shows, this not entirely true.
Summing up, one should use the method that allows for the longest time steps and the least number of truncations (or applications of P D ) and mathematical operations. All methods have an optimal time step which is a compromise between the errors in U n and the rounding and truncation errors made on each step. Regarding performance, the Arnoldi decomposition is competitive with the Trotter decomposition as it reaches the same accuracy for longer time steps. Furthermore, there is a huge potential for parallelizability, roughly O(N v N H /L), which all other presented algorithms lack.
The fact that previous results are model-independent is confirmed by a systematic scanning of all possible Hamiltonians in Eq. (9). A selection is shown in Fig. 1c-d. The Arnoldi method is shown to be the most accurate one, even for gapless problems. When the Arnoldi method fails it is due to truncation errors -the evolved state cannot be accurately represented by MPS-, which affect equally all MPS/DMRG algorithms [ Fig. 1d].
Molecules in the lattice.-We have used the Trotter and Arnoldi methods to simulate the conversion of bosonic atoms into molecules, when confined in an optical lattices and moving through a Feschbach resonance [11]. The goal is to study how correlation properties are transferred from the atoms into the molecules and how this dynamics is affected by atom motion and conversion efficiency.
The effective model combines the soft-core Bose-Hubbard model used to describe the Tonks gas experiments [14], with a coupling to a molecular state [15] Here, a i↑ , a i↓ and b i are bosonic operators for atoms in two internal states and the molecule; n (a) and n (m) i are the total number of atoms and of molecules on each site, and we have the usual two-level coupling with Rabi frequency Ω and detuning ∆ := E m −U m . For simplicity, we will assume that atoms and molecules interact strongly among themselves (U ↑,↑ , U ↓,↓ , U m → ∞), so that we can treat them as hard-core, a 2 i,σ , b 2 i , a iσ b i = 0. Also since molecules are heavier, we have neglected their tunneling amplitude, although that could be easily included.
As the energy of the molecular state is shifted from E m ≫ U m down to E m ≪ U m , the ground state of Eq. (10) changes from a pair of coupled of Tonks gases, to a purely molecular insulator. We want to study the dynamics of this crossover as E m is ramped slowly from one phase to the other.
The simplest situation corresponds to no hopping: isolated atoms experience no dynamics, while sites with two atoms may produce a molecule. The molecular and atomic correlations at the end of the process are directly related to two-body correlations in the initial state [16], a † k a k t=T ∼ a † k a k t=0 − n k↑ n k↓ t=0 . Therefore, this process can thus be used as a tool to probe quantum correlations between atoms. Studying the twolevel system {a † k↑ a † k↓ |0 , b † k |0 }, we conclude that for this process to work with a 90% efficiency, the ramping time should be larger than T ∼ 1.5/Ω [See Fig. 3(b)].
We have simulated numerically the ramping of small lattices, L = 10 to 32 sites, with an initial number of atoms N ↑,↓ = L/2, 3L/4. The value of the molecular coupling has been fixed to Ω = 1 and the interaction has been ramped according to E m = U ↑↓ + 4Ω(1 − 2t/T ) using the ideal ramp time T = 1.5/Ω. We have used two particular values of the inter-species interaction, U ↑,↓ /Ω = 0, 2, and scanned different values of the hopping J/Ω ∈ [0, 0.4]. The initial condition was always the ground state of the model with these values of U ↑↓ /J and no coupling. These states contain the correlations that we want to measure.
The main conclusions is that indeed the correlations of the molecules are almost those of the initial state of the atoms (11), even for J = 0.4Ω when the process has not been adiabatic. An intuitive explanation is that hopping is strongly suppressed as we approach the resonance, due to the mixing between atomic states, which can hop, and molecular states, which are slower. We can say that the molecules thus pin the atoms and measure them. This explanation is supported by a perturbation analysis at J ≪ Ω, where one finds that a small molecular contamination slows the atoms on the lattice. This analysis breaks down, however, for J ∼ Ω, the regime in which the numerical simulations are required.
Conclusions.-We have performed a rather exhaustive comparison of different methods for simulating the evolution of big, one-dimensional quantum systems [3,4,7,8] with two other methods developed in this work. All methods are substantially affected by truncation and rounding errors, and to fight the latter we must choose large integration time-steps. However, the only algorithm which succeeds for large time-steps is an Arnoldi method developed in this work. Finally, this algorithm can be applied to problems with long range interactions.
Using this algorithm, we have simulated the dynamics of cold atoms in a 1D optical lattice when crossing a Feschbach resonance. The main conclusion is that with rather fast ramp times it is possible to map the correlations of the atomic cloud (two Tonks gases in this case) and use this as a measuring tool in current experiments. This result connects with similar theoretical predictions for fermions in Ref. [16]. Simple generalizations of this work will allow us in the future to analyze losses and creation of entanglement with the help of the molecular component.