Multi-Product Formula (MPF) in Hamiltonian Simulation

• Multi-product formula (MPF) is a construction that approximates quantum time evolution by combining weighted, low-order product formulas to cancel error terms.
• MPF employs Vandermonde-type linear cancellation constraints to remove successive error orders, bridging Trotterization and linear-combination-of-unitaries methods.
• Variants like dynamic and dual-channel MPF reduce circuit depth and improve noise robustness, offering efficient scaling with evolution time and precision.

Multi-product formula (MPF) usually denotes a construction in which a target object is approximated or expanded as a linear combination of lower-order product expressions. In quantum Hamiltonian simulation, the term refers to approximations of the time-evolution operator $U(t)=e^{-iHt}$ obtained by combining product-formula circuits evaluated at different effective step sizes so that successive terms in the Trotter error expansion cancel. In this setting, MPFs interpolate between product-formula methods and linear-combination-of-unitaries techniques: they retain the low-order circuit structure of Trotterization, but aim for near-linear dependence on evolution time and poly-logarithmic dependence on precision (Aftab et al., 2024).

1. Formal construction

For a time-independent Hamiltonian decomposed as $H=\sum_{\gamma=1}^{\Gamma}H_\gamma$, an MPF starts from a low-order product formula $U_p(s)$ or $S_\chi(s)$ and forms a weighted sum of repeated short-time evolutions. Two equivalent notational forms appear in the literature: $$U_{\mathrm{MPF}}(t)=\sum_{j=1}^M c_j\,U_p(t/k_j)^{k_j}, \qquad M_{l,\chi}(t)=\sum_{j=1}^l a_j\,S_\chi^{\,k_j}(t/k_j).$$ Here the $k_j$ are positive integers and the coefficients are chosen so that low-order terms in the expansion of the product formula cancel (Aftab et al., 2024).

The cancellation conditions are Vandermonde-type linear constraints. In the hybrid hardware-oriented formulation, one imposes

$$\sum_{j=1}^l a_j=1, \qquad \sum_{j=1}^l \frac{a_j}{k_j^\eta}=0$$

for a prescribed range of exponents $\eta$, which removes the first $l-1$ error terms. For a second-order base formula, solving the linear system with rows $1,k_j^{-2},\ldots,k_j^{-2m+2}$ yields a short-interval global error of order $H=\sum_{\gamma=1}^{\Gamma}H_\gamma$0 (Vazquez et al., 2022).

This algebraic structure is closely related to Richardson extrapolation. In the notation of one later analysis, an MPF of “order” $H=\sum_{\gamma=1}^{\Gamma}H_\gamma$1 is written as

$H=\sum_{\gamma=1}^{\Gamma}H_\gamma$2

with integer dilation factors $H=\sum_{\gamma=1}^{\Gamma}H_\gamma$3 and real weights $H=\sum_{\gamma=1}^{\Gamma}H_\gamma$4 fixed so as to cancel all terms through $H=\sum_{\gamma=1}^{\Gamma}H_\gamma$5 (Mizuta, 9 Jul 2025).

2. Error cancellation and commutator structure

A central question is whether MPF error bounds can reflect the same locality- and commutator-sensitive structure that makes product formulas attractive. In the rigorous analysis of the well-conditioned MPF, the relevant quantities are the nested-commutator sums

$H=\sum_{\gamma=1}^{\Gamma}H_\gamma$6

For one short step $H=\sum_{\gamma=1}^{\Gamma}H_\gamma$7, the MPF error is bounded by a series whose terms are weighted by products of these nested-commutator sums. In the commuting case, all $H=\sum_{\gamma=1}^{\Gamma}H_\gamma$8 vanish for $H=\sum_{\gamma=1}^{\Gamma}H_\gamma$9, and the MPF is exact (Aftab et al., 2024).

That result established explicit commutator scaling, but it also prompted a technical dispute about how system-size dependence should be interpreted. One later paper argues that the nested-commutator bound expressed through arbitrarily large $U_p(s)$0 does not by itself resolve size-efficient complexity, because locality benefits are absent when arbitrarily high nested commutators are required. The same work therefore introduces an alternative analysis based on a truncated BCH or Floquet-Magnus expansion, with truncation order

$U_p(s)$1

