High-Order Product Formulas Overview

• High-order product formulas are structured identities that express composite operators as integrals or ordered products, supporting both harmonic analysis and quantum simulation applications.
• They combine analytical methods—using special functions, Gegenbauer polynomials, and Bessel functions—with recursive Lie–Trotter–Suzuki techniques to achieve precise error cancellation.
• These formulas are applied in diverse contexts such as digital quantum simulation, stochastic integral computations, and quantum chemistry, highlighting trade-offs between order, cost, and simulation accuracy.

High-order product formulas are structured identities or approximation schemes that rewrite a composite object in terms of simpler building blocks. In harmonic analysis, a product formula is an identity of the form $$K(\lambda,x)K(\lambda,y)=\int K(\lambda,z)\,d\mu_{x,y}(z)$$, with the spectral parameter $\lambda$ separated from the measure $d\mu_{x,y}$. In numerical Hamiltonian simulation, a product formula approximates $e^{-iHt}$ by ordered products of exponentials of simpler Hamiltonian fragments. Related constructions also approximate exponentials of commutators and nested commutators. The phrase “high-order” is therefore used in more than one technical sense: in generalized Hankel and Dunkl analysis it refers to the integer deformation parameter $n$, whereas in Lie–Trotter–Suzuki theory it refers to cancellation of BCH terms so that the local error is $\mathcal O(t^{k+1})$ for an order-$k$ formula (Boubatra et al., 2020, Lopez-Cerezo, 14 Jul 2025).

1. Terminology and structural principles

In the harmonic-analytic sense, a product formula expresses the pointwise product of two kernel evaluations as an integral superposition of the same kernel. This structure is central because it allows one to define generalized translation and generalized convolution, and to prove Plancherel and inversion formulas and analyze positivity and hypergroup structures (Boubatra et al., 2020).

In the Hamiltonian-simulation sense, a product formula has the generic form

$$S(t)=\prod_{i=1}^{q} e^{p_{i,1}H_1 t}\cdots e^{p_{i,L}H_L t},$$

for $H=\sum_{j=1}^{L}H_j$. It is called order $k$ when

$\lambda$0

in operator norm. Repetition with step size $\lambda$1 converts local error $\lambda$2 into global error $\lambda$3 (Lopez-Cerezo, 14 Jul 2025).

A recurring structural distinction is between generic and symmetric formulas. In the Suzuki hierarchy, symmetry means $\lambda$4, which forces only odd powers of $\lambda$5 to appear in the kernel or error generator. In corrected formulas, similarity transforms $\lambda$6 preserve the time-reversal symmetry needed for recursive order-raising. In the Jordan–Banach setting, the same symmetry principle survives and implies that an odd-order symmetric approximant is automatically one order higher (Bagherimehrab et al., 2024, Chehade et al., 2024).

A common ambiguity concerns the phrase “high-order.” In the generalized Hankel setting, “high-order product formula” refers to all integers $\lambda$7 and the corresponding higher-degree Gegenbauer/Bessel structure. In quantum simulation, it refers instead to formal approximation order, typically obtained through recursive BCH cancellation. The literature treats both usages as standard, but they are mathematically distinct (Boubatra et al., 2020, Bagherimehrab et al., 2024).

2. Special-function and harmonic-analysis product formulas

For kernel transforms associated with Bessel and Dunkl analysis, the prototype identity is

$\lambda$8

A recent one-dimensional example is the generalized Hankel function $\lambda$9, the kernel of the transform $d\mu_{x,y}$0, with

$d\mu_{x,y}$1

and

$d\mu_{x,y}$2

for $d\mu_{x,y}$3 and $d\mu_{x,y}$4 (Boubatra et al., 2020).

The core analytic input is a new integral representation for

$d\mu_{x,y}$5

