PaperlayerView exact revision on arXiv

Source-linked AI summary

Time-evolution methods for matrix-product states

Sebastian Paeckel, Thomas Köhler, Andreas Swoboda, Salvatore R. Manmana, Ulrich Schollwöck, Claudius Hubig

arXiv:1901.05824v3cond-mat.str-elcond-mat.stat-mechquant-ph

TL;DR

Time evolution of MPS requires combining accurate propagation with efficient truncation, but existing approaches have different strengths and weaknesses. This review compares operator-based, Krylov, and TDVP methods across representative correlated-system problems, finding strong performance for 2TDVP while identifying method- and timestep-dependent trade-offs.

  • Problem

    MPS time evolution must combine accurate time propagation with efficient truncation, while a direct comparison of available methods has been missing.

  • Method

    The review explains and compares TEBD, MPO W^II, global and local Krylov, and one- and two-site TDVP methods across four representative dynamics problems.

  • Results

    2TDVP provided the best numerical data and efficient calculations for dynamical spin structure factors, while TEBD2 offered a runtime–accuracy trade-off and other methods showed timestep-dependent errors.

  • Takeaways & Limitations

    For longer times, 2TDVP or local Krylov appeared most promising, whereas for short times all tested methods worked reasonably well.

  • Takeaways & Limitations

    Local Krylov time evolution includes errors from local TDSE solution, SVD truncation, and sequential Lie-Trotter evolution.

Abstract

from arXiv · show

Matrix-product states have become the de facto standard for the representation of one-dimensional quantum many body states. During the last few years, numerous new methods have been introduced to evaluate the time evolution of a matrix-product state. Here, we will review and summarize the recent work on this topic as applied to finite quantum systems. We will explain and compare the different methods available to construct a time-evolved matrix-product state, namely the time-evolving block decimation, the MPO $W^\mathrm{II}$ method, the global Krylov method, the local Krylov method and the one- and two-site time-dependent variational principle. We will also apply these methods to four different representative examples of current problem settings in condensed matter physics.

1. Introduction

This review addresses how to evolve finite-system matrix-product states accurately and efficiently, comparing methods that approximate the propagator directly or act on the state. It evaluates these approaches across four representative correlated-system dynamics problems.

  • Motivation: MPS time evolution combines efficient Hilbert-space truncation with accurate numerical integration of the time-dependent Schrödinger equation.The review focuses on finite-dimensional quantum states evolved in real or imaginary time.
  • Method classes: Two broad strategies either approximate the action of the time-evolution operator on the state or construct an efficient approximation to the operator itself.The global Krylov method belongs to the first strategy, while Suzuki-Trotter and MPO methods belong to the second.
  • Operator-based methods: Suzuki-Trotter decomposition expresses sufficiently small time steps through smaller matrix exponentials and can be applied directly to MPS time evolution.MPO approximations extend this idea by exploiting matrix-product structure for efficient exponentiation.
  • Scope and contribution: The review fills a missing direct comparison by examining Suzuki-Trotter, MPO W^I,II, TDVP, and Krylov-subspace approaches.The methods are tested on four representative problems involving real- and imaginary-time evolution, correlators, critical and gapped systems, and two-dimensional settings.

2. Matrix-product states and operators

MPS represent one-dimensional quantum states as products of site tensors linked by virtual indices, while tensor notation describes their contractions and index structure. Their efficiency is tied to bond dimensions and weak entanglement.

  • Motivation: The exponentially growing Hilbert space restricts exact diagonalization to roughly 40–50 sites for spin systems, motivating MPS-based methods.MPS methods address this restriction particularly for one-dimensional systems.
  • Tensor notation: A tensor is a multidimensional collection of numbers, with each index corresponding to a dimension and each graphical leg representing one index.Tensor contractions are represented graphically by connecting shared legs.
  • Matrix-product states: An MPS expresses the coefficient tensor of a lattice quantum state as a product of L rank-3 site tensors connected by virtual indices.The end indices m0 and mL are one-dimensional dummy indices, while σj labels local physical states.
  • Matrix-product states: Different virtual bond dimensions determine which quantum states can be represented exactly, and unrestricted growth can represent any state as an MPS.MPS are efficient representations for one-dimensional weakly entangled states.

2.3. Matrix-product operators (MPO)

