Dynamical Simulations of Schrödinger's Equation via Rank-Adaptive Tensor Decompositions

Dynamical Sim ulations of Sc hr¨ odinger’s Equation via Rank-Adaptiv e T ensor Decomp ositions N. Anders P etersson ∗ 1 , Chase Ho dges-Heilmann 2 , and Stefanie G ¨ un ther 3 1,3 Cen ter for Applied Scientific Computing, La wrence Liv ermore National Lab oratory 2 Departmen t of Mathematics, Universit y of New Mexico Marc h 14, 2026 Abstract W e study lo w-rank tensor metho ds for the numerical solution of Sc hr¨ odinger’s equation with time-indep enden t and explicitly time-dep enden t Hamiltonians, motiv ated by large-scale sim ulations of man y-bo dy quan tum systems and quantum computing devices sub ject to time- dep enden t con trol pulses. W e outline the recent application of the ”basis up date and Galerkin” (BUG) metho d for tensor trains, and describe the established TDVP and TDVP-2 algorithms based on the time-dep endent v ariational principle. F or comparison, w e also consider the BUG metho d in the T uc ker format. All these approac hes enable memory efficient represen tations of partially en tangled quantum states and thereb y mitigate the exp onen tial cost of conv entional state-v ector form ulations. The rank-adaptivit y relies on the truncated singular v alue decompo- sition, in whic h the rank of a matrix is reduced by setting its smallest singular v alues to zero, based on a threshold parameter that con trols the truncation error. Numerical exp eriments on represen tative time-indep enden t and time-dep enden t Hamiltonian mo dels quan tify the tradeoff b et w een accuracy and compression across metho ds, with particular attention to the in terplay b et w een the time-step and the truncation threshold, and how the computational effort scales with the num b er of sub-systems in the quan tum system. Keywor ds— Low-rank tensor methods, Quantum dynamics, Matrix pro duct states 1 In tro duction Classical simulations of man y-b ody quan tum systems and quantum computing devices generally b ecome in tractable as the num ber of coupled sub-systems (qubits) increase. This is due to the exp onen tial gro wth of the quantum state v ector and the asso ciated increase in computational effort, as the n umbers of sub-systems increase. F or example, in the case of a quantum system consisting of N tw o-level sub-systems, the state v ector has 2 N elemen ts, and the Hamiltonian is an op erator of size 2 N ˆ 2 N . Storing and manipulating these ob jects becomes intractable, ev en for relatively mo dest n umbers of sub-systems. An alternativ e p erspective is to represen t the quantum state as a high-dimensional tensor, | ψ y P C d 1 ˆ d 2 ˆ¨¨¨ˆ d N , where each mo de corresp onds to a sub-system. While the full tensor represen tation contains the same n umber of entries as the state v ector, it enables the use of tensor decomposition techniques that allo w for a compressed represen tation of the state. Compact parameterizations of the quantum state based on a net w ork of low-rank tensors, e.g., matrix- pro duct states (MPS) (also known as tensor-trains), hav e b een used in several quan tum simulation appli- cations. F or example, MPS metho ds w ere used in gate-based circuit simulations to ev aluate observ ables in ∗ p etersson1@llnl.gov 1 a kick ed Ising mo del on 127 qubits [1]. Another application of MPS is the recen t classical computation of the ground state energy , to c hemical precision, in the 76-orbital/152-qubit F eMo-cofactor model, based on the DMRG metho d [2]. In terms of Hamiltonian simulations, the MPS representation has also emerged as a p o w erful to ol for simulating the dynamics of quan tum many-bo dy systems, providing an efficient w a y to rep- resen t en tangled quan tum states in large Hilb ert space dimensions. Methods such as the time-evolving block decimation metho d (TEBD) [3] and algorithms based on the time-dep enden t v ariational principle [4, 5, 6, 7] ha ve b een dev elop ed to study the dynamical b eha vior of such systems, where the state is represented as an MPS. Existing approaches ha ve focused on simulating the quan tum dynamics by solving the time-dep enden t Sc hr¨ odinger equation for a giv en time-independent Hamiltonian, e.g., a linear vibronic model including 252 normal modes [8]. In this pap er w e broaden the application of tensor decomp osition methods to dynamical simulations of Sc hr¨ odinger’s equation when the Hamiltonian is time-dep endent, e.g., to study quantum computing devices sub ject to time-dep enden t control pulses. Accurate dynamical simulation is a k ey comp onen t for optimal con trol in quan tum computing and quan tum error-correction, e.g., to perform state-to-state transformations, or to implemen t unitary gates. By concatenating the control pulses in time for all individual gates in a quan tum algorithm, dynamical simulations will also allo w the ov erall fidelity of an algorithm to b e predicted when applied to a quan tum device. In the following w e fo cus on time-in tegration of tw o tensor decomp ositions: T uck er tensors and tensor trains. W e describ e the nascen t application of the ”basis update and Galerkin” (BUG) metho d for tensor trains [9], as w ell as the classical TD VP and TDVP-2 algorithms [7]. F or T uck er tensors we use the BUG metho d developed by Ceruti et al. [10]. The rank-adaptivity in the TDVP-2 and BUG algorithms relies on the truncated singular v alue decomp osition (SVD), in whic h the rank of a matrix is reduced by setting its smallest singular v alues to zero. The rank is determined suc h that the error due to the truncation is smaller than a threshold parameter ϵ ! 1. T o minimize the amount of storage required b y the tensor decomp osition, it is desirable to make the ϵ -threshold as large as possible. At the same time, ϵ must b e sufficiently small to main tain the accuracy of the time-integration. In n umerical exp erimen ts we study the interpla y b etw een the time-step and the ϵ -threshold, with respect to both accuracy and storage requiremen ts. W e also ev aluate ho w the computational effort scales with the n umber of sub-systems in the quan tum system. Related w orks include the developmen t of tensor-based methods for solving differential equations using alternativ e decomp ositions, such as tree-tensor netw orks [11, 12], impro vemen ts in the accuracy of existing tensor evolution algorithms [13], and low-rank approaches for open quantum systems gov erned by the Lind- blad equation [14]. In parallel, ongoing research aims to improv e the efficiency , accuracy , and scalability of quan tum optimal con trol methods [15, 16]. The remainder of the paper is structured as follo ws. The go v erning equations are stated in Section 2 and the con v entional matrix-vector approac h for solving Schr¨ odinger’s equation is described in Section 3. The tensor train decomposition and the time-in tegration algorithms TDVP and TDVP-2 are presented in Section 4. The BUG algorithm for T uc ker tensors and tensor trains is outlined in Section 5, follo wed b y n umerical exp erimen ts presen ted in Section 6. Concluding remarks are given in Section 7. 2 Quan tum sim ulations W e consider the dynamics of a comp osite quantum system, consisting of N sub-systems, where the state in eac h sub-system b elongs to a Hilb ert (complex v ector) space H k of dimension d k ě 1, for k P r 1 , N s . The comp osite state then belongs to the Hilb ert space H “ H 1 b ¨ ¨ ¨ b H N “ C d 1 ˆ¨¨¨ˆ d N of dimension N tot “ ś N k “ 1 d k . The dynamics of the system is go verned by Schr¨ odinger’s equation, | 9 ψ p t qy “ ´ i ˆ H p t q| ψ p t qy , (2.1) where ˆ H p t q denotes the Hamiltonian op erator and | ψ p t qy is the quantum state at time t . The Hamiltonian is a Hermitian operator, implying that Schr¨ odinger’s equation is time-reversible and that the norm of the state is conserv ed in time. Given an orthonormal basis in eac h sub-system t| σ k yu , the state is represented 2 b y | ψ y “ d 1 ÿ σ 1 “ 1 ¨ ¨ ¨ d N ÿ σ N “ 1 ψ σ 1 ,...,σ N | σ 1 ¨ ¨ ¨ σ N y , (2.2) where ψ P C d 1 ˆ¨¨¨ˆ d N is an an order- N tensor. The Hamiltonian is a linear op erator that maps a state in H on to another state in the same space, ˆ H : H Ñ H . F ormally , ˆ H P B p H q “ H b H ˚ , where H ˚ is the dual space of H . The Hamiltonian is represented by ˆ H “ d 1 ÿ σ 1 “ 1 ¨ ¨ ¨ d N ÿ σ N “ 1 d 1 ÿ σ 1 1 “ 1 ¨ ¨ ¨ d N ÿ σ 1 N “ 1 H σ 1 ,...,σ N ; σ 1 1 ,...,σ 1 N | σ 1 ¨ ¨ ¨ σ N yx σ 1 1 ¨ ¨ ¨ σ 1 N | , (2.3) where H P C d 2 1 ˆ¨¨¨ˆ d 2 N , noting that b oth v ector and dual indices are subscripts in (2.3). The action of the Hamiltonian operator on a state, | ϕ y “ ˆ H | ψ y follo ws b y tensor contraction, | ϕ y “ d 1 ÿ σ 1 “ 1 ¨ ¨ ¨ d N ÿ σ N “ 1 ϕ σ 1 ,...,σ N | σ 1 ¨ ¨ ¨ σ N y , ϕ σ 1 ,...,σ N “ ÿ σ 1 1 ,...,σ 1 N H σ 1 ,...,σ N ; σ 1 1 ,...,σ 1 N ψ σ 1 1 ,...,σ 1 N . (2.4) In the following w e fo cus on a Hamiltonian describing a quantum device with fixed transition frequencies and fixed couplings, including time-dep endent con trol pulses. In a (global) frame rotating with angular frequency ω d , the Hamiltonian reads ˆ H p t q “ ˆ H s ` ˆ H c p t q , (2.5) ˆ H s “ N ÿ k “ 1 ˜ p ω k ´ ω d q ˆ a : k ˆ a k ´ ξ k 2 ˆ a : k ˆ a : k ˆ a k ˆ a k ` ÿ ℓ ą k J kℓ p ˆ a : k ˆ a ℓ ` ˆ a k ˆ a : ℓ q ¸ , (2.6) ˆ H c p t q “ N ÿ k “ 1 p k p t qp ˆ a k ` ˆ a : k q ` iq k p t qp ˆ a k ´ ˆ a : k q . (2.7) Here, ˆ a k is the low ering operator, ω k and ξ k denote the transition frequency and anharmonicit y , in sub- system k P r 1 , N s . The coupling strength betw een sub-system k and ℓ is J kℓ . Note that this mo del includes the time-independent Heissen b erg and Ising Hamiltonians as special cases. In the curren t study w e assume that the sub-systems (qubits) are arranged in a linear chain and we consider the dynamics due to nearest neigh b or in teractions. 3 Matrix-v ector form ulation The quan tum state can b e represented in v ectorized form,  ψ : “ vec p ψ q ,  ψ α “ ψ k 1 ,...,k N , α “ r L p k 1 , . . . , k N ; d 1 , . . . , d N q , (3.1) where r L maps the tensor tuple indices p k 1 , . . . , k N q , where k j P r 1 , d j s , to a linear index using lexicographical ordering (last tuple index v aries the fastest). The in verse of the index map is denoted p k 1 , . . . , k N q “ r T p α ; d 1 , . . . , d N q , see e.g. Ballard and Kolda [17]. In the following, the dimensions of the sub-systems will b e suppressed to simplify the notation. In v ectorized form, Schr¨ odingers equation is d dt  ψ “ ´ i H  ψ , t ě 0 ,  ψ p 0 q “  ψ 0 , (3.2) where H P C N tot ˆ N tot is the Hamiltonian matrix, with N tot “ ś N k “ 1 d k . Corresp onding to (2.4), the action of the Hamiltonian b ecomes  ϕ α “ N tot ÿ β “ 1 H α,β  ψ β , H α,β “ H σ 1 ,...,σ N ; σ 1 1 ,...,σ 1 N , α P r 1 , N tot s , (3.3) 3 with p σ 1 , . . . , σ N q “ r T p α q and p σ 1 1 , . . . , σ 1 N q “ r T p β q . When ev aluating the Hamiltonian mo del (2.6)-(2.7), the lo wering op erator satisfies ˆ a k “ I d 1 b ¨ ¨ ¨ b I d k ´ 1 b L r k s b I d k ` 1 b ¨ ¨ ¨ b I d N , (3.4) with L r k s “ » — — — — — – 0 1 0 ? 2 . . . . . . 0 ? d k ´ 1 0 fi ffi ffi ffi ffi ffi fl P R d k ˆ d k . (3.5) Note that the lexicographical ordering of  ψ is consistent with defining b as the Kroneck er product. When the Hamiltonian matrix is constan t in time, the solution of Sc hr¨ odinger’s equation in vectorized form (3.2) can alw ays b e obtained by matrix exp onen tiation,  ψ p t q “ exp p´ it H q  ψ 0 . (3.6) The computational complexit y of matrix exp onen tiation is O p N 3 tot q op erations, making this an efficient wa y of directly integrating Schr¨ odinger’s equation o ver long time domains, without time-stepping. How ev er, this approac h is only tractable when N tot is not to o large. F or larger system sizes, or when H dep ends on time, w e must rely on n umerical time-stepping to appro ximately solv e Schr¨ odinger’s equation. F or this purp ose w e define a grid in time, t k “ k δ , for k “ 0 , 1 , 2 , . . . , where δ ą 0 is the time-step. W e denote the numerical appro ximation of the state by  ψ k «  ψ p t k q . An unconditionally stable, norm-conserving and time-reversible sc heme for in tegrating Sc hr¨ odinger’s equation is the implicit mid-point rule (IMR):  ψ k ` 1 “  ψ k ` δ  g k ` 1 { 2 , k “ 0 , 1 , 2 , . . . ,  ψ 0 “  ψ p 0 q , (3.7)  g k ` 1 { 2 “ ´ i 2 H k ` 1 { 2 `  ψ k `  ψ k ` 1 ˘ , (3.8) where H k ` 1 { 2 “ H p t k ` δ { 2 q . By substituting  ψ k ` 1 from (3.7) in to the righ t hand side of (3.8) w e arrive at a linear system for  g k ` 1 { 2 , ` I ` iδ 2 H k ` 1 { 2 ˘  g k ` 1 { 2 “  b k ,  b k “ ´ i H p t k ` 1 { 2 q  ψ k . (3.9) Here, the linear system in (3.9) can b e solv ed iterativ ely by , e.g., Jacobi’s metho d. Note that iterative metho ds only require the action of the Hamiltonian to b e ev aluated. This approach requires less storage and is often significantly more efficient compared to first ev aluating the Hamiltonian matrix and then applying it to a v ector, esp ecially when the Hamiltonian dep ends on time. Ev en though the IMR metho d is unconditionally stable, the time-step must chosen to b e sufficien tly small to giv e an accurate solution, and it is advisable to test the numerical accuracy with, e.g., Richardson extrap olation. As an alternativ e to IMR, Sc hr¨ odinger’s equation can also be solved by an explicit time- in tegration metho d, e.g., the classical RK-4 method. Similar to solving the linear system (3.9) iteratively , explicit methods also only require the action of the Hamiltonian to be ev aluated. F or Sc hr¨ odinger’s equation, the stabilit y region of the explicit method must include the imaginary axis and the time-step m ust be sufficien tly small to ensure stabilit y . The stabilit y restrictions in an explicit metho d, and the accuracy requiremen ts in an implicit method, are closely related because they both dep end on the largest eigenv alue (in magnitude) of the Hamiltonian. Because of its unconditional stabilit y , in this w ork w e use the IMR metho d together with Richardson extrap olation for estimating the accuracy in the solution. 4 T ensor trains A tensor train (TT), also called a matrix pro duct state (MPS), represen ts an order- N tensor as a pro duct of N order-3 tensors [7]. Consider a tensor ψ P C d 1 ˆ d 2 ˆ¨¨¨ˆ d N . The elements in the tensor-train decomp osition 4 Figure 1: Left: A tensor diagram of a tensor train (MPS) of an order-4 tensor of size d 1 ˆ d 2 ˆ d 3 ˆ d 4 , is represen ted by four cores (circles), each of order three. The vertical legs represen t the ph ysical indices and the horizontal legs connecting the cores corresp ond to con tractions o v er virtual indices. The sizes of the first and last cores are 1 ˆ d 1 ˆ b 1 and b N ´ 1 ˆ d N ˆ 1, resp ectiv ely; horizontal legs are suppressed on all singleton dimensions. Right: the con traction of an MPS by an order 4 ˆ 4 matrix pro duct op erator (MPO), represen ted b y four (square) cores of order four. of ψ can be written as ψ σ 1 ,...,σ N “ b 0 ÿ m 0 “ 1 ¨ ¨ ¨ b N ÿ m N “ 1 M r 1 s p m 0 , σ 1 , m 1 q ¨ ¨ ¨ M r N s p m N ´ 1 , σ N , m N q , σ k P r 1 , d k s . (4.1) The factors M r k s P C b k ´ 1 ˆ d k ˆ b k are order-3 tensors that are referred to as cores. F urther, b k ě 1 (with b 0 “ b N “ 1) is called the b ond dimension b et ween the cores M r k s and M r k ` 1 s . The indices σ k are called ph ysical indices and b k are called virtual indices. F or eac h fixed physical index, we may define a matrix M σ k r k s “ M r k s p : , σ k , : q , σ k P r 1 , d k s , k P r 1 , N s . Then, eac h elemen t in the tensor ψ may b e written as a matrix pro duct, ψ σ 1 ,...,σ N “ M σ 1 r 1 s M σ 2 r 2 s ¨ ¨ ¨ M σ N r N s . (4.2) Note that M σ 1 r 1 s is a row vector b ecause b 0 “ 1 and M σ N r N s is a column v ector b ecause b N “ 1. As a result, the matrix product on the righ t hand side of (4.2) ev aluates to a scalar. In the following, indices may b e suppressed if they are clear from the con text, e.g., ψ “ M r 1 s ¨ M r 2 s ¨ ¨ ¨ M r N s , (4.3) where a cen tered dot ( ¨ ) represents con traction of tensors o ver common indices. T o further simplify the notation, it can b e conv enien t to represen t a tensor train in diagrammatical format, see Figure 1. It is kno wn that an y tensor of order N can b e represen ted as a tensor train with N cores, for example, using the TT-SVD algorithm [17]. F or a general tensor the bond dimensions m ust gro w exponentially tow ards the cen ter of the train, resulting in the same exp onen tial growth of storage as in the v ectorized represen tation. Ho wev er, if the b ond dimensions can be b ounded b y b k ď B , the total storage for ψ b ecomes O p N B 2 q , which is linear in N , if B can b e b ounded independently of N . The MPS representation of a given tensor is not unique b ecause it is sub ject to gauge in v ariance. F or example, we can insert the matrix iden tity C ´ 1 C “ I betw een the cores k and k ` 1 in (4.2), resulting in the mo dified cores Ă M σ k r k s “ M σ k r k s C ´ 1 and Ă M σ k ` 1 r k ` 1 s “ C M σ k ` 1 r k ` 1 s , for any non-singular matrix C P C b k ˆ b k . The gauge in v ariance can be exploited to impro ve the n umerical stabilit y of algorithms and simplify tensor con tractions. In the following we will use t w o canonical forms of a tensor train: ψ σ 1 ,...,σ N “ # A σ 1 r 1 s ¨ ¨ ¨ A σ j ´ 1 r j ´ 1 s M σ j r j s B σ j ` 1 r j ` 1 s ¨ ¨ ¨ B σ N r N s , (mixed canonical) , A σ 1 r 1 s ¨ ¨ ¨ A σ j r j s C r j s B σ j ` 1 r j ` 1 s ¨ ¨ ¨ B σ N r N s , (b ond-cen tered canonical) . (4.4) 5 The canonical forms are defined by the left- and right-normalization conditions of the cores, d k ÿ σ k “ 1 ´ A σ k r k s ¯ : A σ k r k s “ I b k , (left-normalized) , d k ÿ σ k “ 1 B σ k r k s ´ B σ k r k s ¯ : “ I b k ´ 1 , (right-normalized) . (4.5) Here and in the following, left-normalized cores are denoted b y A r j s , righ t-normalized cores by B r j s , and general cores b y M r j s . In the mixed canonical form, M r j s is called the orthogonalit y cen ter. In the b ond- cen tered canonical form, C r j s P C b j ˆ b j is called the b ond-matrix. A quan tum state on mixed canonical form can be written as | ψ y “ ÿ σ j ,m j ´ 1 ,m j | ψ L j ´ 1 p m j ´ 1 qy b ´ M σ j r j s p m j ´ 1 , m j q| σ j y ¯ b | ψ R j ` 1 p m j qy , (4.6) where the left- and righ t-state v ectors are defined b y | ψ L j p m j qy “ ÿ σ 1 ,...,σ j A σ 1 r 1 s ¨ ¨ ¨ A σ j r j s p : , m j q| σ 1 ¨ ¨ ¨ σ j y , m j P r 1 , b j s , (4.7) | ψ R j ` 1 p m 1 j qy “ ÿ σ j ` 1 ,...,σ N B σ j ` 1 r j ` 1 s p m 1 j , : q ¨ ¨ ¨ B σ N r N s | σ j ` 1 ¨ ¨ ¨ σ N y , m 1 j P r 1 , b j s . (4.8) Eac h set of vectors has m utually orthogonal elemen ts, as is made precise in the follo wing lemma. Lemma 1. Assume that the b asis ve ctors t| σ k yu N k “ 1 ar e mutual ly ortho gonal and normalize d. L et the set of left state ve ctors t| ψ L j p m j qyu b j m j “ 1 b e define d by (4.7) , wher e the or der-3 tensors A r k s satisfy the left- normalization c ondition in (4.5) . Then the left state ve ctors satisfy the ortho gonality c ondition x ψ L j p α q| ψ L j p β qy “ δ α,β , α , β P r 1 , b j s . (4.9) Similarly, let the set of right state ve ctors t| ψ R j p m 1 j qyu b j m 1 j “ 1 b e define d by (4.8) , wher e the or der-3 tensors B r k s satisfy the right-normalization c ondition in (4.5) . Then the right state ve ctors satisfy the ortho gonality c ondition x ψ R j p α q| ψ R j p β qy “ δ α,β , α , β P r 1 , b j s . (4.10) Pr o of. F ollows directly from the left- and righ t-normalization conditions on A r j s and B r k s . The orthogonalit y center in the mixed canonical form (4.6) can b e mo ved from M r j s to the b ond betw een the cores j and j ` 1. The first step is to factor M r j s in to a left-normalized core and a square bond-matrix C r j s P C b j ˆ b j . F actorization is ac hiev ed by reshaping M r j s in to a matrix M L j P C b j ´ 1 ¨ d j ˆ b j , follow ed by a singular-v alue decomp osition M L j “ : U j S j V : j . A left-orthogonal core is constructed b y in verse reshaping U j P C b j ´ 1 ¨ d j ˆ b j in to the tensor A r j s P C b j ´ 1 ˆ d j ˆ b j . Here, the left-orthogonalit y condition follo ws from U : j U j “ I b j . The bond-matrix is defined by C r j s “ S j V : j P C b j ˆ b j , resulting in M σ j r j s p m j ´ 1 , m j q “ b j ÿ m 1 j “ 1 A σ j r j s p m j ´ 1 , m 1 j q C r j s p m 1 j , m j q . (4.11) This decomposition results in the bond-centered canonical representation of the state, | ψ y “ ÿ m j ,m 1 j C r j s p m j , m 1 j q| ψ L j p m j qy b | ψ R j ` 1 p m 1 j qy . (4.12) The decomposition (4.12) is closely related to the Sc hmidt decomp osition. The rank of the matrix C r j s is called the Sc hmidt rank and determines the en tanglement b et ween the states to the left and right of the b ond b et ween sites j and j ` 1. 6 By contracting the bond-matrix C r j s with the core on the right, we can mo ve the orthogonality center one step to the righ t, M σ j ` 1 r j ` 1 s p m j , m j ` 1 q “ b j ÿ m 1 j “ 1 C r j s p m j , m 1 j q B σ j ` 1 r j ` 1 s p m 1 j , m j ` 1 q . (4.13) This op eration is an essen tial comp onen t in the TD VP algorithm (presen ted b elo w) for solving Schroedinger’s equation. 4.1 Matrix pro duct op erators The tensor-train representation of a state can b e generalized to represen t the elements of the Hamilto- nian (2.3) as a matrix pro duct op erator (MPO), ˆ H σ 1 ,...,,σ N ; σ 1 1 ,...σ 1 N “ b 1 0 ÿ w 0 “ 1 ¨ ¨ ¨ b 1 N ÿ w N “ 1 W r 1 s p w 0 , σ 1 , σ 1 1 , w 1 q ¨ ¨ ¨ W r N s p w N ´ 1 , σ N , σ 1 N , w N q . (4.14) Here, eac h core W r k s P C b 1 k ´ 1 ˆ d k ˆ d k ,b 1 k is an order 4 tensor. The first and last bond dimensions are b 1 0 “ b 1 N “ 1, and the interior b ond dimensions satisfy b 1 k ě 1 for k “ 1 , . . . , N ´ 1. The in terior bond dimensions dep end on the specific Hamiltonian model. F or example, b 1 k “ 2 suffices if the Hamiltonian only has lo cal in teractions within eac h sub-system. F or all Hamiltonians with short-range interactions in one dimension, the required MPO b ond dimension is small p« 5 q and is indep enden t of the num b er of sub-systems [18]. Similar to the MPS, for each pair of ph ysical indices p σ k , σ 1 k q w e ma y define a matrix, W σ k ,σ 1 k r k s “ W r k s p : , σ k , σ 1 k , : q P C b 1 k ´ 1 ˆ b 1 k , σ k , σ 1 k P r 1 , d k s , k P r 1 , N s . (4.15) The MPO in matrix form allows the action of the Hamiltonian on a state to b e expressed as an MPS. This op eration is shown as a tensor diagram in Figure 1. Let | ψ y b e giv en as an MPS with general cores as in (4.2) and consider | ϕ y “ ˆ H | ψ y . The resulting MPS for | ϕ y can b e obtained as the tensor pro duct of the individual tensors, with elemen ts ϕ σ 1 ,...,σ N “ ÿ m 0 ,...m N ÿ w 0 ,...,w N Ă M σ 1 r 1 s ; p w 0 m 0 q , p w 1 m 1 q ¨ ¨ ¨ Ă M σ N r N s ; p w N ´ 1 m N ´ 1 q , p w N m N q (4.16) Ă M σ k r k s ; p w k ´ 1 m k ´ 1 q , p w k m k q “ ÿ σ 1 k W r k s p w k ´ 1 , σ k , σ 1 k , w k q M r k s p m k ´ 1 , σ 1 k , m k q , (4.17) where the indexing of the core tensors in the resulting MPS is Ă M σ k r k s ; p w k ´ 1 m k ´ 1 q , p w k m k q : “ Ă M r k s p m 1 k ´ 1 , σ k , m 1 k q , with m 1 k “ p w k ´ 1 q b k ` m k . The state | ϕ y is also in the MPS format, but with inflated b ond dimensions: b 2 k “ b 1 k ¨ b k . Thus, rep eated application of the Hamiltonian to a state quic kly increases the storage requiremen ts for the MPS. Remark 1. The inflate d MPS in (4.16) - (4.17) is never explicitly forme d in the time-inte gr ation metho ds we describ e b elow. 4.2 The time-dep enden t v ariational principle This section gives an ov erview of the tim e-dependent v ariational principle and the asso ciated pro jection and contraction op erators. These operators are used in the TDVP and TD VP-2 algorithms, which will b e describ ed in the following sections. V ariational approaches for solving the time-dep endent Schroedinger equation date bac k to the works by Dirac [19], F renkel [20], McLac hlan [21], and the time-dep enden t v ariational principle analyzed by Kramer & Saraceno [22] and Broeckho v e et al. [23]. In the v ariational approac h, the ev olution of the state is constrained 7 to a manifold M of matrix pro duct states where the b ond dimensions are fixed, resulting in the pro jected Sc hro edinger equation, | 9 ψ y “ ´ i ˆ P T , | ψ y ˆ H p t q| ψ y . (4.18) Here, the action of the Hamiltonian on the state is pro jected on to the tangent plane of the manifold, where the tangen t plane is centered at the current state | ψ y . The state-dep enden t pro jector can be written as a sum of lo cal pro jection operators [7], ˆ P T , | ψ y “ N ÿ j “ 1 ˆ P L,ψ j ´ 1 b ˆ I j b ˆ P R,ψ j ` 1 ´ N ´ 1 ÿ i “ 1 ˆ P L,ψ j b ˆ P R,ψ j ` 1 , (4.19) where ˆ P L,ψ j ´ 1 and ˆ P R,ψ j are constructed from the gauge-fixed left- and righ t-normalized cores in the mixed- canonical (4.6) and b ond-cen tered (4.12) representations of the curren t state, ˆ P L,ψ j ´ 1 “ d j ´ 1 ÿ m j ´ 1 “ 1 | ψ L j ´ 1 p m j ´ 1 qyx ψ L j ´ 1 p m j ´ 1 q| b ˆ I R j , (4.20) ˆ P R,ψ j ` 1 “ ˆ I L j b d j ÿ m j “ 1 | ψ R j ` 1 p m j qyx ψ R j ` 1 p m j q| . (4.21) Here, ˆ I R j and ˆ I L j are the identit y operators on sites k ě j and k 1 ď j , resp ectiv ely . When the state is on mixed canonical form with the orthogonalit y center at core j , it is straightforw ard to v erify that ˆ P L,ψ j ´ 1 | ψ y “ | ψ y and ˆ P R,ψ j ` 1 | ψ y “ | ψ y . If we instead write the same state in the bond-centered canonical form, w e easily see that ˆ P L,ψ j | ψ y “ | ψ y and ˆ P R,ψ j ` 1 | ψ y “ | ψ y . Since the same state can b e written in either mixed canonical form (with the orthogonalit y center at any core), or in b ond-cen tered form (where the b ond-matrix is b et ween any t wo consecutiv e cores), we conclude that every term in (4.19) is a pro jection op erator. W e pro ceed b y defining a core-centered con traction op erator for extracting the core when the state is on the mixed canonical form (4.6) with the orthogonality center at site j . Note that x ψ L j ´ 1 p m j ´ 1 q| b x σ j | b x ψ R j ` 1 p m j q| ψ y “ M σ j r j s p m j ´ 1 , m j q . W e can therefore obtain all e lemen ts of M σ j r j s b y defining the con traction op erator as an outer pro duct of left- and right-con traction operators, ˆ Q σ j j “ » — – x ψ L j ´ 1 p 1 q| . . . x ψ L j ´ 1 p b j ´ 1 q| fi ffi fl b x σ j | b ” x ψ R j ` 1 p 1 q| , ¨ ¨ ¨ , x ψ R j ` 1 p b j q| ı . (4.22) The core is then extracted by applying the con traction to the state, ˆ Q σ j j | ψ y “ » — – M σ j r j s p 1 , 1 q ¨ ¨ ¨ M σ j r j s p 1 , b j q . . . . . . . . . M σ j r j s p b j ´ 1 , 1 q ¨ ¨ ¨ M σ j r j s p b j ´ 1 , b j q fi ffi fl “ M σ j r j s , ô ˆ Q j | ψ y “ M r j s . (4.23) Another important property follo ws b y applying the core-cen tered contraction op erator to the j -th term in the first sum of (4.19), ˆ Q σ j j p m j ´ 1 , m j q ´ ˆ P L,ψ j ´ 1 b ˆ I j b ˆ P R,ψ j ` 1 ¯ “ ˆ Q σ j j p m j ´ 1 , m j q ÿ m 1 j ´ 1 ,m 1 j ´ | ψ L j ´ 1 p m 1 j ´ 1 qyx ψ L j ´ 1 p m 1 j ´ 1 q| b ˆ I j b | ψ R j ` 1 p m 1 j qyx ψ L j ` 1 p m 1 j q| ¯ “ ` x ψ L j ´ 1 p m j ´ 1 q| b x σ j | b x ψ R j ` 1 p m j q| ˘ “ ˆ Q σ j j p m j ´ 1 , m j q , (4.24) that is, ˆ Q j ´ ˆ P L,ψ j ´ 1 b ˆ I j b ˆ P R,ψ j ` 1 ¯ “ ˆ Q j . (4.25) 8 Assuming that the state is on the b ond-cen tered canonical form (4.12), with the b ond-matrix b et ween sites j and j ` 1, we hav e x ψ L j p m j q| b x ψ R j ` 1 p m 1 j q| ψ y “ C j p m j , m 1 j q . The b ond-cen tered contraction op erator then follo ws from the outer product, ˆ Q B j “ » — – x ψ L j p 1 q| . . . x ψ L j p b j q| fi ffi fl b ” x ψ R j ` 1 p 1 q| , ¨ ¨ ¨ , x ψ R j ` 1 p b j q| ı , ˆ Q B j | ψ y “ C r j s . (4.26) F urthermore, by applying the b ond-cen tered con traction op erator to the j -th term in the second sum of (4.19), w e get ˆ Q B j ´ ˆ P L,ψ j b ˆ P R,ψ j ` 1 ¯ “ ˆ Q B j . (4.27) 4.3 The TD VP algorithm The TD VP algorithm [7], which is closely related to the pro jector-splitting metho d developed b y Lubich et al. [24], is based on the time-dep enden t v ariational principle. By inserting the explicit form of the pro jection op erator (4.19) in to the pro jected Schr¨ odinger equation (4.18), w e arrive at | 9 ψ y “ ´ i N ÿ j “ 1 ´ ˆ P L,ψ j ´ 1 b ˆ I j b ˆ P R,ψ j ` 1 ¯ ˆ H p t q| ψ y ` i N ´ 1 ÿ j “ 1 ´ ˆ P L,ψ j b ˆ P R,ψ j ` 1 ¯ ˆ H p t q| ψ y . (4.28) In the TD VP algorithm, a splitting approac h is employ ed in which eac h term on the righ t hand side of (4.28) is in tegrated individually and sequen tially , $ & % | 9 ψ y “ ´ i ´ ˆ P L,ψ j ´ 1 b ˆ I j b ˆ P R,ψ j ` 1 ¯ ˆ H p t q| ψ y , | 9 ψ y “ ` i ´ ˆ P L,ψ j b ˆ P R,ψ j ` 1 ¯ ˆ H p t q| ψ y , j “ 1 , . . . , N ´ 1 , (4.29) | 9 ψ y “ ´ i ´ ˆ P L,ψ N ´ 1 b ˆ I N ¯ ˆ H p t q| ψ y . (4.30) Here w e ha ve defined ˆ P L,ψ 0 “ 1 to allow for a uniform formulation of all terms. In TD VP , the state is calculated at discrete times, t k “ k δ , where δ ą 0 is the time-step and k “ 1 , . . . , N t . Starting from the initial state written on mixed canonical form with the orthogonality center at the leftmost site ( j “ 1), the state is evolv ed in t w o half-steps based on Strang splitting. In the first half-step, the individual cores are up dated to time t 0 ` δ { 2 in a left-to-right sweep b y solving (4.29)-(4.30), resulting in a state on mixed canonical form with the orthogonalit y cen ter at the righ tmost site ( j “ N ). This is follow ed b y the second half-step in which the individual cores are ev olv ed to time t 0 ` δ in a righ t-to-left sweep b y solving by solving (4.29)-(4.30) in rev ersed order, resulting in a state in mixed canonical form with the orthogonalit y center back at the leftmost site. The time-stepping pro cedure is then rep eated un til the final time is reached. By applying the core-centered contraction (4.23) to the first equation in (4.29), w e obtain the site-local Sc hro edinger equation for core M r j s , 9 M r j s p t q “ ´ i ˆ Q j ¨ ˝ ÿ σ 1 1 ,...,σ 1 N ˆ H σ , σ 1 p t q A σ 1 1 r 1 s ¨ ¨ ¨ A σ 1 j ´ 1 r j ´ 1 s M σ 1 j r j s p t q B σ 1 j ` 1 r j ` 1 s ¨ ¨ ¨ B σ 1 N r N s | σ 1 y ˛ ‚ (4.31) Here, ˆ H p t q is expressed in MPO format as in (4.16). On the right hand side, note that only the core M r j s dep ends on time; all other factors are fixed and can b e pre-computed. W e can therefore define an effectiv e Hamiltonian acting on M r j s and write the site-lo cal equation as 9 M r j s “ ´ i ˆ H ef f j ¨ M r j s , ˆ H ef f j “ L j ´ 1 p t 0 ` δ { 2 q ¨ W r j s ¨ R j ` 1 p t 0 q , t P r t 0 , t 0 ` δ { 2 s , (4.32) 9 Figure 2: TDVP: ev aluating the action of the effectiv e Hamiltonian on the state in a 4-site system. The MPS of the state is in the top row, the MPO of the Hamiltonian in the middle row, and the resulting action is represen ted b y the open legs in the b ottom row. Left: the single-site Hamiltonian acting on the (green) core M r 2 s , connected with dashed legs. Righ t: the b ond-cen tered Hamiltonian acting on the (ligh t blue) b ond-matrix C r 1 s , connected with dashed legs b et w een the first and second core. using M r j s p t 0 q as initial condition. The effectiv e Hamiltonian ˆ H ef f j follo ws as a contraction betw een the con tributions from core W r j s in the MPO, as well as the environmen ts on the left ( L j ´ 1 ) and right ( R j ` 1 ) of the orthogonalit y center. This decomp osition is depicted as a tensor diagram in Figure 2 (left). W e note that ev aluating the action of the effectiv e Hamiltonian only requires con tractions of lo w-order tensors [7]. F urther, the right environmen ts p R N , R N ´ 1 , . . . R 2 q can b e recursively precomputed in a right-to-left sweep. The left en vironment p L j ´ 1 q is up dated based on L j ´ 2 during the left-to-right sweep. The resulting differential equation can b e solved b y IMR com bined with an Jacobi iteration solving the associated linear system, cf. (3.9). After the core M r j s has b een evolv ed to time t 0 ` δ { 2, the orthogonalit y center is mov ed one step to the righ t as in (4.13), defining the left-normalized tensor A r j s p t 0 ` δ { 2 q and the b ond-matrix C r j s p t 0 ` δ { 2 q . The TD VP metho d then pro ceeds by solving the b ond-matrix differen tial equation 9 C r j s “ ` i ˆ K ef f j ¨ C r j s , t 0 ` δ ě t ě t 0 , (4.33) bac kwards in time. Again, the effectiv e Hamiltonian ˆ K ef f j “ L j p t 0 ` δ { 2 q ¨ R j ` 1 p t 0 q can be formed b y con tractions of low-order tensors as is illustrated in Figure 2 (righ t). Given C r j s p t 0 q , the core at the next orthogonalit y center then follows from M r j ` 1 s p t 0 q : “ C r j s p t 0 q ¨ B r j ` 1 s p t 0 q . The left-to-right sweep is rep eated un til the rightmost site is reac hed, resulting in a mixed-canonical representation of the state at time t 0 ` δ { 2, with orthogonalit y cen ter at site N . The second half-step of TDVP proceeds in a similar w ay but up dates the cores in the opp osite order, starting from the righ tmost site and sweeping left, resulting in all cores being up dated to time-lev el t k ` δ . The TDVP algorithm conserves b oth the norm and the energy of the state during the time evolution [7]. Ho wev er, the accuracy of the solution dep ends on the bond dimensions that are fixed throughout the time ev olution, i.e., they need to be assigned before the state ev olution is known. Hence, the main difficulty with TD VP is to choose the b ond dimensions such that they are large enough to accurately capture the state dynamics, but small enough to av oid making the problem computationally intractable. 4.4 The TD VP-2 algorithm The TDVP-2 algorithm ov ercomes the main shortcoming of the TDVP algorithm by adaptiv ely changing the bond dimensions based on a threshold parameter, ϵ , in a truncated SVD that will b e described b elo w. Similar to TD VP , the TD VP-2 algorithm is based on a Strang splitting of the pro jected Sc hro edinger equation that evolv es the state in t wo half-steps. T o initialize, the state is first transformed on mixed 10 Figure 3: TD VP-2: ev aluating the action of the effectiv e t w o-site Hamiltonian on th e merged (green rectangular) core. The MPS of the state is in the top row, the MPO of the Hamiltonian in the middle row, and the resulting action is represen ted b y the four op en legs in the b ottom row. canonical form with the orthogonality center at the leftmost site p j “ 1 q . The first half-step in TDVP-2 pro ceeds by contracting M r j s and B r j ` 1 s to form the merged order-4 tensor M σ j ,σ j ` 1 r j,j ` 1 s p m j ´ 1 , m j ` 1 q “ ÿ m 1 j M σ j r j s p m j ´ 1 , m 1 j q B σ j ` 1 r j ` 1 s p m 1 j , m j ` 1 q , M r j,j ` 1 s “ M r j s ¨ B r j ` 1 s . (4.34) This tensor is then ev olv ed forwards in time under an effectiv e t wo-site Hamiltonian, 9 M r j,j ` 1 s “ ´ i ˆ H ef f j,j ` 1 ¨ M r j,j ` 1 s , t 0 ď t ď t 0 ` δ { 2 , (4.35) ˆ H ef f j,j ` 1 “ L j ´ 1 p t 0 ` δ { 2 q ¨ W r j s ¨ W r j ` 1 s ¨ R j ` 2 p t 0 q . (4.36) The action of the effectiv e Hamiltonian on the merged core M r j,j ` 1 s is illustrated in Figure 3. A truncated SVD is then applied to split the order-4 tensor bac k in to t w o order-3 tensors, in whic h the threshold parameter ( ε ) controls the accuracy of the splitting and determines how the bond dimension b j needs to b e mo dified (see Appendix A for details). This results in M r j,j ` 1 s p t 0 ` δ { 2 q « A r j s p t 0 ` δ { 2 q ¨ C r j s p t 0 ` δ { 2 q ¨ B r j ` 1 s p t 0 ` δ { 2 q , (4.37) M r j s p t 0 ` δ { 2 q : “ C r j s p t 0 ` δ { 2 q ¨ B r j ` 1 s p t 0 ` δ { 2 q , (4.38) i.e., the orthogonality cen ter is now at site j ` 1, but at time-level t 0 ` δ { 2. The next step is to ev olve M r j ` 1 s bac kwards in time b y solving the single-site differen tial equation 9 M r j ` 1 s “ ´ i ˆ H ef f j ` 1 ¨ M r j ` 1 s , ˆ H ef f j ` 1 “ L j p t 0 ` δ { 2 q ¨ W r j ` 1 s ¨ R j ` 2 p t 0 q , t 0 ` δ { 2 ě t ě t 0 , (4.39) resulting in M j ` 1 p t 0 q . Note that this equation is identical to (4.32) in TDVP , but it is evolv ed backw ards in time. The merged core, M r j ` 1 ,j ` 2 s p t 0 q , can then formed and the pro cess is repeated until the rightmost site is reac hed. Similar to TDVP , the second half-step updates the cores in the opp osite order in a right-to-left sw eep resulting in an appro ximation of the state at time t 0 ` δ . 5 The basis up date and Galerkin (BUG) algorithm The BUG algorithm is a time integration metho d that can be used to solv e Schr¨ odinger’s equation when the state is represen ted by a tensor decomp osition. The most basic case o ccurs for a tensor of order 2, i.e., a matrix differen tial equation. The BUG algorithm is based on a splitting approac h that av oids stiffness due to small singular v alues in the evolv ed matrix [10], making explicit time-in tegration metho ds tractable and improving the conv ergence of iterativ e solvers within an implicit metho d. The basic BUG algorithm has b een generalized to evolv e tensor differential equations where the state is represen ted b y an arbitrary 11 hierarc hical tree tensor decomp osition, see Lubich et al. [25]. While this algorithm also applies to T uck er tensors and tensor trains as sp ecial cases, it is non-trivial to decipher a practical implemen tation from this v ery general description. In this section we first outline the rank-adaptive BUG algorithms for matrices and T uc ker tensors, follo wed by the BUG algorithm for tensor trains (MPS-BUG). 5.1 BUG for matrices and T uck er tensors The basic BUG algorithm applies to matrix differen tial equations. Consider the order-2 tensor (matrix) T P C d 1 ˆ d 2 go verned by the differential equation, 9 T “ F p t, T q , t ą 0 , T p 0 q “ T 0 , (5.1) where the function F P t R ˆ C d 1 ˆ d 2 Ñ C d 1 ˆ d 2 u is assumed to be sufficien tly regular. The order-2 tensor is represented in the form T “ U S V : , with U : U “ I and V : V “ I . This resembles a singular v alue decomp osition, except that the matrix S is not assumed to b e diagonal. As before, w e denote the time-step b y δ ą 0. Starting from the T 0 “ T p t 0 q , factored as T 0 “ U 0 S 0 V : 0 , one time-step of the rank-adaptive BUG sc heme calculates T 1 « T p t 1 q , with t 1 “ t 0 ` δ , as describ ed in Algorithm 1. The algorithm is rep eated until the desired final time is reac hed. Algorithm 1 One time-step with rank-adaptive BUG for matrix differential equations. 1: Inputs: Initial state at time t 0 : T 0 “ U 0 S 0 V : 0 , SVD threshold ϵ ě 0, time-step δ . 2: Outputs: State at time t 1 “ t 0 ` δ : T 1 “ U 1 S 1 V : 1 . 3: K-step: In tegrate the differen tial equation: 9 K p t q “ F ` t, K p t q V : 0 ˘ V 0 , t 0 ă t ď t 1 , K p t 0 q “ U 0 S 0 . 4: QR-step #1: QR-factorize the concatenated matrix [ K p t 1 q| U 0 s “ : ˆ U 1 R 1 . 5: L-step: In tegrate the differen tial equation: 9 L p t q “ ´ F ` t, U 0 L p t q : ˘ ¯ : U 0 , t 0 ă t ď t 1 , L p t 0 q “ V 0 S : 0 . 6: QR-step #2: QR-factorize the concatenated matrix r L p t 1 q| V 0 s “ : ˆ V 1 R 2 . 7: S-step: In tegrate the differen tial equation: 9 ˆ S p t q “ ˆ U : 1 F ` t, ˆ U 1 ˆ S p t q ˆ V : 1 ˘ ˆ V 1 , t 0 ă t ď t 1 , ˆ S p t 0 q “ ˆ U : 1 U 0 S 0 V : 0 ˆ V 1 . 8: T runcate b ond dimension: Calculate ϵ -truncated SVD : ˆ S p t 1 q « : U 1 S 1 V : 1 . 9: State at time t 1 : T 1 “ U 1 S 1 V : 1 . W e pro ceed by discussing the BUG algorithm for T uck er tensors. An y tensor Y P C d 1 ˆ ... ˆ d N can b e represen ted in T uc ker format as a factorization of one order- N tensor G and N matrices U i [17], Y p i 1 , . . . , i N q “ b 1 ÿ j 1 “ 1 ¨ ¨ ¨ b N ÿ j N “ 1 G p j 1 , . . . , j N q U 1 p i 1 , j 1 q ¨ ¨ ¨ U N p i N , j N q . (5.2) Here, G P C b 1 ˆ ... ˆ b N is called the core, and U k P C d k ˆ b k for k P r 1 , N s , are the factor matrices. Note that the T uc ker format can only giv e a compressed representation of Y when the bond dimensions are smaller than the ph ysical dimensions, i.e., b k ă d k for k P r 1 , N s . A T uc ker-tensor decomp osition can be written as a sequence of T ensor-Times-Matrix (TTM) pro ducts [17], denoted as Y “ G ˆ 1 U 1 ¨ ¨ ¨ ˆ N U N “ G ˆ N j “ 1 U j , or Y “ v G ; U 1 , . . . , U N w . (5.3) 12 Let the T uck er tensor defined ab o v e satisfy the tensor differential equation, 9 Y p t q “ F p t, Y p t qq , t 0 ă t ď T , Y p t 0 q “ v G p t 0 q ; U 1 p t 0 q , . . . , U N p t 0 qw “ : Y 0 . (5.4) As b efore, the function F P t R ˆ C d 1 ˆ ... ˆ d N Ñ C d 1 ˆ ... ˆ d N u is assumed to be sufficien tly regular. The rank-adaptiv e BUG algorithm [26] computes a T uck er tensor Y 1 « Y p t 1 q , where Y 1 “ v G p t 1 q , U 1 p t 1 q , . . . , U N p t 1 qw . (5.5) Here, the K and L steps in Algorithm 1 are replaced b y N sub-steps for solving the corresp onding matricised differen tial equations. Here we only remark on how the algorithm is used for solving Schr¨ odinger’s equation. In this case w e need to ev aluate F p t, Y q “ ´ i ˆ H p t q Y . Using lexicographical ordering, v ec p Y q “ p U 1 b ¨ ¨ ¨ b U N q v ec p G q . In general the Hamiltonian can b e written as ˆ H “ ř j ˆ H j , where each term satisfies, ˆ H j “ S 1 b ¨ ¨ ¨ b S N , ˆ H j v ec p Y q “ p S 1 U 1 b ¨ ¨ ¨ b S N U N q vec p G q . The action of eac h term in the Hamiltonian on a T uck er tensor th us results in another T uck er tensor, with the same core and same dimensions of the factor matrices. Suc h T uc k er tensors can b e added efficien tly , e.g., b y calculating a site-lo cal basis for the factor matrices at eac h site. W e conclude that the action of the Hamiltonian on a T uck er tensor can be performed efficien tly through site-local operations. F or further details, w e refer to Ceruti et al. [26]. 5.2 BUG for tensor trains Compared to the TD VP-2 metho d [7] the BUG algorithm for tensor trains has one significant adv antage: eac h sub-step in BUG only in volv es forw ards in time evolution. While this may not b e imp ortan t for the Sc hr¨ odinger equation, whic h is time-reversible, the backw ards sub-step within TDVP-2 could lead to time- stepping instabilities for dissipativ e differen tial equations. F or example, in the Monte-Carlo wa ve function metho d where the ev olution is go v erned by Schr¨ odinger’s equation with a non-Hermitian effective Hamilto- nian [27]. In contrast to the TD VP and TD VP-2 algorithms, where the orthogonalit y center is mov ed from one core to the next, the MPS-BUG metho d keeps the orthogonality center fixed throughout the time-stepping, t ypically somewhere near the middle of the train. In the following, assume the state to b e given as an MPS at time t 0 , with the orthogonalit y cen ter at site j c , ψ σ 1 ,...,σ N p t 0 q “ A σ 1 r 1 s p t 0 q ¨ ¨ ¨ A σ j c ´ 1 r j c ´ 1 s p t 0 q M σ j c r j c s p t 0 q B σ j c ` 1 r j c ` 1 s p t 0 q . . . B σ N r N s p t 0 q , (5.6) where the tensor cores ha v e the sizes A r k s p t 0 q P C b k ´ 1 ˆ d k ˆ b k , M r j s p t 0 q P C b j ´ 1 ˆ d j ˆ b j , B r m s p t 0 q P C b m ´ 1 ˆ d m ˆ b m , (5.7) for k P r 1 , j c ´ 1 s and m P r j c ` 1 , N s , with b 0 “ b N “ 1. In each timestep of the MPS-BUG algorithm, the cores to the left of the orthogonality cen ter, A σ k r k s for k P r 1 , j c ´ 1 s , are updated in a left-to-right sweep, while the cores to the righ t of the orthogonalit y center, B σ m r m s for m P r j c ` 1 , N s , are up dated in a right-to-left sweep. These sweeps are indep enden t of each another and can b e p erformed concurrently . The core at the orthogonality center, M σ j c r j c s , is up dated after all A r k s p t 1 q and B r m s p t 1 q hav e been calculated. During these operations, all bond dimensions are temporarily doubled to accoun t for potential increases in entanglemen t. A t the end of eac h timestep, the bond dimensions are re-ev aluated based on a truncated singular v alue decomp osition of eac h core, where singular v alues smaller than ϵ are remo ved. Ev olving the state by one time-step results in another MPS, ψ σ 1 ,...,σ N p t 1 q « A σ 1 r 1 s p t 1 q ¨ ¨ ¨ A σ j c ´ 1 r j c ´ 1 s p t 1 q M σ j c r j c s p t 1 q B σ j c ` 1 r j c ` 1 s p t 1 q ¨ ¨ ¨ B σ N r N s p t 1 q , t 1 “ t 0 ` δ, (5.8) 13 Figure 4: MPS-BUG: application of ˆ H ef f 2 to r A r 2 s (green, top ro w), giving the slope ˆ H ef f 2 r A r 2 s (3 solid legs in to empty space, bottom ro w). where the cores temp orarily are assigned doubled b ond dimensions, A r k s p t 1 q P C 2 b k ´ 1 ˆ d k ˆ 2 b k , M r j c s p t 1 q P C 2 b j c ´ 1 ˆ d j ˆ 2 b j c , B r m s p t 1 q P C 2 b m ´ 1 ˆ d m ˆ 2 b m , (5.9) for k P r 2 , j c ´ 1 s and m P r j c ` 1 , N ´ 1 s . Note that the b ond dimensions at b oth ends of the tensor train are not doubled. Instead, A r 1 s p t 1 q P C 1 ˆ d 1 ˆ 2 b 1 and B r N s p t 1 q P C 2 b N ´ 1 ˆ d N ˆ 1 . Similar to the site-lo cal differen tial equation in the TD VP algorithm (4.32), the sites to the left and righ t of the orthogonalit y center are evolv ed by sequen tially by solving lo cal Sc hr¨ odinger equations based on effectiv e Hamiltonians. Sites to the left of the orthogonalit y center are up dated b y solving the site-local Sc hr¨ odinger equation: 9 Ă A r k s “ ´ i ˆ H ef f k ¨ r A r k s , ˆ H ef f k “ L k ´ 1 p t 1 q ¨ W r k s ¨ R k ` 1 p t 0 q , k P r 1 , j c ´ 1 s , (5.10) where r A r k s P C 2 b k ´ 1 ˆ d k ˆ b k (except for k “ 1 where the first mode alw ays has dimension 1). A tensor diagram for ev aluating the action of the effectiv e Hamiltonian is sho wn in Figure 4. In a corresponding wa y , sites to the righ t of the orthogonalit y center are updated b y solving the site-lo cal differen tial equation: 9 r B r m s “ ´ i ˆ H ef f m ¨ r B r m s , ˆ H ef f m “ L m ´ 1 p t 0 q ¨ W r m s ¨ R m ` 1 p t 1 q , m P r j c ` 1 , N s , (5.11) where r B r m s P C b m ´ 1 ˆ d m ˆ 2 b m (except for m “ N where the third mo de alw ays has dimension 1). Before these site-lo cal equation can be solved, we note t wo differences from the TDVP algorithm. First, the bond dimensions of the initial conditions for r A r k s p t 0 q and r B r m s p t 0 q are differen t from the b ond dimensions of A r k s p t 0 q and B r m s p t 0 q , respe ctiv ely . Secondly , the b ond dimension of r A r k s p t 1 q is differen t from A r k s p t 1 q , and the b ond dimension of r B r m s p t 1 q is different from B r m s p t 1 q . The algorithms for accurately consolidating the tensor dimensions is b ey ond the scope of this presentation, but will b e discussed in detail in a forthcoming pap er by Sulz [9]. Once A r k s p t 1 q and B r m s p t 1 q ha ve b een calculated for all k and m , the site at the orthogonality cen ter is up dated by solving a third site-lo cal Schr¨ odinger equation: 9 Ă M r j c s “ ´ i ˆ H ef f j c ¨ Ă M r j c s , ˆ H ef f j c “ L j c ´ 1 p t 1 q ¨ W r j c s ¨ R j c ` 1 p t 1 q . (5.12) After all sites of the tensor train hav e b een up dated, the state at time t 1 is represented in the form (5.8)- (5.9), where all b ond dimension hav e b een doubled. T o limit the storage requiremen ts while main taining the prescrib ed solution accuracy , a truncated singular v alue decomp osition is then performed on each core, in whic h all singular v alues smaller than ϵ are set to zero, thus allowing the bond dimensions to b e reduced accordingly . 14 6 Numerical Examples In this section w e numerically ev aluate the computational performance for solving Sc hr¨ odinger’s equation using the T uck er tensor and tensor train (MPS) decomp ositions, and compare them to the conv entional matrix-v ector approac h. In the matrix-v ector formulation, Schr¨ odinger’s equation can b e solved by matrix exp onen tiation when the Hamiltonian is time-independent and the dimension of the state v ector is sufficiently small. F or time-dep endent Hamiltonians, Sc hr¨ odinger’s equation must be solved n umerically; here w e use the C++ code Quandary [28] that implements the approac h in Section 3. The tensor train and T uck er tensor algorithms from Sections 4 and 5 are implemented in the Julia language, based on the iT ensor pack age [29]. All calculations are p erformed on a MacBo ok Pro with an Apple M4 chip. T o ev aluate the effectiv eness of the TDVP-2 and the rank-adaptiv e BUG metho ds, we apply the al- gorithms to b oth a time-independent and a time-dep enden t Hamiltonian mo del. In the time-indep enden t case, we study the transverse Ising mo del. In the time-dependent case, w e consider a simplified mo del of a quan tum device sub ject to time-dep enden t con trol pulses. 6.1 Time-indep enden t Hamiltonians W e start by considering the one-dimensional transverse-field Ising mo del for N qubits, with Hamiltonian H “ ´ J N ´ 1 ÿ j “ 1 S z j S z j ` 1 ´ g N ÿ j “ 1 S x j , (6.1) whic h has a matrix pro duct op erator (MPO) representation with b ond dimension 3, see [30]. The op erators S z j and S x j are defined by , S z j “ I bp j ´ 1 q b S z b I bp N ´ j q , S x j “ I bp j ´ 1 q b S x b I bp N ´ j q , where S z “ 1 2 „ 1 0 0 ´ 1 ȷ , S x “ 1 2 „ 0 1 1 0 ȷ . Unless otherwise noted, all n umerical sim ulations in this section use a random product state as the initial state ψ p 0 q , the error in the states is measured in L 2 norm at the final time T “ 10 . 0, and the exact state is obtained via matrix exp onen tiation. Without truncation of the SVD ( ϵ “ 0), our implemen tations of the TD VP-2, T uc ker-BUG, and MPS-BUG algorithms hav e all b een v erified to b e second order accurate in the time-step. A ma jor adv antage of tensor decomp ositions is that adding sub-systems (qubits) to the comp osite quan- tum system do es not necessarily lead to an exponentially increasing computational cost. This is demon- strated in Figure 5, whic h shows run times and b ond dimensions for ev olving the Ising mo del to T “ 5 . 0, for increasing n umbers of sub-systems ( N ). Here, w e compare TDVP-2 and MPS-BUG, with and without transv erse field terms. Without transv erse field ( g “ 0), the Hamiltonian is diagonal in the computational basis, and en tanglemen t can only b e generated through phase accum ulation, and only if the initial state is in a sup erposition state. Ho wev er, these sim ulations start from the ground state, and the solution therefore ev olves as a pro duct state. As a result, the b ond dimensions are b “ 1 throughout the simulation, explaining the basically linear gro wth in run time as N increases, for b oth metho ds. In the case with transv erse field terms, g “ 0 . 5, superp ositions are generated for an y initial condition, leading to increasing en tanglemen t b y the coupling terms, S z j S z j ` 1 . In this case, the maximum b ond dimensions for TDVP-2 increases initially , but then plateaus at b “ 6 for most N . This b ehavior leads to an almost linear growth in run time for TD VP-2. F or MPS-BUG, it’s maxim um bond dimension is identical to TD VP-2 for N ď 26. F or larger n umbers of qubits, the bond dimension for MPS-BUG keeps increasing with N , leading to longer run times than TD VP-2. The reason for the increase in b ond dimension is curren tly unclear. Figure 6 (left panel) sho ws the state error at final time T “ 10 . 0, as function of the SVD threshold ϵ , for T uc k er-BUG, TD VP-2, and MPS-BUG, for a transverse Ising model with N “ 10 qubits. Here each color corresp onds to a differen t time-step. Compared to TDVP-2, twice as man y time-steps are needed b y 15 Figure 5: T ransv erse Ising mo del with J “ 1 and g “ t 0 , 0 . 5 u , starting from the ground state, in tegrated with TD VP-2 and MPS-BUG to time T “ 5. Left: run times as functions of the num b er of qubits. Right: Maxim um b ond dimensions. b oth BUG methods to get similar solution errors. This is b ecause TDVP-2 takes tw o sub-steps within eac h full time-step. As the n umber of time-steps increases, the SVD truncation error accumulates more rapidly . Tw o regimes can be iden tified. F or sufficiently small ϵ , the accuracy is alwa ys dominated by the O p δ 2 q time-stepping error. F or each n umber of time-steps , the colored squares and circles indicate the largest ϵ “ ϵ trunc where the solution error is dominated by the O p δ 2 q time-stepping error. Hence, using ϵ ă ϵ trunc only increases the computational effort, without impro ving the accuracy in the solution. F or fixed v alues of ϵ ą ϵ trunc , the error can increase when the num ber of time-steps increases ( δ decreases). F rom these results w e conclude that ϵ m ust decrease with δ to ac hieve a second order accurate solution. F or MPS-BUG, we also note an in termediate range for ϵ ą ϵ trunc , where the solution is approximately first order accurate. The maxim um n umber of en tries stored in the tensor decompositions during the time-ev olution of the Ising mo del are shown in Figure 6 (righ t panel). Tw o trends emerge. F or a fixed SVD threshold, increasing the n umber of time-steps reduces the b ond dimensions due to more frequent SVD truncations, leading to reduced storage. Ho wev er, for a fixed n um b er of time-steps, the storage requirements increase when the SVD threshold is reduced. This highlights the trade-off b et ween accuracy and memory usage inheren t to tensor-decomp osition metho ds. W e no w analyze the same exp eriment using the T uck er-BUG in tegrator, see Figure 6 (bottom ro w). As b efore, when the SVD threshold is sufficiently small, the state error is dominated by the time-stepping error and con verges to the same lev el as for TD VP-2 and MPS-BUG. W e also note that the error in the T uck er- BUG metho d remains small for a larger ϵ compared to TD VP-2 and MPS-BUG. How ever, the T uc ker-BUG metho d do es not show the same gradual reduction in storage as TDVP-2 and MPS-BUG. This b ehavior reflects the binary all-or-nothing truncation options for T uck er tensors when applied to tw o-lev el (qubit) sub-systems. T ensor decomposition metho ds can b e used to study the dynamics of observ ables, e.g., the magnetization at eac h site, m j “ x ψ | M j | ψ y , j P r 1 , N s , (6.2) where M j “ 2 I b j ´ 1 b S z b I b N ´ j . In Figure 7, we consider the magnetization in an Ising mo del with N “ 20 qubits that is evolv ed to final time T “ 15. The initial data has magnetization ` 1 at site N , and ´ 1 at all other sites. P anels (a)-(c) sho w the magnetization using TD VP-2 for different v alues of the SVD threshold ϵ . By visual insp ection, the magnetization do es not change b et w een ϵ “ 1 . 36 ¨ 10 ´ 5 (panel (a)) and ϵ “ 2 . 28 ¨ 10 ´ 4 (panel (b)). Ho wev er, the largest v alue of ϵ “ 8 . 89 ¨ 10 ´ 3 (panel (c)) give a v ery differen t result, indicating that the SVD truncation lead to a significan t error. In this case, the state has N tot « 10 6 elemen ts, making it intractable to calculate an exact solution via matrix exp onen tiation on the av ailable hardw are. T o illustrate ho w the ϵ -threshold in 16 Figure 6: Accuracy and storage requiremen ts as function of the SVD truncation parameter ( ϵ ) for TD VP-2 (top ro w), MPS-BUG (middle ro w), and T uc k er-BUG (bottom ro w). Here, J “ 1 and g “ t 0 . 5 , 1 . 0 u in the Ising mo del with N “ 10 qubits. Colors indicate different n umbers of time-steps, solid lines correspond to g “ 1 . 0 and dashed lines to g “ 0 . 5. Left column: State error at final time T “ 10 . 0. Righ t column: Maximum n um b er of en tries stored for eac h tensor decomp osition. 17 (a) TD VP-2, ϵ 1 (b) TD VP-2, ϵ 2 (c) TD VP-2, ϵ 3 (d) Magnetization difference at t “ 5 . 0. Figure 7: The magnetization observ able in a comp osite system with N “ 20 qubits ev olving under the Ising Hamiltonian with J “ 1 . 0 and g “ 0 . 5, for different SVD ϵ -thresholds. The initial state has magnetization m j “ ´ 1 for j P r 1 , 19 s and m 20 “ ` 1, and the system is ev olved to final time T = 15.0 using 256 time-steps. In the b ottom-righ t panel, we ev aluate the magnetization vector at time t “ 5 . 0 for ϵ k “ 10 ´ 3 ¨ 10 ´ k { 4 for k P r 1 , 20 s , and rep ort the norm of the difference }  m k ´  m k ´ 1 } as function of ϵ k . 18 Figure 8: Real and imaginary parts of the control pulse applied to qubit #1, in the case of N “ 6 qubits in the system. the SVD truncation influences the magnetization, we ev aluate the magnetization vector  m “ r m 1 , . . . m N s at time t “ 5 as function of ϵ . Starting with ϵ 0 “ 10 ´ 3 w e define ϵ k “ ϵ k ´ 1 ¨ 10 ´ 1 { 4 , for k P r 1 , 20 s . This giv es ϵ k “ ϵ 0 ¨ 10 ´ k { 4 and ϵ 20 “ 10 ´ 8 . W e denote the final magnetization, calculated with ϵ “ ϵ k , by  m k and then study the norm of differences in magnetization, ∆ k : “ }  m k ´  m k ´ 1 } , for k “ 1 , 2 , . . . , 20, see Figure 7(d). The slop e is « 1 . 25 indicating that the error in magnetization satisfies  m k ´  m exact “ O p ϵ p k q , with p “ 1 . 25. 6.2 Time-Dep enden t Hamiltonians Time-dep enden t Hamiltonians are common in quan tum control applications, where they are used to mo del quan tum dynamics due to external con trol pulses applied to the system. In the following we consider the Hamiltonian mo del for a linear chain of N coupled sup erconducting transmon qubits stated in (2.5)-(2.7). The dimension of the Hilbert (vector) space of eac h sub-system is d k “ 2, resulting in a state vector | ψ y P C 2 N . The MPO representation of this Hamiltonian has bond dimension equal to 4, for all sites. In the examples b elo w, w e tak e the frequency of rotation, ω d , to be the arithmetic mean of the transition frequencies t ω k u N k “ 1 . The time-dep endence of the control functions p k p t q and q k p t q is assumed to be a sum of harmonic carrier w av es with frequency Ω k , where the en velope function is parameterized b y B-spline w av elets. As a test case, we consider the state-to-state problem where the state of each qubit is transformed from the N -qubit ground state into the sup erposition state, | 0 y b N Ý Ñ ´ | 0 y`| 1 y ? 2 ¯ b N . (6.3) Due to the inheren t coupling betw een the qubits, it w ould b e c hallenging to analytically determine the con trol functions for the state-to-state transformation. Instead, w e numerically determine the control pulses in the con trol Hamiltonian (2.7) using the Quandary code [28]; for details of the n umerical approach, see [31, 32]. In this approac h, Quandary sim ultaneously optimizes for the con trol pulses in all qubits. F or example, in the case of N “ 6, the optimized con trol pulse for qubit #1 is shown in Figure 8. In the following, we test the numerical accuracy of the TDVP-2 and T uck er-BUG algorithms by compar- ing the state and infidelity with those obtained with Quandary . Due to the time dependence in the con trol Hamiltonian, the evolution of the quan tum system cannot b e ev aluated analytically . Instead, a reference solution is obtained using Quandary . 6.3 Coupled systems with optimized con trol pulses W e consider the case of short-range interaction in the Hamiltonian and set the coupling strength betw een qubits k and ℓ to b e J kℓ 2 π “ # 5 . 0 MHz , ℓ “ k ` 1 , 0 . 0 , otherwise . 19 Figure 9: Solution accuracy (left column) and storage requirements (righ t column) at the final time T “ 40 for systems with N “ 6 qubits (top row) and N “ 10 qubits (bottom ro w). Left column sho ws the norm of the differences b et w een the states from the tensor-decomp osition methods and Quandary , as function of the ϵ -threshold in the truncated SVD. Also shown are the infidelities from TD VP-2, T uck er-BUG, and Quandary . Right column sho ws the total num ber of en tries stored in the MPS and the T uc ker-tensor as function of ϵ . The dashed green-line indicates the storage needed for one state vector in Quandary . Here, we choose the transition frequencies to b e ω k { 2 π “ 4 . 64 ` 0 . 06 p k ´ 1 q GHz, for k P r 1 , N s . Because of the coupling betw een qubits, we use three carrier frequencies for each qubit (tw o for the first and last qubits in the chain). The carrier frequencies are Ω k “ t ω k ´ 1 , ω k , ω k ` 1 u ´ ω d . T o establish a reference solution, w e ev olve this system with Quandary to final time T “ 40 . 0 and used Ric hardson extrap olation to estimate the error in the final state vector. F or N “ 6 qubits, with 4152 time- steps, the norm of the error was « 2 . 97 ¨ 10 ´ 5 . F or N “ 10 qubits, with 10,480 time-steps, the norm of the error became « 7 . 82 ¨ 10 ´ 5 . Figure 9 (left column) shows the difference b et w een the states computed with tensor decomposition metho ds (TD VP-2, T uc ker-BUG) and Quandary , as function of the ϵ -threshold in the truncated SVD. W e note that the infidelit y computed with the TDVP-2 and T uc k er-BUG metho ds is at least one order of magnitude smaller than the norm of the state-differences. Also note that the fidelit y obtained b y both tensor decom positions agree with Quandary when ϵ ď 10 ´ 5 ( N “ 6) and ϵ ď 5 ¨ 10 ´ 6 ( N “ 10), resp ectiv ely . Figure 9 (right column) sho ws the total num b er of en tries stored in the tensor decomp ositions compared to the vector storage in Quandary . In this figure we note the all-or-nothing truncation prop ert y of the T uck er- BUG algorithm as the b ond dimension in the factor matrices drops from 2 to 1 at ϵ « 10 ´ 4 . Compared to Quandary , TDVP-2 uses less storage when ϵ ą 5 . 5 ¨ 10 ´ 5 ( N “ 6) and when ϵ ą 10 ´ 9 ( N “ 10), resp ectively . F or the N “ 10 case, w e conclude that TDVP-2 needs less storage than Quandary to calculate b oth the final state and the fidelit y with the same accuracy . 20 Figure 10: Comparing Quandary and TD VP-2 on a state-to-state transformation problem for differen t n um b ers of qubits with fixed SVD threshold, ϵ “ 10 ´ 6 . Left: run time. Right: infidelit y . 6.4 Larger systems of qubits with optimized control pulses Consider a system with N qubits in a linear chain where the transition frequencies oscillate with p eriod 8: r p k q “ p k ´ 1 q mod 4 , ω k “ # 5 . 18 ´ 0 . 06 r p k q , k “ 1 , 2 , 3 , 4 , 9 , 10 , 11 , 12 , . . . 5 . 18 ` 0 . 06 r p k q ´ 0 . 15 , k “ 5 , 6 , 7 , 8 , 13 , 14 , 15 , 16 , . . . Because it b ecomes difficult to numerically optimize the con trol pulses as the num b er of qubits b ecomes larger, w e here consider the case with J kℓ “ 0. This allows the control pulses for the state-to-state transformation (6.3) to b e optimized indep endently for each qubit. As a result, only 1 carrier frequency is needed p er qubit, Ω k “ ω k ´ ω d . Figure 10 sho ws a run time comparison b etw een TDVP-2 and Quandary for an increasing num b er of qubits. T o obtain a fair comparison with TDVP-2, which is implemented as a serial co de, w e also execute Quandary on a single core. The T uc ker-BUG integrator was omitted from these test b ecause the algorithm b ecomes very slow for N ě 10 qubits. As is shown in Figure 10, Quandary is faster than TDVP-2 for systems with 6-13 qubits. P ast 13 qubits, TD VP-2 b ecomes faster than Quandary . In Quandary , the exponential increase of the state-v ector leads to an exp onen tial increase in run time. With TDVP-2, the b ond dimension can be k ept small and indep enden t of N b ecause there is no coupling betw een the qubits ( J kℓ “ 0). This results in run times that only increase linearly in the n um b er of qubits. In Figure 10 (right), w e compare the infidelities from TD VP-2 and Quandary . W e conclude that both tec hniques result in almost identical results. 7 Conclusions In this pap er w e ha v e deriv ed and n umerically ev aluated tensor train and T uc ker tensor decomp osition metho ds for solving Schr¨ odinger’s equation in the time domain. F or time-in tegration of tensor trains, we review ed the classical TDVP and TDVP-2 algorithms, and also outlined the MPS-BUG algorithm. F or T uc ker decompositions, we sk etc hed the T uc ker-BUG algorithm. In n umerical experiments w e ha ve ev aluated the accuracy , memory usage and run times of the rank-adaptiv e TDVP-2, MPS-BUG and T uck er-BUG metho ds. The rank-adaptivity in all the methods is based the ϵ -truncated singular v alue decomp osition (SVD). Solution accuracy and run times ha ve b een compared with the conv en tional matrix-v ector form ulation implemen ted in the Quandary [31, 28] softw are. Numerical experiments based on composite quan tum systems consisting of t w o-level sub-systems (qubits) w ere presen ted for t w o Hamiltonian models: the time-indep enden t transv erse Ising mo del, and a mo del of transmon qubits sub ject to time-dependent control pulses. F or the tw o-lev el cases considered here, T uc ker tensors exhibit an inheren tly binary b eha vior in the ϵ -truncated SVD, where the factor matrices either hav e 21 full size (2 ˆ 2), or the truncated size (2 ˆ 1). As a result, each sub-systen is either treated as fully entangled, or completely decoupled, from the rest of the system. T uck er tensors are therefore of limited use for t wo-lev el systems, but may b e more effective for quantum systems with higher lo cal dimensions, where more gradual truncation is p ossible. F or the transverse Ising mo del with short-range interactions, we ha v e demonstrated that quan tum sim- ulations with tensor trains can b e tractable for systems with 100 qubits (or more), using a laptop computer. This is possible b ecause of the limited gro wth in en tanglement during the ev olution, leading to run times that gro w as lo w-order polynomials, instead of exp onen tially with con ven tional matrix-vector form ulations. F or mo dest num bers of sub-systems, where the accuracy can be ev aluated via exact solutions based on matrix exp onen tiation, we hav e analyzed the effect of the ϵ -threshold in the truncated SVD, and how it needs to b e related to the time-step. W e also studied how the ϵ -truncation influences the magnetization observ able in a 20-qubit system. When the state is entangled, the MPS-BUG metho d is observed to b e slo wer than TD VP-2. A part of this can b e attributed to the SVD factorizations during the truncation of the inflated b ond dimensions at the end of eac h time-step. Ho wev er, w e hav e also observ ed that the b ond dimension in MPS-BUG sometimes b ecomes larger than in TDVP-2, esp ecially when the n umber of coupled sub-systems is large. The reason for this b eha vior is curren tly po orly understo o d. The performance of the tensor decomp osition methods for time-dep enden t Hamiltonian models was ev aluated on a state-to-state transformation problem. Here, we considered a quantum system of N coupled transmon qubits sub ject to time-dependent con trol pulses. Reference solutions w ere generated with the Quandary code and used to estimate the accuracy of the tensor decomposition metho ds. W e studied ho w the ϵ -threshold in the truncated SVD affects the accuracy and storage requirements in the TD VP-2 and T uc ker-BUG algorithms. F or an ϵ that giv es the same accuracy as Quandary , the TD VP-2 algorithm uses less storage than Quandary , if N is sufficien tly large. F or the same ϵ , the T uck er-BUG algorithm ga ve similar accuracy as TDVP-2, but w as slow er and used significan tly more storage. Finally , we compared the run time p erformance of TD VP-2 and Quandary as the n umber of qubits increases. In this case, w e observ ed that TD VP-2 outp erforms Quandary when the n umber of qubits exceeds N « 13, without degrading the fidelit y of the state-to-state transformation. One natural extension of the current work is to consider qibits that are connected in a tw o-dimensional lattice. Another interesting extension is to apply tensor train metho ds to mo del quantum decoherence by solving Sc hr¨ odingers equation with a non-Hermitian Hamiltonian, e.g., within the Mon te-Carlo w av efunction metho d. Ac kno wledgmen ts W e gratefully ac knowledge financial support from the office of Adv anced Scien tific Computing Researc h (ASCR) within the U.S. Departmen t Of Energy , aw ard #SCW1886. The authors thank Dr. Dominik Sultz and Dr. Jonas Kusc h for helpful discussions on the BUG algorithm for tensor trains. This work was p erformed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Con tract DE-A C52-07NA27344. This is contribution LLNL-JRNL-2016933. A The singular-v alue decomp osition Consider a matrix M P C m ˆ n . The (compact) Singular V alue Decomposition (SVD) of M is given by M “ U Σ V : , U : U “ I r , V : V “ I r , Σ “ diag r σ 1 , σ 2 , ..., σ r s , r ď min p m, n q , (A.1) where U P C m ˆ r , V P C n ˆ r , Σ P R r ˆ r , and r is the rank of M . The elements of Σ are real and p ositive and are called the singular v alues of M . In the following we assume them to b e sorted in descending order, σ 1 ě σ 2 ě ¨ ¨ ¨ ě σ r . Sev eral algorithms for constructing lo w-rank appro ximate representation s of tensors are based on low- rank approximations of the SVD. Denote b y Ă M an approximation of the matrix M , formed by setting the smallest singular v alues to zero, r Σ “ diag r σ 1 , ..., σ k , 0 , ..., 0 s . (A.2) 22 Because the matrices U and V hav e orthonormal columns, the F rob enius norm difference can b e calculated exactly , } M ´ Ă M } F “ } U Σ V : ´ U r Σ V : } F “ } Σ ´ r Σ } F “ g f f e r ÿ i “ k ` 1 σ 2 i ď ϵ. (A.3) Here, k ě 1 is the largest index such that the inequalit y holds for a giv en ϵ ě 0. The matrix Ă M has k non-zero singular v alues and therefore has rank k . One can show that Ă M is the rank- k matrix that is the closest to M in F ro ebenius norm, th us being the optimal approximation of M . A.1 Order-4 T ensor SVD In TDVP-2 a SVD is p erformed on an order-4 tensor to decompose it in to tw o order-3 tensors and mo ve the orthogonality center to the right. This is done by reshaping the tensor into a matrix, p erforming a SVD on that matrix, follo wed by inv erse reshape op erations. F or example, let T P C b i ˆ d i ˆ d i ` 1 ˆ b i ` 2 represen t the merged core for sites i ` 1 and i ` 2. Reshap e T in to a matrix, ˆ T : “ reshap e p T q P C b i d i ˆ d i ` 1 b n ` 2 . P erform a SVD on ˆ T , ˆ U ˆ S ˆ V : : “ ˆ T , ˆ U P C b i d i ˆ s i , ˆ S P C s i ˆ s i , ˆ V : P C s i ˆ d i ` 1 b i ` 2 . (A.4) Here w e hav e s i “ min t b i d i , d i ` 1 b i ` 2 u , and w e set the in termediate b ond dimension to b i ` 1 “ s i . The in v erse reshap e op erations is then p erformed on the matrix ˆ U and on the matrix ˆ R : “ ˆ S ˆ V : , where ˆ R P C s i ˆ d i ` 1 b i ` 2 , A r i ` 1 s : “ reshap e ´ 1 p ˆ U q P C b i ˆ d i ˆ b i ` 1 , M r i ` 2 s : “ reshap e ´ 1 p ˆ R q P C b i ` 1 ˆ d i ` 1 ˆ b i ` 2 (A.5) Because ˆ U : ˆ U “ I s i , the tensor A r i ` 1 s is left-normalized, and the orthogonality cen ter has mov ed to site i ` 2. By contraction ov er the joint index b i ` 1 w e can verify T “ A r i ` 1 s ¨ M r i ` 2 s , if the SVD is performed without truncation. In the case of ϵ -truncation in the SVD, the bond dimension is reduced to b i ` 1 ă s i , where b i ` 1 ě 1 is set as the largest index that satisfies the inequality (A.3) for a giv en ϵ . The discarded singular v alues lead to the approximation r T : “ A r i ` 1 s ¨ M r i ` 2 s « T , where } T ´ r T } F “ } ˆ U ˆ Σ ˆ V : ´ ˆ U r Σ ˆ V : } F “ } ˆ Σ ´ r Σ } F ď ϵ. (A.6) Here, r Σ is the appro ximate singular v alue matrix where the b i ` 1 largest singular v alues of ˆ Σ are retained. References [1] T omisla v Begu ˇ si ´ c, Johnnie Gra y , and Garnet Kin-Lic Chan. F ast and conv erged classical simulations of evidence for the utilit y of quan tum computing before fault tolerance. Scienc e A dvanc es , 10(3):eadk4321, 2024. [2] Huanc hen Zhai, Chenghan Li, Xing Zhang, Zhendong Li, Seungho on Lee, and Garnet Kin-Lic Chan. Classical solution of the F eMo-cofactor mo del to chemical accuracy and its implications, 2026. arXiv- 2601.04621. [3] A.J. Daley , C. Kollath, U. Sc hollw¨ ock, and G. Vidal. Time-dependent density-matrix renormalization- group using adaptive effective Hilb ert spaces. J. Stat. Me ch. , 2004. [4] I. V. Oseledets. T ensor-train decomp osition. SIAM Journal on Scientific Computing , 33(5):2295–2317, 2011. 23 [5] Jutho Haegeman, J. Ignacio Cirac, T obias J. Osborne, Iztok Pi ˇ zorn, Henri V ersc helde, and F rank V erstraete. Time-dep enden t v ariational principle for quantum lattices. Phys. R ev. L ett. , 107:070601, Aug 2011. [6] Jutho Haegeman, Christian Lubic h, Iv an Oseledets, Bart V andereyc ken, and F rank V erstraete. Unifying time ev olution and optimization with matrix pro duct states. Phys. R ev. B , 94:165116, Oct 2016. [7] Sebastian Paec k el, Thomas K¨ ohler, Andreas Swobo da, Salv atore R. Manmana, Ulrich Schollw¨ oc k, and Claudius Hubig. Time-evolution metho ds for matrix-product states. A nnals of Physics , 411:167998, 2019. [8] F.A.Y.N. Schr¨ oder, D.H.P . T urban, A.J. Musser, et al. T ensor netw ork simulation of m ulti- en vironmental op en quantum dynamics via mac hine learning and entanglemen t renormalisation. Nat. Commun. , 10:1062, 2019. [9] Dominik Sulz. p ersonal comm unication. [10] Gianluca Ceruti and Christian Lubic h. An uncon ven tional robust in tegrator for dynamical lo w-rank appro ximation. BIT Numeric al Mathematics , 62(1):23–44, Mar 2022. [11] Gianluca Ceruti, Christian Lubic h, and Hanna W alac h. Time integration of tree tensor net works. SIAM Journal on Numeric al Analysis , 59(1):289–313, 2021. [12] Gianluca Ceruti, Daniel Kressner, and Dominik Sulz. Lo w-rank tree tensor net w ork op erators for long- range pairwise interactions. SIAM Journal on Scientific Computing , 47(4):A2248–A2271, 2025. [13] Jonas Kusc h. Second-order robust parallel in tegrators for dynamical lo w-rank appro ximation. BIT Numeric al Mathematics , 65(3):31, Jun 2025. [14] Daniel Appel¨ o and Yingda Cheng. Kraus is king: High-order completely p ositiv e and trace preserv- ing (CPTP) lo w rank method for the Lindblad master equation. Journal of Computational Physics , 534:114036, 2025. [15] Spencer Lee and Daniel Appel¨ o. High-order Hermite optimization: F ast and exact gradien t computation in open-lo op quan tum optimal con trol using a discrete adjoin t approac h. Journal of Computational Physics , 552:114697, 2026. [16] N. Anders Petersson, Stefanie G ¨ un ther, and Seung Whan Ch ung. A time-parallel multiple-shooting metho d for large-scale quantum optimal control. Journal of Computational Physics , 524:113712, 2025. [17] Grey Ballard and T amara G. Kolda. T ensor De c omp ositions for Data Scienc e . Cam bridge Universit y Press, 2025. [18] C. Hubig, I. P . McCullo c h, and U. Schollw¨ oc k. Ge neric construction of efficien t matrix pro duct opera- tors. Phys. R ev. B , 95:035129, Jan 2017. [19] P . A. M. Dirac. Note on exchange phenomena in the Thomas atom. Mathematic al Pr o c e e dings of the Cambridge Philosophic al So ciety , 26(3):376–385, 1930. [20] J. F renkel. Wave me chanics: advanc e d gener al the ory. Clarendon Press Oxford, 1964. [21] A.D. McLac hlan. A v ariational solution of the time-dep enden t sc hro dinger equation. Mole cular Physics , 8(1):39–44, 1964. [22] P .H. Kramer and M. Saraceno. Ge ometry of the Time-Dep endent V ariational Principle in Quantum Me chanics . Springer, 1981. [23] J. Bro ec khov e, L. Lathou wers, E. Kestelo ot, and P . V an Leuven. On the equiv alence of time-dep enden t v ariational principles. Chemic al Physics L etters , 149(5):547–550, 1988. [24] Christian Lubich, Iv an V. Oseledets, and Bart V andereyck en. Time integration of tensor trains. SIAM Journal on Numeric al Analysis , 53(2):917–941, 2015. 24 [25] Gianluca Ceruti, Jonas Kusc h, Christian Lubic h, and Dominik Sulz. A parallel basis up date and Galerkin integrator for tree tensor net works. SIAM Journal on Scientific Computing , 48(1):A27–A48, 2026. [26] Gianluca Ceruti, Jonas Kusc h, and Christian Lubic h. A rank-adaptive robust integrator for dynamical lo w-rank appro ximation. BIT , 62(4):1149–1174, Decem b er 2022. [27] Klaus Mølmer, Yv an Castin, and Jean Dalibard. Monte Carlo w av e-function method in quan tum optics. J. Opt. So c. Am. B , 10(3):524–538, Mar 1993. [28] N. Anders P etersson Stefanie G¨ un ther. Quandary: Optimal control for open and closed quantum systems. https://github.com/LLNL/quandary/ , 2025. [29] Matthew Fishman, Stev en R. White, and E. Miles Stoudenmire. The iT ensor Soft ware Library for T ensor Net w ork Calculations. SciPost Phys. Co deb ases , page 4, 2022. [30] Maarten V an Damme, Jutho Haegeman, Ian McCullo c h, and Laurens V anderstraeten. Efficien t higher- order matrix pro duct op erators for time evolution. SciPost Phys. , 17:135, 2024. [31] Stefanie Gunther, N. Anders Petersson, and Jonathan L. DuBois. Quandary: An op en-source C++ pac k age for high-p erformance optimal con trol of open quantum systems . In 2021 IEEE/ACM Se c ond International Workshop on Quantum Computing Softwar e (QCS) , pages 88–98, Los Alamitos, CA, USA, No vem b er 2021. IEEE Computer So ciet y . [32] N. Anders P etersson and F ortino Garcia. Optimal con trol of closed quantum systems via B-splines with carrier w a ves. SIAM Journal on Scientific Computing , 44(6):A3592–A3616, 2022. 25