derived from Gegenbauer’s addition theorem, orthogonality and recurrence relations of Gegenbauer polynomials, Sonine’s integral and differentiation formulas for Bessel functions, and the three-term recurrence relation for $d\mu_{x,y}$6. This yields an explicit product formula for $d\mu_{x,y}$7: $d\mu_{x,y}$8 where $d\mu_{x,y}$9 for $e^{-iHt}$0, and the kernel $e^{-iHt}$1 is explicit in terms of the classical Bessel kernel $e^{-iHt}$2 and a Gegenbauer polynomial $e^{-iHt}$3 evaluated at a geometric parameter $e^{-iHt}$4 (Boubatra et al., 2020).

The measure $e^{-iHt}$5 is real-valued, compactly supported in a segment, has total mass $e^{-iHt}$6, and satisfies

$e^{-iHt}$7

It is not necessarily positive; explicit negative values of the kernel occur. This non-positivity is technically important because it distinguishes the resulting structure from positive hypergroup convolution, even though many Fourier-type properties remain valid (Boubatra et al., 2020).

The “high-order” aspect here is indexed by $e^{-iHt}$8. For $e^{-iHt}$9, the formula reduces to Rösler’s product formula for the rank-one Dunkl kernel; for $n$0, it reproduces the modified Hankel case of Ben Saïd; for general $n$1, the kernel is a linear combination of $n$2 and $n$3, the product formula involves a Gegenbauer polynomial of degree $n$4, and the geometry is controlled by $n$5. The associated translation operator

$n$6

is symmetric in $n$7, bounded on $n$8 with norm at most $n$9, and induces a convolution diagonalized by $\mathcal O(t^{k+1})$0 (Boubatra et al., 2020).

3. Lie–Trotter–Suzuki hierarchies for Hamiltonian simulation

For a bipartition $\mathcal O(t^{k+1})$1, the basic formulas are

$\mathcal O(t^{k+1})$2

Using BCH,

$\mathcal O(t^{k+1})$3

so PF1 has local error $\mathcal O(t^{k+1})$4, while PF2 satisfies

$\mathcal O(t^{k+1})$5

hence has local error $\mathcal O(t^{k+1})$6. For $\mathcal O(t^{k+1})$7 steps of size $\mathcal O(t^{k+1})$8, a $\mathcal O(t^{k+1})$9-th order formula has global error $k$0 (Bagherimehrab et al., 2024).

Suzuki’s recursive construction generates arbitrarily high even orders from $k$1: $k$2 These formulas satisfy

$k$3

so they are order $k$4 with local error $k$5 (Bagherimehrab et al., 2024).

For a general decomposition $k$6, the symmetric second-order base formula is

$k$7

and the number of exponentials in the recursive order-$k$8 Suzuki formula satisfies

$k$9

Rigorous bounds show that, after optimizing over the order $$S(t)=\prod_{i=1}^{q} e^{p_{i,1}H_1 t}\cdots e^{p_{i,L}H_L t},$$0, the total number of exponentials can be made near-linear in $$S(t)=\prod_{i=1}^{q} e^{p_{i,1}H_1 t}\cdots e^{p_{i,L}H_L t},$$1 and subpolynomial in $$S(t)=\prod_{i=1}^{q} e^{p_{i,1}H_1 t}\cdots e^{p_{i,L}H_L t},$$2, although the constants in such bounds are loose in practice (Lopez-Cerezo, 14 Jul 2025).

Negative coefficients are a persistent feature of high-order Suzuki formulas. The practical recursion

$$S(t)=\prod_{i=1}^{q} e^{p_{i,1}H_1 t}\cdots e^{p_{i,L}H_L t},$$3

has $$S(t)=\prod_{i=1}^{q} e^{p_{i,1}H_1 t}\cdots e^{p_{i,L}H_L t},$$4 and a refined “fractal” structure, but it still includes negative internal time segments. This is mathematically legitimate because the factors are unitary, but it remains one of the standard caveats in high-order splitting methods (Lopez-Cerezo, 14 Jul 2025).