Matrix-product operators represent operators as contractions of rank-4 site tensors, with virtual bonds encoding operator structure. Canonical tensor gauges provide normalized forms that simplify contractions and stabilize calculations.

  • MPO structure: An MPO represents an operator as a contraction of L rank-4 tensors with physical domain and image indices.Its virtual bonds have dimensions wj, with maximal bond dimension w.
  • MPO construction: For operators built from local terms, MPO construction separates terms within each partition from terms connecting the two partitions at a bond.Operator-valued matrices define recursion relations that can be implemented using finite-state machines.
  • Canonical forms: Gauge freedom allows resolutions of the identity to be inserted between MPS tensors without changing the represented state.Left- and right-normalized tensors contract with their adjoints to yield identities, enabling canonical forms.

2.5. Normalizing an MPS

MPS normalization is obtained by sequential QR decompositions that move transfer factors between neighboring tensors. Sweeping from either edge produces a fully left- or right-normalized state.

  • QR normalization: QR decomposition reshapes a site tensor into a matrix, producing a left-normalized tensor and a transfer matrix for the neighboring site.The transfer matrix Rj is multiplied into Mj+1 during a left-to-right sweep.
  • QR normalization: Right normalization instead groups the right and physical legs as matrix rows and the left tensor leg as columns.This produces a right-normalized tensor Bj.
  • Canonical sweeps: A complete left-normalized state results from sweeping left to right, while a complete right-normalized state results from sweeping right to left.The normalization direction determines the canonical form obtained.

2.6. Truncating an MPS

MPS truncation seeks a smaller-bond-dimension state that approximates the original, balancing local optimality against global quality. SVD truncation is efficient, while iterative variational sweeps can improve the approximation but depend on initialization and may require two-site updates.

  • Time evolution increases entanglement and therefore generally requires larger MPS bond dimensions, making accurate approximation methods essential.
  • Direct truncation via SVD: Sequential bondwise SVDs are locally optimal but may fail to produce a globally optimal approximation when truncation errors are large.
  • Direct truncation via SVD: SVD truncation keeps the m′ largest singular values, with the discarded weight quantifying the approximation error.
  • Direct truncation via SVD: A target discarded weight such as 10^-10 can correspond to an error of approximately 10^-5 · L.
  • Iterative variational truncation: Variational truncation optimizes site tensors through repeated sweeps, but convergence depends strongly on the initial guess and can stall at a locally optimal state.
  • Iterative variational truncation: Two-site optimization permits changing the bond dimension and redistributing quantum-number sectors, unlike a single-site update.

2.7. Finite temperatures

Finite-temperature mixed states are represented within the MPS framework by purifying them with auxiliary degrees of freedom. Imaginary-time evolution prepares the thermal state, while real-time evolution can act only on the physical space without changing the density matrix.

  • MPS represent pure states by default, so finite-temperature mixed states require purifications or minimally entangled typical thermal states.
  • At infinite temperature, the grand-canonical state can be initialized as a product of maximally entangled physical-auxiliary pairs.
  • Finite-temperature states are prepared by evolving the infinite-temperature purification along imaginary time over β/2, yielding ρ ∝ e^-βH after tracing out auxiliaries.
  • Real-time evolution can act on the physical degrees of freedom while leaving the auxiliary degrees untouched.
  • Auxiliary-space unitaries leave the density matrix invariant and may reduce purification entanglement, but their computational benefit depends on the system.

2.8. Application of an MPO to an MPS

Applying an MPO directly to an MPS multiplies bond dimensions, so practical methods combine operator application with truncation or variational compression. Direct, variational, and zip-up approaches trade simplicity, accuracy, and computational cost differently.

  • Direct MPO-MPS application produces an MPS with bond dimension m′ = m · w, which usually exceeds what is needed for an efficient representation.
  • Variational application: The variational update uses a mixed-canonical guess state and recursively constructed boundary tensors Lj−1 and Rj+1.
  • Subsequent direct SVD truncation costs m^3w^3d per site, whereas variational compression costs m′^2mwd per site.
  • Variational application: Variational MPO application minimizes the distance between a guess MPS and the MPO-applied source state through local tensor updates.
  • Zip-up method: Zip-up truncates during contraction under the assumption that the MPO only slightly disrupts the MPS canonical form.
  • Zip-up method: The zip-up method has leading SVD cost O(m^3σw), linear in the MPO bond dimension w, and may be followed by variational sweeps for improved accuracy.