and derives an MPF error bound involving only commutators up to order $U_p(s)$2 (Mizuta, 9 Jul 2025).

A common simplification is to describe MPF as merely “higher-order Trotter.” That characterization is incomplete. The defining feature is not a single higher-order product formula, but a linear combination of low-order product formulas whose error terms are arranged to cancel. This distinction matters because the resulting bounds and implementation strategies differ from those of ordinary Suzuki recursion.

3. Complexity and asymptotic scaling

The main asymptotic promise of MPF is a combination of near-optimal dependence on evolution time $U_p(s)$3 and target precision $U_p(s)$4 with commutator-sensitive dependence on the Hamiltonian decomposition. In the explicit-commutator analysis, one defines an MPF commutator parameter

$U_p(s)$5

If these stabilize to $U_p(s)$6, then the MPF query cost is

$U_p(s)$7

which is almost linear in $U_p(s)$8 and poly-logarithmic in $U_p(s)$9 (Aftab et al., 2024).

This contrasts with conventional Trotter–Suzuki simulation. A $S_\chi(s)$0th-order product formula has cost

$S_\chi(s)$1

which is super-linear in $S_\chi(s)$2 and polynomial in $S_\chi(s)$3. In the comparison drawn in the same analysis, post-Trotter methods such as qubitization and quantum signal processing have cost

$S_\chi(s)$4

which is optimal in $S_\chi(s)$5 but scales with the spectral norm $S_\chi(s)$6 and cannot exploit commutator cancellations (Aftab et al., 2024).

The later truncated-commutator analysis seeks to recover size-efficient locality dependence rigorously. For a $S_\chi(s)$7-local, $S_\chi(s)$8-extensive Hamiltonian on $S_\chi(s)$9 qubits, it derives a total query cost

$$U_{\mathrm{MPF}}(t)=\sum_{j=1}^M c_j\,U_p(t/k_j)^{k_j}, \qquad M_{l,\chi}(t)=\sum_{j=1}^l a_j\,S_\chi^{\,k_j}(t/k_j).$$0

and interprets this as polynomial in $$U_{\mathrm{MPF}}(t)=\sum_{j=1}^M c_j\,U_p(t/k_j)^{k_j}, \qquad M_{l,\chi}(t)=\sum_{j=1}^l a_j\,S_\chi^{\,k_j}(t/k_j).$$1 with exponent $$U_{\mathrm{MPF}}(t)=\sum_{j=1}^M c_j\,U_p(t/k_j)^{k_j}, \qquad M_{l,\chi}(t)=\sum_{j=1}^l a_j\,S_\chi^{\,k_j}(t/k_j).$$2 together with polylogarithmic dependence on $$U_{\mathrm{MPF}}(t)=\sum_{j=1}^M c_j\,U_p(t/k_j)^{k_j}, \qquad M_{l,\chi}(t)=\sum_{j=1}^l a_j\,S_\chi^{\,k_j}(t/k_j).$$3 (Mizuta, 9 Jul 2025). This suggests that the precise form of “commutator scaling” for MPF is now an active subject of refinement rather than a closed matter.

4. Implementation models and conditioning

The original MPF proposal is naturally implemented through a linear combination of unitaries. One prepares an ancilla superposition encoding the weights, applies controlled product-formula unitaries, and uncomputes the ancilla so that the desired linear combination is obtained on success. In this formulation, the success probability scales as $$U_{\mathrm{MPF}}(t)=\sum_{j=1}^M c_j\,U_p(t/k_j)^{k_j}, \qquad M_{l,\chi}(t)=\sum_{j=1}^l a_j\,S_\chi^{\,k_j}(t/k_j).$$4, which makes conditioning of the coefficient vector central (Vazquez et al., 2022).

A hybrid implementation avoids coherent LCU altogether. Instead of constructing the linear combination unitarily, one runs the distinct product-formula circuits independently, measures the same observable on each circuit, and classically combines the expectation values. The paper introducing this approach states that it has the same approximation bounds as the fully quantum MPFs, but requires no additional qubits, no controlled operations, and is not probabilistic (Vazquez et al., 2022).