4. Correctors, extrapolation, randomization, and typical-case theory

Corrected product formulas (CPFs) modify a standard product formula by inserting correctors that cancel selected BCH terms. The basic templates are similarity-corrected formulas

$$S(t)=\prod_{i=1}^{q} e^{p_{i,1}H_1 t}\cdots e^{p_{i,L}H_L t},$$5

and symmetric-corrected formulas

$$S(t)=\prod_{i=1}^{q} e^{p_{i,1}H_1 t}\cdots e^{p_{i,L}H_L t},$$6

For perturbed systems $$S(t)=\prod_{i=1}^{q} e^{p_{i,1}H_1 t}\cdots e^{p_{i,L}H_L t},$$7 with $$S(t)=\prod_{i=1}^{q} e^{p_{i,1}H_1 t}\cdots e^{p_{i,L}H_L t},$$8, a CPF2 built from

$$S(t)=\prod_{i=1}^{q} e^{p_{i,1}H_1 t}\cdots e^{p_{i,L}H_L t},$$9

and the corrector

$H=\sum_{j=1}^{L}H_j$0

satisfies

$H=\sum_{j=1}^{L}H_j$1

Recursive CPFs then achieve

$H=\sum_{j=1}^{L}H_j$2

while in the non-perturbed case they improve $H=\sum_{j=1}^{L}H_j$3 to $H=\sum_{j=1}^{L}H_j$4. Symplectic correctors have additive overhead because

$H=\sum_{j=1}^{L}H_j$5

so only two extra exponentials are required per full simulation (Bagherimehrab et al., 2024).

Multi-product formulas (MPFs) replace a single higher-order step by a linear combination of lower-order steps with different Trotter exponents: $H=\sum_{j=1}^{L}H_j$6 For symmetric base formulas $H=\sum_{j=1}^{L}H_j$7, the coefficients satisfy

$H=\sum_{j=1}^{L}H_j$8

The conditioning factor

$H=\sum_{j=1}^{L}H_j$9

controls both LCU success probability and noise amplification. Ill-conditioned arithmetic sequences have $k$0, whereas well-conditioned sequences can satisfy $k$1 for practical $k$2 in $k$3-based MPFs. The hybrid realization computes expectation values from separate shallow circuits and combines them classically, so it requires no additional qubits, no controlled operations, and is not probabilistic. In hardware demonstrations on TFIM, well-conditioned MPFs achieved up to an order of magnitude algorithmic error reduction and up to $k$4 circuit depth reduction (Vazquez et al., 2022).

Randomization enters at two levels. First, randomized MPFs implement a linear combination of unitary product-formula circuits by sampling from a distribution over branches rather than using oblivious amplitude amplification; this reduces the circuit depth and gives simulation error that shrinks exponentially with the circuit depth, with rigorous concentration bounds for the sampling estimator (Faehrmann et al., 2021). Second, average-case analyses of higher-order Suzuki formulas show that product-formula error can scale much better for the vast majority of input states than in worst-case operator norm. For general $k$5-local Hamiltonians, the typical-state gate complexity is governed by local and global $k$6-norm quantities rather than the worst-case local $k$7-norm, and analogous improvements extend to fermionic Hamiltonians and Gaussian-coefficient models such as SYK (Chi-Fang et al., 2021).

5. Commutators, non-associative algebras, and Lévy chaoses

Exponentials of commutators require a distinct high-order theory because the basic target is

$k$8

A foundational recursive framework begins from the group commutator

$k$9

and constructs arbitrarily high-order approximations to exponentials of commutators and nested commutators, with nearly linear scaling in the total evolution time and subpolynomial scaling in $\lambda$00 (Childs et al., 2012). A later refinement gives a direct six-exponential third-order commutator formula

$\lambda$01

higher-order recursive families labeled $\lambda$02-, $\lambda$03-, $\lambda$04-, and $\lambda$05-copy constructions, and formulas for