2.9. Expectation values

MPS expectation values are evaluated by contracting the tensor network for ⟨φ|Ô|ψ⟩ rather than treating the operator application and overlap as separate dense operations. Iterative left-to-right or right-to-left contractions exploit reusable boundary tensors.

  • In dense linear algebra, evaluating ⟨φ|Ô|ψ⟩ involves applying Ô to |ψ⟩ and then computing the overlap with ⟨φ|.
  • The MPO-MPS expectation-value network can be contracted from left to right using Lj or from right to left using Rj.
  • Concurrent evaluation of left and right boundary contractions provides a simple form of two-fold parallelization.

3. Overview of time-evolution methods

The review compares five MPS time-evolution methods, emphasizing that their differing strengths and weaknesses make problem-specific method selection necessary. It also discusses practical techniques and tests the methods on representative settings.

  • Five methods are reviewed: TEBD, MPO W I,II, global Krylov, local Krylov, and one- and two-site TDVP.Each has different strengths and weaknesses, so the most promising approach depends on the individual problem.
  • TEBD and MPO methods approximate the time-evolution operator before repeatedly applying it to the MPS.These methods construct an approximation to U(δ) and use it to generate successive time-evolved states.
  • Krylov methods directly approximate U(δ)|ψ⟩ without explicitly constructing U(δ) in the full Hilbert space.They offer precise short-step evolution, while global Krylov may require MPS representations of highly entangled Krylov vectors.
  • The local Krylov method avoids representing highly entangled global Krylov vectors by solving the TDSE locally on lattice-site pairs.Its trade-off is that it cannot evaluate observables at arbitrarily small intermediate time steps and may incur uncontrolled projection error.
  • The review also covers techniques for longer-time evolution, adaptive time-step selection, and other practical tricks, alongside comparisons on prototypical examples.

4. Approximations to ˆU(δ)

This section introduces two approaches that approximate the time-evolution operator and apply the resulting operator to an MPS to produce the next time-evolved state. It notes that more generic constructions remain possible.

  • The two approaches approximate U(δ), which is applied to |ψ(t)⟩ to obtain |ψ(t + δ)⟩.The same procedure produces subsequent states at later time steps.

4.1. Time-evolving block decimation (TEBD) or Trotter decomposition

TEBD constructs approximate time-evolution operators through Trotter-Suzuki decompositions and applies them as MPOs, with accuracy controlled by time-step order and MPS truncation. Higher-order schemes reduce step counts but require more exponentials per step, while swap gates extend the approach to long-range terms.

  • TEBD decomposes the Hamiltonian into internally commuting parts whose exponentials form an approximate time-evolution operator.For nearest-neighbor terms, even and odd bond contributions can be exponentiated efficiently as MPOs.
  • Errors: O(δ) and O(δ2) are the total-interval errors for first- and second-order TEBD, respectively, when N = T/δ steps cover a fixed interval.The per-step errors are O(δ2) and O(δ3), respectively.
  • Second-order TEBD symmetrizes the first-order decomposition, while fourth-order TEBD further increases the decomposition order.The second-order operator has second-order error per time step; fourth-order construction requires many more terms.
  • Errors: TEBD time-step error preserves unitarity, but truncation error affects both unitarity and conserved quantities.Truncation convergence can be estimated by increasing the MPS bond dimension and monitoring discarded weight.
  • Errors: 10^-8 requires 10^8, 10^4, and 10^2 steps for TEBD1, TEBD2, and TEBD4, respectively.TEBD4 performs approximately five times as much work per step while reducing the number of steps by a factor of 100 relative to TEBD2.
  • Long-range terms: Swap gates bring long-range interacting sites together for exponentiation and can have nearly all overhead removed through suitable ordering.

4.2. The MPO W I,II method

