Error bounds for splitting methods in unitary problems

Abstract: Splitting methods constitute a widely used class of numerical integrators for ordinary and partial differential equations, particularly well suited to problems that can be decomposed into simpler subproblems. High-order splitting schemes are available that achieve high accuracy while preserving key qualitative properties of the underlying dynamical system, and are successfully used across a broad range of fields. In this work, we present a systematic analysis of both local and global errors arising from arbitrary splitting methods applied to unitary problems. Two complementary types of error estimates are derived. The first is expressed in terms of operator norms, while the second is formulated using norms of commutators and can, under suitable assumptions, be extended to certain classes of unbounded operators. Special attention is devoted to the case where only two operators are involved. The theoretical results are illustrated by deriving explicit error bounds for some representative schemes.

• The paper presents a rigorous analysis of local and global error bounds for splitting methods applied to unitary evolution operators.
• It establishes both operator norm-based and commutator-based error estimates, with special improvements shown for symmetric high-order and two-operator splittings.
• The results provide actionable insights for optimizing integrator designs in quantum simulation and Hamiltonian systems by quantifying error accumulation over multiple steps.

Error Bounds for Splitting Methods in Unitary Problems

Overview and Motivation

Splitting methods are an essential class of geometric numerical integrators for ordinary and partial differential equations (ODEs, PDEs) that allow for efficient, structure-preserving simulations across quantum dynamics, Hamiltonian systems, stochastic processes, and sampling algorithms. This work provides a rigorous and systematic analysis of local and global error bounds for arbitrary splitting methods applied to unitary problems, i.e., problems whose evolution operators are unitary due to the underlying skew-adjoint generators. The analysis delivers two complementary error estimates: one based on operator norms and the other utilizing norms of nested commutators, with explicit attention to cases involving both bounded and unbounded operators. Special focus is placed on symmetric high-order methods and two-operator splittings.

Formal Framework for Splitting Algorithms

Consider the linear initial value problem for $u(t) \in \mathcal{H}$,

$$\frac{du}{dt} = A u = \sum_{j=1}^N A_j u, \qquad u(0) = u_0$$

where $A_j$ are (typically) skew-adjoint operators on a complex Hilbert space. The evolution operator $U(t) = \exp(tA)$ is unitary.

In practice, $A$ often decomposes additively into subproblems whose flows $\exp(tA_j)$ are easier to compute or simulate. The first-order Lie–Trotter splitting,

$u_{k+1} = \exp(hA_1)\exp(hA_2)\cdots\exp(hA_N)u_k,$

and the second-order Strang symmetric splitting,

$$S_2(h) = \exp\left(\tfrac{h}{2}A_1\right)\exp\left(\tfrac{h}{2}A_2\right)\cdots\exp\left(hA_N\right)\cdots\exp\left(\tfrac{h}{2}A_2\right)\exp\left(\tfrac{h}{2}A_1\right),$$

are the prototypical schemes. Higher-order schemes are constructed by recursively composing basic splittings with carefully chosen coefficients.

The central technical challenge is the quantification and control of the error between the exact evolution and the splitting approximation. This is nontrivial, particularly at high order, for unbounded operators, or over many integration steps.

Operator Norm-Based Error Bounds

The first family of results bounds the local and global error in terms of the operator norms $\|A_j\|$ and the splitting coefficients.

Let $\Psi(h)$ be a product-formula splitting of order $$\frac{du}{dt} = A u = \sum_{j=1}^N A_j u, \qquad u(0) = u_0$$0, i.e., $$\frac{du}{dt} = A u = \sum_{j=1}^N A_j u, \qquad u(0) = u_0$$1. Define $$\frac{du}{dt} = A u = \sum_{j=1}^N A_j u, \qquad u(0) = u_0$$2 and set $$\frac{du}{dt} = A u = \sum_{j=1}^N A_j u, \qquad u(0) = u_0$$3. For general (not necessarily symmetric) splittings, the local error satisfies

$$\frac{du}{dt} = A u = \sum_{j=1}^N A_j u, \qquad u(0) = u_0$$4

For Lie–Trotter and Strang compositions, the error bounds admit explicit dependence on the $$\frac{du}{dt} = A u = \sum_{j=1}^N A_j u, \qquad u(0) = u_0$$5-norm of the coefficient vector.

