Beyond Trotterization: Variational Product Formulas for Quantum Simulation

TL;DR

This work tackles the high-depth barrier of digital quantum simulation by replacing conventional Trotter–Suzuki time steps with a variational, unitary product-form for time evolution derived from an action principle. The method optimizes a parameterized unitary $U_a(t)$ to satisfy the operator Schrödinger equation globally, yielding equations of motion for the parameters and enabling efficient, gate-friendly circuit implementations. Across a simple two-level system, the quantum Ising model, and the XXZ chain, the variational approach consistently outperforms TS (and even Ruth’s higher-order schemes) in accuracy while often reducing gate counts; analytic cubic corrections and Krylov-based partial traces further enhance practicality. The framework offers a flexible, hardware-conscious pathway to accurate long-time quantum dynamics, with broad applicability to general time-evolution operators and potential as a tool for error mitigation and efficient gate design.

Abstract

We propose a variational alternative to the Trotter-Suzuki decomposition that provides greater control over errors while preserving the unitary structure of time evolution. The variational parameters in our ansatz are derived from a global action principle, where Euler-Lagrange equations govern their optimal dynamics. Unlike conventional wavefunction-based variational methods, our approach specifically targets the time evolution operation and this allows a single set of optimized parameters to be applied to any initial state for a fixed Hamiltonian avoiding costly optimization procedures. Our method outperforms the standard Trotter-Suzuki formulas, typically achieving higher accuracy than higher-order Suzuki schemes. This translates directly to quantum computing applications, where it enables the design of quantum circuits with fewer gates which reduces noise and improves precision. Although we focus on quantum dynamics, the method is broadly applicable to problems involving general time-evolution operators. Applied to various model Hamiltonians, our approach reduces errors by factors of 2 to 5 compared to Trotter-Suzuki decompositions, demonstrating its promise for accurate quantum simulation with improved efficiency. In certain cases, the variational ansatz achieves higher accuracy than more complex higher-order Suzuki formulas while reducing the gate count by nearly half within a single circuit layer. Furthermore, we derive approximate analytical expressions for the variational parameters up to cubic order in time, valid for generic Hamiltonians. These approximations enable long-time quantum simulations with improved accuracy over equivalent Suzuki decompositions, providing ready-to-use evolution formulas that match Suzuki's gate complexity while delivering better performance.

Figures (13)

• Figure 1: Flowchart outlining the main steps of our variational algorithm. Initially, we start with a parameterize-time evolution ansatz that reflects the required complexity of the single layer in the quantum circuit. Next, the classical stage determines the variational parameters based on the action principle. Afterwards, we perform error analysis to determine the time interval where the ansatz is valid. Lastly, we generate the optimal parameters that will be then used as an input for running quantum simulations.
• Figure 2: (a) A plot of the normalized Frobenius norm of the difference between the exact and the approximate time-evolution, $\mathcal{E}_F$ (Eq. \ref{['eqn:error_metric']}), for the two-level model Eq. \ref{['eqn:2levelsH']} using the variational ansatz $U_a(t)=e^{ic_0A}e^{ic_1B}$ in comparison with the first order Trotter-Suzuki (TS) decomposition for both operator orderings AB and BA. (b) Similar to (a), we employ the 3-exponential variational ansatz $U_a(t)=e^{ic_0A}e^{ic_1B}e^{ic_2A}$ and compute the normalized Frobenius error against the symmetric TS formula Eq. \ref{['eqn:symmTS']} for the operator orderings ABA and BAB. Parameters: $h_x=5$, $h_z=2$.
• Figure 3: (a,c) Dynamical magnetization for the QIM [Eq. \ref{['eqn:QIM']}] comparing exact calculations, the symmetric TS formula [Eq. \ref{['eqn:symmTS']}], and our variational ansatz [Eq. \ref{['eqn:ansatz_V2']}]. (b,d) Relative error in magnetization comparing TS and variational methods with operator ordering BAB for stroboscopic time steps of sizes $\tau=0.2$ and $\tau=0.4$, respectively. Parameters: $h_x = h_z = J = 1.0$, $N=5$, initial state $|\psi_0\rangle = |{\uparrow}\rangle^{\otimes N}$.
• Figure 4: (a) Normalized Frobenius norm difference between exact time evolution and approximate methods for the QIM [Eq. \ref{['eqn:QIM']}]: variational ansatz [Eq. \ref{['eqn:ansatz_V3']}] and Ruth's formula [Eq. \ref{['eqn:ruth']}]. Operator orderings: variational (ABAB, BABA) and Ruth (ABABABA, BABABAB). (b) The stroboscopic error when using the variational ansatz and the Ruth decomposition with time steps of sizes $\tau_{\rm var}=0.6$, and $\tau_{\rm suz}=0.5$, respectively. Parameters: $h_x=h_z=J=1.0$, $N=5$
• Figure 5: (a,c) The dynamical magnetization for the QIM [Eq. \ref{['eqn:QIM']}] obtained via exact time evolution, Ruth formula [Eq. \ref{['eqn:ruth']}], and four exponentials variational ansatz [Eq. \ref{['eqn:ansatz_V3']}] for different operator orderings. For the variational ansatz with ordering BABA; for Ruth with operator ordering $(BA)_3B$ (denoting BABABAB). (b,d) The corresponding relative error in the dynamical magnetization at $\tau=0.5$ and $\tau=1.0$, respectively. Parameters: $h_x=h_z=J=1.0$, and $N=5$, with the initial state $|\psi_0\rangle=|\uparrow\rangle^{\otimes N}$.