Efficient time-evolution of matrix product states using average Hamiltonians

Technical Summary: Efficient Time-Evolution of Matrix Product States Using Average Hamiltonians

Problem Statement
Simulating the dynamics of quantum many-body systems (QMBS) governed by time-dependent Hamiltonians remains a significant challenge due to the exponential growth of the Hilbert space. While tensor networks, particularly Matrix Product States (MPS), offer an efficient framework for one-dimensional systems, standard time-evolution algorithms (such as Time-Evolving Block Decimation, TEBD, and MPO $W_{II}$) are primarily optimized for time-independent Hamiltonians. When applied to time-dependent scenarios, standard implementations often rely on the "Riemann stepper," which substitutes the instantaneous Hamiltonian at the beginning of a time interval ($H(0)$) into the evolution operator. This approach degrades the accuracy of high-order integrators, reducing the convergence rate of the evolution operator from third-order (for a single step) to second-order for the time-dependent case, thereby limiting simulation precision for a given step size.

Methodology
The authors propose a modification to existing MPS time-steppers to restore higher-order accuracy for time-dependent Hamiltonians without altering the fundamental architecture of tensor network libraries. The core methodology involves replacing the instantaneous Hamiltonian used within a single-step propagator with a carefully constructed "average Hamiltonian" derived from numerical integration theory.

Specifically, the authors substitute the standard Riemann approximation with a Hamiltonian calculated using Simpson's rule:
$$H_{\text{Simp}}(\Delta t) \equiv \frac{1}{6} \left( H(0) + 4H(\Delta t/2) + H(\Delta t) \right)$$
This average Hamiltonian is then fed into standard second-order TEBD or MPO $W_{II}$ algorithms. Theoretical analysis demonstrates that while the Riemann stepper ($H(0)$) yields an evolution operator with an error bound of $O(\Delta t^2)$, the Simpson stepper retains the native third-order local accuracy ($O(\Delta t^3)$) of the underlying integrators. The proof relies on bounding the difference between the exact time-evolution operator and the approximate operator using the Hölder-$\infty$ norm, showing that the Simpson approximation reduces the integration error term significantly compared to the Riemann approximation.

Key Contributions and Results
The proposed method was validated through numerical simulations of a chain of spin-1 sites representing nitrogen-vacancy (NV) color centers in diamond, driven by microwave pulses. The system was modeled using a Hamiltonian comprising drift, interaction, and control terms, with parameters derived from experimental data.

1. Accuracy Improvement: The simulations compared the "Simpson stepper" against the standard "Riemann stepper" across varying system sizes (3 to 8 NV centers) and integration step counts.

   • For small systems (e.g., 3 NV centers), the Simpson stepper reduced the average error by a factor of approximately 1000 for moderate step sizes compared to the Riemann stepper.
   • For larger systems (up to 8 NV centers), the method consistently yielded lower propagation errors, with improvements in average error ranging from factors of 50 to over 4000 depending on the step count and system size.
   • The error ratio ($E_{\text{Riemann}} / E_{\text{Simpson}}$) remained high across all tested configurations, confirming the robustness of the method.

2. Computational Cost: The increased accuracy comes with a modest increase in runtime. The Simpson stepper requires three Hamiltonian evaluations per time step (at $t=0$, $t=\Delta t/2$, and $t=\Delta t$) compared to one for the Riemann stepper. The authors note that this runtime increase is slight and attributed solely to the additional function evaluations, making the trade-off highly favorable for accuracy.

3. Scalability: The method demonstrated steady performance improvements as the system size increased, maintaining its advantage over the standard approach without introducing new truncation errors, provided the bond dimension was sufficiently large to keep truncation errors negligible.

Significance and Limitations
The paper claims that this approach provides a simple, yet highly efficient, augmentation to existing MPS solvers. By leveraging classical numerical integration ideas (Simpson's rule), the method restores the higher-order convergence properties of integrators that are otherwise lost when simulating time-dependent dynamics. This is particularly relevant for applications in quantum control, spectroscopy, and Floquet engineering where time-dependent driving fields are ubiquitous.

The authors explicitly note several limitations and practical considerations:

• Smoothness Requirement: The method relies on the existence and boundedness of the first four derivatives of the Hamiltonian $H(t)$ within each interval.
• Hamiltonian Construction Cost: If the time-dependent Hamiltonian is computationally expensive to construct, the requirement for three evaluations per step may significantly increase computational costs.
• Entanglement Limits: While the method improves integration error, it does not eliminate limitations imposed by entanglement growth in large-scale many-body systems, which are governed by truncation errors in the MPS representation.

The authors suggest that the Simpson averaging concept is adaptable and could be combined with higher-order quadrature rules (e.g., Gauss–Legendre) or multi-stage composition methods to construct even higher-order effective Hamiltonians in future work.