$\lambda$06

that include linear terms in the commutator simulation at no extra asymptotic cost. Applications include digital counterdiabatic driving, one-dimensional fermion chains with nearest- and next-nearest-neighbor hopping terms, and truncated Kapit–Mueller Hamiltonians (Chen et al., 2021).

High-order product formulas also extend beyond associative operator algebras. In Jordan–Banach algebras, the paper “Error Estimates and Higher Order Trotter Product Formulas in Jordan-Banach Algebras” constructs first-, second-, and third-order formulas and then recursive higher-order formulas with explicit norm bounds. For two elements $\lambda$07, the third-order Jordan formula

$\lambda$08

satisfies

$\lambda$09

with an explicit bound

$\lambda$10

The recursive order-raising conditions

$\lambda$11

and their symmetric counterparts reproduce the Suzuki mechanism in a non-associative setting (Chehade et al., 2024).

A different generalization appears in stochastic analysis. For multiple integrals $\lambda$12 with respect to the Brownian–Poisson random measure associated with a Lévy process, products admit an explicit finite expansion

$\lambda$13

provided suitable $\lambda$14-integrability conditions hold for all contracted kernels. In the Brownian case this reduces to the classical Wiener–Itô contraction formula; in the presence of jumps, higher-order diagonal interactions appear, and additional kernel integrability becomes necessary. The same framework yields explicit expectations, moments, cumulants, and a central limit theorem for normalized first-order integrals (Tella et al., 2023).

6. Quantum chemistry, architecture, and formula selection

In molecular ground-state energy estimation via phase estimation, deterministic higher-order product formulas are evaluated by balancing their per-step cost against the eigenvalue error model

$\lambda$15

For one-dimensional hydrogen chains $\lambda$16 through $\lambda$17, the benchmark paper “Evaluating higher-order product formulae for molecular ground-state energy estimation” compares total gate count $\lambda$18 and $\lambda$19-layer depth. Among previously considered deterministic formulas, the eighth-order construction introduced by Morales et al. minimizes both cost metrics at a chemically relevant target error. However, increasing the formal order does not automatically reduce the total cost: near chemical accuracy, the tenth-order formula introduced in the same work can be less efficient than the eighth-order one. Motivated by this, the paper constructs a new fourth-order Yoshida-type formula with

$\lambda$20

and finds that it achieves the lowest total gate count among the formulae considered for all H-chain instances near chemical accuracy and over much of the $\lambda$21-$\lambda$22 mHa target-error window for most instances, while also reducing the $\lambda$23-layer depth (Abe et al., 29 May 2026).

For realistic electronic-structure Hamiltonians, recent work re-examines Trotter methods under fault-tolerant architectures in which non-Clifford operations may be generated more locally and cheaply. The SPRINT framework—Symmetry-Protected Randomized near-Integrable Trotter—combines near-integrability, randomization, symmetry protection, use of QROM, and a Generalized Rank Decomposition (GRADE) of electronic Hamiltonians. Applied to the X-ray absorption spectrum of $\lambda$24, it reduces the Toffoli gate cost by a factor of $\lambda$25 relative to the previous state of the art, with a gate cost only $\lambda$26 higher than qubitization while requiring a dramatic $\lambda$27 fewer logical qubits (Casares et al., 29 Jun 2026).

Across these quantum-chemistry studies, formula choice is explicitly regime-dependent. The benchmark evidence shows that order alone is not an adequate guide: the decisive quantities are the number of exponentials per step, the fitted error prefactor $\lambda$28, the target error, and the compilation model. This suggests that “high-order product formula” is best understood not as a single asymptotic prescription, but as a design space in which recursive order-raising, correctors, extrapolation, near-integrability, and hardware-aware compilation are combined to match the structure of the Hamiltonian and the architecture (Abe et al., 29 May 2026, Casares et al., 29 Jun 2026).