For palindromic (time-symmetric) compositions, sharper bounds are proven, improving the leading constant by a factor $$\frac{du}{dt} = A u = \sum_{j=1}^N A_j u, \qquad u(0) = u_0$$6, and reflecting the cancellation of error terms due to symmetry.

Global error after $$\frac{du}{dt} = A u = \sum_{j=1}^N A_j u, \qquad u(0) = u_0$$7 steps is bounded as $$\frac{du}{dt} = A u = \sum_{j=1}^N A_j u, \qquad u(0) = u_0$$8 times the local bound, exploiting the unitarity of the numerical and exact propagators.

Refinements: Expansion coefficients are precisely computed using multi-index notation, yielding sharper bounds, particularly for high-order methods, when truncating at higher orders or for special symmetric compositions.

Explicit Numerical Comparison of Bounds

A key highlight is the substantial improvement attained by symmetric and refined error formulas over generic bounds. For instance, in the specific case of a 6th-order, 13-stage symmetric composition (method of Sofroniou and Spaletta [sofroniou05dos]), the local error bound constants achieve improvements by more than three orders of magnitude compared to the generic estimate for a wide range of normalized time-steps.

Figure 1: Error bounds for the time-symmetric 6th-order scheme, displaying the marked improvement of time-symmetric and refined estimates over nonsymmetric bounds for varying normalized step $$\frac{du}{dt} = A u = \sum_{j=1}^N A_j u, \qquad u(0) = u_0$$9.

Commutator-Norm-Based Bounds and Lie Algebraic Structure

The second family of results provides error bounds in terms of commutators, enabling sharp analysis especially when the operators nearly commute. The central tool is the Baker–Campbell–Hausdorff (BCH) expansion and constructions based on the free Lie algebra generated by $A_j$0.

For instance, for the Strang splitting with two operators,

$A_j$1

the optimal local error bound is

$A_j$2

The coefficient constants are shown to be optimal, inherited from the BCH expansion.

For higher-order and multi-stage symmetric compositions with two operators (Runge–Kutta–Nyström-type), the error bound is expanded in a basis of the free Lie algebra, in which each term can be explicitly tracked. This yields particularly tight bounds in applications where structural commutators vanish: for example, when $A_j$3, as in classical Hamiltonian systems with quadratic kinetic energy or the time-dependent Schrödinger equation.

The analysis also covers the case where $A_j$4 or $A_j$5 are unbounded (under precise invariance and regularity assumptions), extending the rigorous applicability of the error bounds.

Numerical Implications and Optimization of Schemes

The explicit analytic form of the error as a function of operator norms and structural commutators facilitates the a priori comparison and selection of integrator schemes, as well as their robust application in applications demanding strict error control, including:

• Quantum simulation via product formulas, where minimizing the error bound directly translates to resource savings in quantum circuit complexity.
• Hamiltonian systems and molecular dynamics, where structure preservation and accuracy are critical.
• The design of new integrators, where the error constant minimization can be used as an objective in the optimizer over admissible composition coefficients.

In practice, the smallest effective error (for fixed order) is obtained by optimizing the 1-norm of the splitting coefficients, a principle that is now rigorously justified in the theoretical error bounds.

Directions for Further Development

The methodology in this work suggests multiple avenues for further research:

• Extending the sharp commutator norm bounds to a wider class of unbounded operator problems, with quantitative tracking of domain invariance and regularity.
• Analysis of splitting methods employing complex coefficients, especially in parabolic and unitary settings.
• Extension to non-linear evolution equations, leveraging Lie formalism for rigorous error analysis.
• Application to concrete quantum simulation models for physically relevant Hamiltonians, with problem-tailored error and resource estimates.

Conclusion

This work establishes a comprehensive mathematical foundation for the quantitative error analysis of splitting and composition schemes in finite- and infinite-dimensional unitary settings. By providing both operator-norm and commutator-based error bounds, it offers precise tools for method design, selection, and rigorous error certification—especially in time-symmetric high-order methods or when nested commutators vanish structurally. The explicit numerical illustration demonstrates the strong practical impact and the importance of algebraic structure in error behavior.

Figure 1: Error bounds for the time-symmetric 6th-order scheme, demonstrating the superiority of refined symmetric and commutator-based estimates compared to generic bounds.