The MPO W I,II method constructs efficient MPO approximations to the time-evolution operator by exploiting the factorization and local structure of Hamiltonian MPOs. W I is efficient but restricted, while W II includes overlapping interaction terms with higher-order accuracy.

  • MPO W I,II approximations exploit operator factorization and can handle long-ranged interactions, making them suitable for broader geometries.The construction is designed so the error per site is independent of system size.
  • Errors: The naive W I approximation has O(Lδ2) total error from O(L) missed overlapping-term combinations, giving constant error per site.A local version is introduced because another Euler-step approximation can become more unstable as system size increases.
  • W I: W I uses the finite-state-machine structure of the Hamiltonian MPO and has MPO bond dimension w−1, enabling efficient numerical application.
  • W I: W I is restricted to non-overlapping local operator terms and fails to reproduce even purely on-site time evolution.This limitation motivates allowing operator strings that overlap on one site.
  • W II: W II retains terms with overlapping support while discarding terms that overlap at more than one site.Single-site terms are treated exactly to arbitrary powers, and the resulting approximation has O(δ3) error that does not affect the second-order approximation of W II.
  • W II: The W II site tensors are constructed from the MPO recursion for H by factorizing the exponential and evaluating operator-valued matrix elements.The derivation uses auxiliary degrees of freedom, complex Gaussian integrals, and coherent-state path integrals to obtain discrete MPO bond indices.

5. The global Krylov method

The global Krylov method directly approximates the action of the time-evolution operator by constructing an orthonormal Krylov subspace from repeated Hamiltonian applications. Its errors can be made very small, but maintaining orthogonality under MPS truncation can be costly and unstable.

  • Krylov-subspace construction: The method constructs a Krylov subspace from the initial state and successive applications of the Hamiltonian, orthonormalizing the resulting vectors.The first vector is the normalized initial state; subsequent vectors are generated recursively.
  • Projected evolution: The projected Hamiltonian is represented by an N × N effective matrix, which is exponentiated before mapping the result back to the Krylov basis.For Hermitian Hamiltonians, the effective matrix is ideally tridiagonal, enabling efficient exponentiation.
  • Observable evaluation: Expectation values can be evaluated without explicitly constructing the time-evolved MPS by combining projected observables with the Krylov coefficient vector.The observable can also be evaluated at intermediate times within the time step.
  • Error control: For L = 100 sites and δ = 0.1, the condition Wδ ≤ N is satisfied at approximately N ≈ 3 Krylov vectors.Here W denotes the Hamiltonian’s spectral range, roughly of the same order as system size.
  • Error control: The inherent Krylov error can often reach O(10^-10) or smaller, while MPS truncation error can also be monitored through the discarded weight.These small errors are attainable at finite time-step size, although at relatively large numerical cost.
  • Loss of orthogonality: Finite-precision arithmetic and MPS truncation can destroy basis orthogonality, and Gram-Schmidt re-orthogonalization may introduce further truncation errors.The resulting degradation can require successive orthogonality checks and re-orthogonalization.
  • Dynamic step sizing: A Krylov subspace computed for one time step can be recycled for interpolation or extrapolation at another step length.Dynamic step sizing is presented as a major flexibility of the method.

6. MPS-local methods

MPS-local methods evolve reduced one- or two-site problems within sweep-based local bases, then update the global MPS representation. They trade local solver and truncation errors against adaptive bond dimensions and conservation properties.

  • Local Krylov methods: Local Krylov methods construct reduced one- or two-site Schrödinger problems using left and right environment bases.The effective Hamiltonian and state are formed in smaller local spaces, where exponential evolution can be solved accurately.
  • Local Krylov methods: A sweep incrementally updates basis transformations and site tensors so the final state is represented in a basis optimized for t + δ.The local Krylov and TDVP schemes diverge in how they obtain the old state in the updated basis.
  • Two-site schemes: Two-site evolution uses an SVD after forward evolution, allowing the MPS bond dimension to adapt to entanglement growth.This flexibility distinguishes the two-site scheme from fixed-bond-dimension single-site evolution.
  • TDVP: Exact local Schrödinger solves preserve the energy and norm of the state, and TDVP achieves this conservation in practice.Other global observables may also be conserved under suitable conditions.
  • Errors: The standard two-site local Krylov method has four error sources, including local-solver error, SVD truncation error, and sequential Lie-Trotter error.The local-solver error can be made very small with a short Krylov basis, while truncation error is measurable during the calculation.
  • Errors: TDVP introduces projection error from restricting the TDSE to a limited-bond-dimension MPS manifold, while its symmetric integrator has O(δ3) per-step error.The projection error is zero at maximal bond dimension, and 2TDVP avoids it for nearest-neighbor Hamiltonians.
  • TDVP: Energy and norm are exactly conserved in 1TDVP and affected only by truncation error in 2TDVP.This conservation is useful for long-time hydrodynamic observables, although using only 1TDVP requires care.

