Time evolution of Matrix Product States

In this work, we develop several new simulation algorithms for one-dimensional (1D) many-body quantum mechanical systems combining the Matrix Product State variational ansatz (vMPS) with Taylor, Padé and Arnoldi approximations to the evolution operator. By comparing with previous techniques based on vMPS and Density Matrix Renormalization Group (DMRG), 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.


Introduction
The Density Matrix Renormalization Group (DMRG) method is a successful technique for simulating large low-dimensional quantum mechanical systems [1]. Developed for computing ground states of one-dimensional (1D) Hamiltonians, it is equivalent to a variational ansatz known as Matrix Product States (MPS) [2,3]. This relation has been recently exploited to develop a much wider family of algorithms for simulating quantum systems, including time evolution [4,5], finite temperature [5,6] and excitation spectra [7]. Some of these algorithms have been translated back to the DMRG language [8]- [10] using optimizations developed in that field and introducing other techniques such as Runge-Kutta or Lanczos approximations of the evolution operator [11]- [14].
The MPS are a hierarchy of variational spaces, S D , (see equation (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 [15], 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. Throughout this paper, we will distinguish between the underlying mathematical structure of MPS and the variational algorithms (Matrix Product State variational ansatz (vMPS)) developed in [3,5] and in this new study.
We will take a pragmatic approach. First of all, most algorithms in this work 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 always the potential accuracy of a given MPS space for representing our states.
The outline of the paper is as follows. In section 2, we briefly introduce MPS and review some of their properties. In subsection 2.2, we present the optimal projection onto a MPS space, which is the keystone of most evolution algorithms. In section 3, we introduce for completeness the DMRG algorithm, focusing on the concepts which are essential for time evolution. In particular, we concentrate on the difficulties faced when implementing DMRG simulations and how those techniques relate to vMPS. In section 4, we review almost all recently developed simulation algorithms under the common formalism based on the optimal truncation operator. Additionally, we introduce three new methods: two of them are based on Taylor and Padé expansions of the evolution operator while the other one uses an 'Arnoldi' basis of MPS to increase the accuracy. Section 5 is a detailed comparison of vMPS and DMRG algorithms using spin− 1 2 models. Our study shows that all methods are strongly limited by truncation and rounding errors. However, among all techniques, vMPS methods and in particular our Arnoldi method perform best for fixed resources, measured by the size of matrices D or size of basis in DMRG. In the last part of our paper, section 6, we present a real-world application of the Arnoldi evolution algorithm, which is to study a model of hard-core bosonic atoms going through a Feschbach resonance. Current experiments [16]- [18] with such systems have focused on the number and stability of the formed molecules. In this study, 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. Finally, in section 7, we summarize our results and open lines of research.

MPS
In this first section, we introduce the notion of MPS, together with some useful properties and an important operator, the projection onto a MPS space. This section is a brief review of the concepts found in [3,5].

Variational space
The techniques in this work are designed for the study of 1D or quasi-1D quantum mechanical lattice models. If N is the number of lattice components and d the number of degrees of freedom of each site, the Hilbert space of states will have a tensor product structure, H = H ⊗N 1 , where d = dim H 1 are the degrees of freedom of a single site. We will consider two examples here: one with spin-1 2 particles, where d = 2 (section 5), and later on a study of atoms and molecules in an optical lattice where d = 25 (section 6).
Given those prerequisites, the space of MPS of size D is defined as the set of states of the form Throughout the paper, we will use different notation for the matrices {A s k k }. An index k or l will always label the site they belong to and, whenever the expression is not ambiguous, the site index will be dropped: A s k := A s k k . The MPS components can also be regarded as tensors, A s k αβ , the Greek indices denoting virtual dimensions α, β = 1 . . . D. Finally, at some point, we will rearrange all values forming a vector A t k := (A 1 11 , A 2 11 , . . . , A d DD ) in a complex space C d×D×D . The first important property of the MPS is that they do not form a vector space. In general, a linear combination of M states with size D requires matrices of size MD to be represented. It is easy to do a constructive proof of the previous fact, but the reader may convince himself with a simple example, made of the two product states |0 ⊗N and |1 ⊗N , which live in S 1 , and the Greenberger-Horne-Zeilinger (GHZ) state |0 ⊗N + |1 ⊗N ∈ S 2 .
The previous remark leads us to another property, which is that the dimension D required to represent a state faithfully 1 is related to the amount of entanglement in the system. More precisely, that dimension is equal to the maximum Schmidt number of the state with respect to any bipartition and thus an entanglement monotone [19]. Indeed, creating entanglement forces us to use bigger and bigger MPS and this is the reason why some problems cannot be simulated efficiently using MPS.
The third important property is that we can efficiently compute scalar products, distances and in general expectation values of the form ψ|O 1 are local operators acting on the individual qubits or components of our tensor product Hilbert space. For instance, given that we know the matrices of those states, {A s k k } for ψ and {B s k k } for φ, the previous expectation value is made of a product of N matrices, where the 'transfer matrices' are defined as follows: Since all usual Hamiltonians and correlators can be decomposed as sum of products of local terms, the previous formulae are very useful. An important remark is that when computing equation (3) one should not directly multiply the matrices E, but cleverly contract the A and B tensors separately, so as to achieve a performance O(dD 3 ).
The last property is that expectation values, distances ψ − φ 2 , fidelities | ψ|φ | and norms ψ 2 , are quadratic forms with respect to each of the matrices in the MPS. Regarding the matrices of the state as elements of a complex vector, we can rewrite equation (2) for the kth site as where the quadratic form Q is built as follows:  This formula is used on all algorithms, from the computation of ground states [3] to the time evolution [5]. An important optimization is to avoid computing the full matrix Q, but to use the structure in (5) together with the sparseness of the transfer matrices.

Projector operator
Even if the MPS do not form a vector space, they are embedded in a bigger Hilbert space and provided with a distance. It is therefore feasible, given an arbitrary state vector, to ask for the best approximation in terms of MPS of a fixed size. The optimal projection onto S D is a highly nonlinear operation and it will be denoted in this study by P D . Following [5] 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 formula is simple to apply and behaves well for a singular operator X. Its use will become evident when studying time evolution algorithms in section 4. One may quickly devise a procedure for computing the minimum of equation (7) based on the definition of distance: All scalar products in equation (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, 2 sweeping over all lattice sites until convergence to a small value which will be the truncation error [5].
In short the algorithm looks as follows: (i) compute some initial guess for the matrices of the optimized state ψ. (ii) Focus on the site k = 1. (iii) Use equations (4) and (5) to find out the quadratic form associated with equation (8) for the first matrix (iv) The stationary points of the error are given by the equation Solve this equation and use the outcome as the new value of A k . (v) Estimate the error using equation (9). If is small enough or does not improve significantly, stop. Otherwise move to another site, k = k ± 1, and continue on step (iii).

MPS and DMRG
Even though the numerical techniques for dealing with MPS variationally seem very different from those of DMRG [1], both methods are intimately connected. First of all, the DMRG produces MPS at each state of its algorithms. Secondly, particular algorithms such as the search for ground states in open boundary condition problems are equivalent in DMRG and vMPS [3]. Thirdly, other concepts, such as basis adaptation and state transformation from DMRG, are analogous to the vMPS ones, even though they are less powerful and less accurate. In this section, we will elaborate on these statements.

DMRG builds MPS
The DMRG algorithms are based on the key idea that interesting states can be expressed using a basis with a small number of vectors. Take for instance the chain in figure 1(a). Any state of this chain can be decomposed in the form where we sum over repeated indices. While for describing individual spins we use all possible states, s k , s k+1 = 1 . . . d, in DMRG the states of the left and right blocks are expressed in a finite basis, α k , β k = 1 . . . M, where M, the number of states kept, is the basic control parameter. DMRG algorithms build those basis states recursively, by taking a smaller block, adding a site and truncating the basis of the bigger block 'optimally' in a way to be specified later. Thus we have the relations It is trivial to see, substituting those equations into equation (11), that all DMRG states are MPS [2,3].

7
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT

Targeting
An important question in DMRG is how many states we have to keep for the left and right blocks and how to optimize them. The criterion is that the basis describing those blocks has to represent accurately a family of target states, |φ n . The algorithm consists of a series of sweeps over the lattice [11,13,14,20] with a recipe to achieve an optimal representation. For instance, let us say that we have an approximate basis around sites k and k + 1 and we will improve this sweeping from left to right, to k + 1 and k + 2 (figure 1). The first step is to build the weighted density matrix of the left piece of the chain with some normalized weights, n w n = 1. Secondly, take the M most significant eigenvectors of this matrix These vectors become the improved new basis for the enlarged left block and are related by a transformation matrix to the basis elements of the smaller block (12a). Matrix elements of observables and of the initial and target states have to be recomputed using this isometry and one continues until the end of the lattice. A similar procedure is employed to create the vectors |β k recursively from right-to-left. Multiple sweeps can be performed this way.
For static problems one uses as target states, |φ n , the ground state and a number of excitations, computed with the Hamiltonian in the truncated basis. In some time evolution algorithms [11]- [14], [20] the targeting is done with respect to approximations of the time evolved state |ψ(t) . In this study, we have used |φ k := |ψ(k t/(N v − 1)) , where the intermediate and final states, |ψ(τ) , are computed using the Hamiltonian on the truncated basis (see subsection 4.2).

Targeting versus projection
There are a number of complicated subtleties and facts that are rarely mentioned in the DMRG literature. The first one is that in most implementations of the targeting algorithm, the state itself is updated and rewritten in the new basis after each step However, this leads to wrong answers because the initial state ψ(0) deteriorates during the algorithm, and it only works when the basis of the left and right block are already large enough. A more accurate procedure consists in keeping the initial state in the original basis, and on each step of the algorithm compute its matrix elements in the new basis. As an example, in figure 2, we plot the error of an evolution algorithm with the trivial algorithm and with the correct one, for an initial product state |1 ⊗L in the σ z i basis and using a simple Hamiltonian H = i σ x i . In both cases, we begin with a basis of size M = 1 vectors per block, letting the basis grow up to 2 L/2 . However only the second method is capable of doing the update.
The algorithm used either a single basis (dashed) or two sets of basis states (solid, dotted) for targeting. The dashed line shows errors due to using a single basis; the dotted line fails because the Hamiltonian does not have nonzero elements on the initial DMRG basis.
An even more important issue is that in current literature this method is formulated in an intuitive way and little is known about why it works and how accurate it is. In this respect, a close look at the matrices (12a) and (12b) created by the targeting with multiple vectors, |φ t shows that this procedure is quite similar to the computing projection operator for some linear combination P M ( k c k |φ k ), but the steps for computing the optimal approximation are wrong if there are more than one target state.
Another consequence of being tied to the notion of 'optimal basis' is that the targeting procedure does not work when the operators O n and the target states |φ n have no elements in the initial basis. A trivial example consists of the same product state |ψ(0) = |1 ⊗N as before, beginning with M = 1 states per block and evolved with H = i σ x i . As shown in figure 2, a conventional DMRG sweep cannot enlarge the basis in the appropriate way and the method fails to follow the evolution of the state. A vMPS algorithm, on the other hand, even if it begins with a state with D = 1, grows the state up to the maximal size D = 2 and makes absolutely no error in the simulation.
The final practical difference is that in current DMRG works [11]- [14] the powers of the Hamiltonian acting on the original state, H n |φ , are approximated by the truncation of H to the current basis. This is another potential source of errors which we cleverly avoid in our Arnoldi and Runge-Kutta methods shown below (subsection 4.3 and 4.5) and its effects are evident in some of the simulations (section 5).
The previous paragraphs can be rephrased as follows: DMRG algorithms need not only an initial input state, but also a suitable and large enough basis for doing the time evolution.

9
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT How to construct this basis is largely heuristic, and contrasts with the systematic way in which vMPS work.

Time evolution
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 like Here, P D is the optimal projection operator defined before and U n ( is itself an approximation of nth order to the evolution operator. Even though this formulation applies qualitatively to all recently developed vMPS and DMRG algorithms, there are subtle differences that make some methods more accurate than others. The actual implementation of time evolution is thus the topic of the following subsection.

Trotter decomposition
While there had been previous works on simulating time evolution using DMRG algorithms [21]- [23], an important breakthrough happened with the techniques in [4], which was later translated to the DMRG language [8]- [10] and has since been applied to a number of interesting problems [24]- [27]. Vidal's seminal paper [4] suggested doing the time evolution with a mixture of two-qubit quantum gates and truncation operations. The idea is to split a Hamiltonian with nearestneighbour interactions into terms acting on even and odd lattice edges Here, h i operates on the sites i and i + 1, and in our case it will correspond to a spin-spin interaction, h i = cos(θ)(s x i s x i+1 + s y i s y i+1 + s z i s z i+1 ) + sin(θ)s z i /2. Finally, H E and H O group the interaction terms acting on even and odd pairs, h 2k and h 2k+1 respectively.
The previous decomposition leads naturally to a Suzuki-Trotter decomposition of the time evolution operator, which is now approximated by a product of unitaries acting on even and odd lattice bonds: Inserting optimal truncations in between the applications of these unitaries, U j := exp(−ih j t), Vidal's algorithm can be recasted in the form of equation (16) |ψ(t + t) := Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Applying each of the unitaries U j is a relatively costless task, and in particular, for open boundary condition problems, the combination P D U j can be done in a couple of steps which amount to contracting neighbouring matrices and performing a singular value decomposition [4,8,9]. In [5], we developed a variant of this method which does not insert so many truncation operators, but waits until the end This procedure has a similar computational cost, O(ND 3 ), but the solution is then expected to be optimal for a given Trotter decomposition and can be generalized to problems with periodic boundary conditions. The accuracy of either method may be increased by an order of magnitude using a different second order decomposition but it is better to apply a Forest-Ruth formula [8,28] with the constant θ = 1/(2 − 2 1/3 ) and which has a Trotter error of order O( t 5 ). Note that better Suzuki-Trotter formulae can be designed but they do not provide a big improvement in the number of exponentials or accuracy [28].

Runge-Kutta and Lanczos with DMRG
After the development of the Trotter methods, in [11] we find a new algorithm in the field of DMRG. The idea now is to use not a Trotter formula, but a Runge-Kutta iteration: These three vectors, together with |ψ(0) , are used to find an optimal basis using the targeting procedure explained in subsection 3.2. Once the basis is fixed, the time evolution is performed with the same formulae but a smaller time step, t/10. There are variants of this technique which approximate the time evolved state using a Lanczos method [29,30] with the truncated matrix of the Hamiltonian in the DMRG basis. The different submethods differ on whether the preparation of the basis is done only using the final and initial states [14] or multiple intermediate time steps [12,13].

Runge-Kutta-like method with MPS
Our implementation of a fourth order method with Runge-Kutta-like formulae uses several simplifications. First of all, since our Hamiltonian is constant, a Runge-Kutta expansion becomes equivalent to a fourth order Taylor expansion of the exponential Here, we have rewritten the fourth order polynomial in terms of its roots, z 1 and z 2 . Using the fact that we know an efficient algorithm to compute P D Y i acting on a MPS, we can write which is our vMPS Runge-Kutta-like algorithm. There are multiple reasons to proceed this way. On the one hand, we do not want to approximate higher powers of the Hamiltonian using a truncated basis (see subsection 3.2). On the other hand, if we expand the exponential to all powers, there will be too many operators and the complexity of the state will increase enormously.
The previous decomposition has proven to be a good compromise.

Padé approximations
Runge-Kutta and in general polynomial approximations to the exponential are not unitary. Indeed, if we look at the eigenvalues of the evolution operator in either equation (24c) or (25), we will see that they are of the form λ n = 4 n=0 (iE N t) n /n!, where E n are the eigenvalues of the Hamiltonian H. From this equation, we see that |λ n | = 1 and some eigenmodes may grow exponentially, which is another way to say that Runge-Kutta algorithms are numerically unstable.
There exist multiple implicit methods that eliminate the lack of unitarity and produce stable approximations. They receive the name 'implicit' because the value of the state at a later time step t is obtained by solving an equation or inverting an operator. We will focus on Padé approximations to the exponential

12
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT which are computed with same order polynomials in the numerator and denominator [31]. The lowest order method is known as the Crank-Nicholson scheme, it arises from a second order discretization of the Schrödinger equation and has the well known form It is easy to verify that the eigenvalues of this operator are just phases and that the total energy is a conserved quantity. The other method that we have used and which we compare in this study is a fourth order expansion Applying either U CN2 or U CN4 on a matrix product state is equivalent to solving the problem in equation (7), where the operators X and Y are the denominator and the numerator in the previous quotients.

Arnoldi method
The last and most important method that we present in this paper combines many of the ideas explained before. First of all, we will use the fact that a linear combination of MPS, such as N v k=1 c k |φ (k) D , resides in a space of bigger matrix product states, S N v D , and it is thus a more accurate representation of the evolved state than a single vector of size D. In the language of DMRG, N v vectors each of size D provide us with an effective basis of size N v D. This optimistic estimate is only possible when the vectors are indeed linearly independent. Our choice for an optimal decomposition will be therefore a Gramm-Schmidt orthogonalization of the Krylov subspace, {ψ, Hψ, H 2 ψ, . . .}, performed using MPS with initial condition |φ 0 := |ψ(0) . 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 all of which can be controlled. First, the error due to using only N v basis vectors is proportional to the norm of the vector φ N v as in ordinary Lanczos algorithms [29,30]. Truncation errors arising from P D can be also computed during the numerical simulation. Out of these errors, in our experience, the final truncation in equation (31) is the most critical one, since the other ones may be compensated by adding more and more vectors.
Finally, for completeness, in this study we have also implemented a Lanczos method. It differs from the previous one in that the basis is built orthogonalizing only with respect to the two previous vectors, so that equation (30) contains only three summands, as in ordinary Lanczos iterations. One then assumes that, due to the Hermiticity of the Hamiltonian, orthogonality to the rest of the Krylov basis is preserved. Furthermore, if this is true, the effective matrices for H and N are tridiagonal and can be constructed with a simple recurrence. This method has a potential gain of O(1/N v ) in speed due to the simplifications in equation (30). However, as we will see later, truncation errors spoil the orthogonality of the Lanczos vectors and makes the method useless for small matrices.

Comparison
We tested all algorithms by simulating the evolution of the same state under a family of spin- 1 2 Hamiltonians with nearest-neighbour interactions As initial state we take the product |ψ(0) ∝ (|0 + |1 ) ⊗L , where |0 and |1 are the eigenstates of s z . By restricting ourselves to 'small' problems (L 20), we can compare all algorithms with accurate solutions based on exact diagonalizations and the Lanczos algorithm [29,30]. Notice that we measure the error in the full wavefunction, ε := ψ D (T) − U(T)ψ(0) 2 , and not on the expectation values of simple correlators whose exact evolution is known [14]. The outcome of some of the simulations is shown in figures 3-5, which we will discuss in the following paragraphs.
The first set of simulations was done for eight spins using matrices of size D = 16 and a DMRG basis of similar size. Since S D contains all possible states, P D = 1, we expect no truncation errors in any of the algorithms. As shown in figure 3 for the XY model with θ = 0.35 and = 0, for medium to long time steps most errors show the expected behaviour. Thus, the errors of the Trotter methods of second and fourth order follow the laws O( t 2 ) and O( t 4 ). Runge-Kutta, Taylor and Padé approximations have an error of O( t 2 ), and for the Arnoldi and Lanczos methods with N v = 4 and 8 vectors we have a qualitative behaviour O( t N v −1 ). Since the size of the matrices and of the DMRG basis is very big, all these laws are followed, irrespective of whether the implementation uses DMRG or vMPS. The only exception seems to be the DMRG Lanczos when implemented with only two target states. This method behaves more poorly than the counterpart with N v target states given by |φ k : Out of the methods that work as expected, all of them break the ideal laws at some point, acquiring an error of order O( t −2 ). This error is exponential in the number of steps T/ t and it signals the finite accuracy of the optimization algorithms, due to the limited precision of the computer. Roughly, since current computers cannot compute the norms of vectors, ψ 2 , with a relative error better than 10 −16 , a worst case estimate is = (T/ t) 2 10 −16 , which perfectly fits these lines.
Theoretically, the performance measured in computation time and memory use of all algorithms is of order O (N v N H D 3 /N), where N H is related to the number of operators in {X, Y, X † X, Y † Y, . . .}, N is the size of the problems and the additional factor N v only appears for Arnoldi and Lanczos methods. For periodic boundary conditions the cost increases by O(D 2 ). However, we have explicitly avoided PBC problems because DMRG cannot handle them so accurately. In practice, out of all methods, the Padé expansions are the slowest ones, while the Trotter formulae, being local, are the fastest. In between we find the Arnoldi and Lanczos ∆ ∆t  methods, which are nevertheless competitive if we consider their accuracy and the fact that they allow for longer time steps.
In the remaining simulations we dropped the methods from subsections 4.3 and 4.4, because they have a similar computational cost and worse performance than the vMPS Arnoldi method. We keep, on the other hand, all Trotter and DMRG methods and continue our study with bigger problems in which the MPS spaces and the DMRG basis are smaller than the limit required to represent all states accurately, D = d N/2 and M = d N/2−1 . More precisely, we choose N = 16 and 20 spins and try with M, D = 32, 64, 80 and 128. Now that the methods are potentially inexact, our previous error laws are modified by the introduction of truncation errors. The first thing we notice in figures 4(a)-(c) is that the Trotter error of second order is rather stable, its error being the same for vMPS and DMRG. However, when we go for higher orders and increase the number of exponentials, truncation errors affect more strongly the Vidal and DMRG implementations than the vMPS one. This shows the difference between making local truncations after each bond unitary (DMRG and Vidal) versus delaying truncations until the end and using the optimal projection [5].
Another important conclusion is that all DMRG methods are very sensitive to truncation errors. As shown in figures 4(a)-(c), there is not a very big difference in accuracy between the DMRG Runge-Kutta and the Lanczos implementation, and both methods are not more accurate than the DMRG Forest-Ruth formula. The main reason why higher order methods from subsection 4.2 do not improve the results is due to estimating the evolution with the matrix of the Hamiltonian on the truncated DMRG basis. Another reason is that under truncation the Lanczos recurrence no longer produces an orthogonal set of vectors.
To prove those statements we compared with a Lanczos method implemented with MPS as explained in subsection 4.5. This method does indeed have a better accuracy than the DMRG ones, which can be attributed to the use of the full Hamiltonian. The errors of the Lanczos are however larger than those of the Arnoldi method for a similar D and we have checked that under truncation the vectors of the Lanczos basis are not truly orthogonal. These small errors accumulate for shorter time steps, and only the Arnoldi method can correct them.
As for the Arnoldi method, it has the greatest accuracy and seems to be stable as D becomes small. For very small matrices the error remains constant as we decrease the time step, but this is only because there is a lower bound in the approximation error given by P D |ψ(t) − |ψ(t) 2 . Figure 6(a) 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. This plot illustrates the fact that all errors in this method are due to the final truncation.
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 and accuracy, the two winning methods are variational MPS algorithms using either the fourth order Forest-Ruth decomposition or the Arnoldi basis. The last method, however, has two advantages. One is that it can deal with nonlocal interactions and the second one is its potential for parallelizability, roughly O(N v N H /L), which most other presented algorithms lack. This can make it competitive with, for instance,  increasing the size of the matrices in the Forest-Ruth method. Regarding DMRG methods, we find that they give comparable results only for a big basis. When truncation errors pop in, their behaviour is less predictable and it does not seem worth going with more elaborate algorithms (Lanczos, Runge-Kutta) versus an ordinary Trotter formula.
The previous results have been confirmed by a systematic scanning of all possible Hamiltonians in equation (32). A selection is shown in figure 7. The Arnoldi method is shown to be accurate, even for gapless problems. When the Arnoldi method fails it is due to truncation errors. In those situations the evolved state cannot be accurately represented by MPS (and for that matter also not by DMRG), simply because (1 − P D )|ψ(t) 2

is finite and large
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT ( figure 3(d)). Increasing the number of Arnoldi vectors will also not help (see figure 6(b)) for the same reason.

Simulation of Feschbach resonances
As a real-world application, we have used the vMPS 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 [16]- [18]. 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 [32], with a coupling to a molecular state [33] 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 mm → ∞), 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 tunnelling 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 equation (33) 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 [34,35], 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 two-level 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 figure 8(b)).
We have simulated numerically the ramping of small lattices, L = 10-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 ramptime 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.  (34). The case J = 0 uses an initial condition J = 0.1 and then switches off tunnelling before ramping. (b) Correlations of the molecular state, m † i m j t=T (circles) after the ramp, and those of the initial atomic state, a † i↓ a † i↑ a j↑ a j↓ t=0 (solid line). The plot for U ↑↓ = 5J has been shifted up by 0.2.
The main conclusions is that, as shown in figure 9, indeed the correlations of the molecules are almost those of the initial state of the atoms (34), 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.
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT

Conclusions
We have performed a rather exhaustive comparison of different methods for simulating the evolution of big, 1D quantum systems [4,5], [8]- [14] with three other methods developed in this study. We find the vMPS methods to be optimal both in accuracy and performance within the formulae of similar order. All procedures 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 very large time steps is an Arnoldi method developed in this study. 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 [34,35]. Simple generalizations of this study will allow us in the future to analyse losses and creation of strongly correlated states with the help of the molecular component.
As a possible outlook, we envision the possibility of developing new algorithms in which the state is approximated by a linear combination of MPS at all times. This should be more efficient than increasing the size of the matrices, and could support distributed computations in a cluster.