Conditioning is quantified by the coefficient 1-norm, $$U_{\mathrm{MPF}}(t)=\sum_{j=1}^M c_j\,U_p(t/k_j)^{k_j}, \qquad M_{l,\chi}(t)=\sum_{j=1}^l a_j\,S_\chi^{\,k_j}(t/k_j).$$5. Large $$U_{\mathrm{MPF}}(t)=\sum_{j=1}^M c_j\,U_p(t/k_j)^{k_j}, \qquad M_{l,\chi}(t)=\sum_{j=1}^l a_j\,S_\chi^{\,k_j}(t/k_j).$$6 amplifies both sampling error and hardware noise in the final post-processed estimate, so well-conditioned MPFs are constructed by searching over integer tuples $$U_{\mathrm{MPF}}(t)=\sum_{j=1}^M c_j\,U_p(t/k_j)^{k_j}, \qquad M_{l,\chi}(t)=\sum_{j=1}^l a_j\,S_\chi^{\,k_j}(t/k_j).$$7 and discarding those whose $$U_{\mathrm{MPF}}(t)=\sum_{j=1}^M c_j\,U_p(t/k_j)^{k_j}, \qquad M_{l,\chi}(t)=\sum_{j=1}^l a_j\,S_\chi^{\,k_j}(t/k_j).$$8 exceeds a chosen threshold. In the hardware-friendly study, the authors report that for $$U_{\mathrm{MPF}}(t)=\sum_{j=1}^M c_j\,U_p(t/k_j)^{k_j}, \qquad M_{l,\chi}(t)=\sum_{j=1}^l a_j\,S_\chi^{\,k_j}(t/k_j).$$9 one can in all practical cases achieve $k_j$0 for $k_j$1-based MPFs or $k_j$2 for $k_j$3-based MPFs (Vazquez et al., 2022).

That same work integrates MPF with three mitigation layers intended for noisy devices: Pauli Twirling, pulse-efficient transpilation, and zero-noise extrapolation based on scaled cross-resonance pulses. The intended role of MPF in this setting is not only asymptotic improvement, but also circuit-depth reduction under restricted hardware resources.

5. Dynamic and dual-channel variants

One line of development treats the MPF coefficients as time-dependent variables rather than fixed extrapolation weights. Dynamic MPF introduces coefficients $k_j$4 chosen to minimize a computable proxy for the Frobenius-norm projection error. Because the exact coupling term depends on the unknown target state, it is replaced by an approximate update derived from previous time steps, yielding a discrete-time linear system for the coefficient vector (Zhuk et al., 2023).

To stabilize this procedure in the presence of uncertainty, the same work proposes Minimax MPF. At each time step, the coefficients are obtained from a convex program with an $k_j$5 penalty,

$k_j$6

where $k_j$7 models bounded hardware and sampling noise. The resulting estimator is accompanied by a rigorous error bound and is described as robust to both algorithmic Trotter errors and bounded sampling and hardware noise (Zhuk et al., 2023).

A different recent variant is the dual-channel multi-product formula (DCMPF). It combines a regular product formula $k_j$8 with the reversed-sequence formula

$k_j$9

and defines

$$\sum_{j=1}^l a_j=1, \qquad \sum_{j=1}^l \frac{a_j}{k_j^\eta}=0$$0

Because the forward-plus-reversed combination is symmetric under $$\sum_{j=1}^l a_j=1, \qquad \sum_{j=1}^l \frac{a_j}{k_j^\eta}=0$$1, the error expansion contains only even powers of $$\sum_{j=1}^l a_j=1, \qquad \sum_{j=1}^l \frac{a_j}{k_j^\eta}=0$$2, and the method cancels twice as many orders with the same number of channels $$\sum_{j=1}^l a_j=1, \qquad \sum_{j=1}^l \frac{a_j}{k_j^\eta}=0$$3 (Park et al., 2 Feb 2026).

For a conventional MPF built from a regular $$\sum_{j=1}^l a_j=1, \qquad \sum_{j=1}^l \frac{a_j}{k_j^\eta}=0$$4th-order product formula, the leading error scales as $$\sum_{j=1}^l a_j=1, \qquad \sum_{j=1}^l \frac{a_j}{k_j^\eta}=0$$5; for DCMPF it scales as $$\sum_{j=1}^l a_j=1, \qquad \sum_{j=1}^l \frac{a_j}{k_j^\eta}=0$$6. The same analysis states that, for the same algorithmic error, the dual-channel construction uses approximately half the circuit depth of a symmetric-PF MPF, while the sampling error remains essentially unchanged (Park et al., 2 Feb 2026).