7. Additional tricks

The review describes implementation tricks that extend accessible times, reduce time-step error, or accelerate post-processing. These gains come with trade-offs including observable-contraction cost, loss of unitarity, and numerical instability in prediction.

  • Combining pictures: Combining Schrödinger- and Heisenberg-picture evolution reaches times t1 + t2 with MPS and MPO bond dimensions typical of t1 and t2.The limiting operation then becomes evaluating observables between large-bond-dimension tensor networks rather than time evolution itself.
  • Complex time steps: Complex intermediate time steps raise a first-order propagator to third-order error per step and second-order overall.The construction requires multiple evolution operators, with cost growing linearly in their number.
  • Complex time steps: Complex-step schemes lose unitarity at each individual time step, which can be disadvantageous.For purely imaginary evolution with real Hamiltonians, avoiding complex arithmetic can reduce memory use by 50% and speed matrix multiplications approximately four-fold, but cannot reduce time-step error.
  • Linear prediction: Linear prediction assumes Green’s functions are composed of exponentially decaying and oscillating contributions, then fits recurrence coefficients to extend finite time series.The workflow discards an initial interval, fits on a selected interval, and verifies predictions against additional calculated data before extension.
  • Linear prediction: Linear-prediction fits can fail when the coefficient matrix R is singular, requiring a shift or a reduced parameter count.Choosing a suitable nonsingular parameter count may itself be difficult.
  • OTOCs: Direct OTOC evolution in physical degrees of freedom requires O(N^2) time steps at t = Nδ, motivating schemes that distribute the evolutions more evenly.The stated goal is linear scaling of effort in time.

8. Examples

Across representative benchmarks, 2TDVP generally provides the most accurate and stable time evolution, while other methods involve trade-offs between step size, runtime, conservation, and spectral fidelity. Method performance remains problem-dependent, with limitations from truncation, projection, and accumulated errors becoming important at longer times.

  • Dynamical structure factor: Outside the light cone, 2TDVP and global Krylov produce small homogeneous correlations, whereas local Krylov generates large contributions and spectral leakage.TEBD2 reproduces the light cone well but retains slow dynamics and spectral weight near ω = 0.
  • Dynamical structure factor: Only 2TDVP gives the correct spectral function at δ = 0.1, while global Krylov joins it at δ = 0.01.All methods locate the primary peak near ω ≈0.5, but several produce shifts or unphysical low-frequency spectral weight.
  • Dynamical structure factor: 2TDVP provides the best dynamical spin-structure-factor data and remains stable at the larger time step δ = 0.1.It is the only method reported to generate a stable evolution at δ = 0.1 in this benchmark.

9. Future developments

The paper identifies three directions for improving real-time MPS evolution: easier MPO construction, less entangled Krylov representations, and energy-conserving methods for hydrodynamics.

  • Future directions: The paper argues that direct real-time observables will continue to have important applications alongside alternative methods for excitation spectra and Green’s functions.It frames these applications as part of the future development agenda for real-time evolution methods.
  • Future directions for hydrodynamics: Local Krylov and TDVP methods may be promising for hydrodynamics, while 1TDVP enforces complete energy conservation but remains limited by bond dimension.The paper also proposes combining 2TDVP with a truncation scheme that preserves energy conservation.
  • Easier and more accurate construction of U(δ): Future MPO methods should simplify and improve the accuracy of constructing U(δ), whose current TEBD2 and MPO W I,II forms have relatively large time-step errors.The paper suggests constructions based on (1 − iδ′H)^N with δ′N = δ and large N, combined with MPO compression.
  • “Better” Krylov vectors: The global Krylov method’s main drawback is rapid entanglement growth in its Krylov vectors, motivating less-entangled MPS-representable subspaces or basis transformations.Such developments could retain the method’s flexible time-step choice and usefulness when Hamiltonian decompositions are not straightforward.
Loading 1901.05824v3…