6. Representative applications and benchmarks

Several studies use MPF to quantify trade-offs among system size, evolution time, precision, circuit depth, and hardware noise. The examples span asymptotic analyses, numerical benchmarks, and small-scale experiments.

Setting	Reported MPF result	Source
Plane-wave electronic structure on $$\sum_{j=1}^l a_j=1, \qquad \sum_{j=1}^l \frac{a_j}{k_j^\eta}=0$$7 spin-orbitals	$$\sum_{j=1}^l a_j=1, \qquad \sum_{j=1}^l \frac{a_j}{k_j^\eta}=0$$8, $$\sum_{j=1}^l a_j=1, \qquad \sum_{j=1}^l \frac{a_j}{k_j^\eta}=0$$9, and MPF cost $\eta$0	(Aftab et al., 2024)
1D Heisenberg chain benchmark	Gate counts scaling roughly $\eta$1–$\eta$2, improving over the $\eta$3-scaling of Trotter	(Aftab et al., 2024)
Spin-boson resource study	Deepest circuit has $\eta$4; reductions up to a factor of 8 for modest systems	(Vazquez et al., 2022)
Five-qubit transverse-field Ising model	Final errors $\eta$5 with only 2–4 Trotter steps; a single PF with $\eta$6 needed for similar accuracy; up to $\eta$7 depth reduction	(Vazquez et al., 2022)
Heisenberg spin chain with dynamic MPF analysis	MPF can reduce required circuit depth by orders of magnitude in the near-term regime $\eta$8	(Zhuk et al., 2023)

In the plane-wave example, the explicit-commutator analysis compares second-order Trotter, post-Trotter methods, and MPF. The reported scalings are $\eta$9 for second-order Trotter, $l-1$0 for post-Trotter methods such as qubitization, and $l-1$1 for MPF. The stated conclusion is a polynomial speedup in $l-1$2 and an exponential speedup in $l-1$3 over second-order Trotter, with competitiveness against qubitization when commutator structure is favorable (Aftab et al., 2024).

For near-term hardware, the five-qubit transverse-field Ising demonstration is significant because it pairs well-conditioned MPF coefficients with explicit error-mitigation procedures. The reported well-conditioned pairs are $l-1$4, $l-1$5, and $l-1$6, with weights $l-1$7, $l-1$8, and $l-1$9. These combinations produced errors of order $1,k_j^{-2},\ldots,k_j^{-2m+2}$0 using only shallow circuits, whereas a single product formula required substantially more Trotter steps to reach similar accuracy (Vazquez et al., 2022).

7. Other mathematical uses of the term

Although current usage on arXiv is dominated by Hamiltonian simulation, “multi-product formula” is not exclusive to quantum algorithms. In stochastic analysis, the term denotes an explicit product expansion for multiple stochastic integrals with respect to the compensated random measure of a Lévy process. Under an integrability condition on contracted kernels, the product

$1,k_j^{-2},\ldots,k_j^{-2m+2}$1

is expanded as a finite sum over contraction data $1,k_j^{-2},\ldots,k_j^{-2m+2}$2: $1,k_j^{-2},\ldots,k_j^{-2m+2}$3 That formula is used to derive moments, cumulants, mixed moments, and central limit theorems (Tella et al., 2023).

In algebraic combinatorics, the term also appears in the theory of multivariate Rogers–Szegö polynomials. There the product

$1,k_j^{-2},\ldots,k_j^{-2m+2}$4

is expanded as

$1,k_j^{-2},\ldots,k_j^{-2m+2}$5

with recursively defined coefficients $1,k_j^{-2},\ldots,k_j^{-2m+2}$6. In that setting, the formula is reinterpreted in symmetric-function terms as a statement about structure constants (Cameron et al., 2013).

These non-quantum usages are terminologically related but mathematically independent. In contemporary quantum-computing literature, MPF almost always refers to the linear-combination construction for Hamiltonian simulation; in other fields, it names a product-expansion identity specific to the objects under study.