High-order splitting of non-unitary operators on quantum computers

1. Problem Statement

Quantum computers are inherently designed to simulate unitary dynamics (reversible time evolution of Hamiltonians). However, most practical problems in computational science involve dissipative and irreversible processes (e.g., friction, diffusion, viscosity), which are described by non-unitary dynamics.

• The Challenge: Standard high-order operator splitting methods (used to improve accuracy beyond simple first or second-order approximations) rely on negative time steps (negative coefficients) to achieve higher orders.
• The Limitation: While negative time steps are manageable for unitary dynamics (by reversing the evolution), they are numerically unstable for dissipative dynamics. Reversing a dissipative process (e.g., un-diffusing heat) is physically impossible and cannot be block-encoded into unitary operators without rescaling, leading to instability.
• Current State: Existing quantum simulations of dissipative systems are largely limited to second-order accuracy (Strang splitting) because higher-order methods have not been successfully adapted for non-unitary operators.

2. Methodology

The authors propose a stable quantum algorithm that utilizes complex-coefficient product formulas to decompose dissipative dynamics into a sequence of simple Hamiltonian evolutions.

Core Concept: Complex Coefficients

Instead of using real coefficients (which force negative time steps for high orders), the authors employ product formulas where:

• Unitary Substages: Use real, positive coefficients ($b_j > 0$). These correspond to standard unitary evolutions ($e^{iH_2 b_j \Delta t}$).
• Dissipative Substages: Use complex coefficients ($a_j \in \mathbb{C}$) where the real part is positive ($\Re(a_j) > 0$).
  • The real part ($\Re(a_j)$) drives the dissipative evolution in imaginary time ($e^{H_1 \Re(a_j) \Delta t}$), which preserves the contractive nature of the dynamics.
  • The imaginary part ($\Im(a_j)$) drives an additional unitary evolution ($e^{iH_1 \Im(a_j) \Delta t}$).

This decomposition allows the total evolution $e^{(H_1 + iH_2)t}$ to be approximated by a sequence of simple, stable steps:
$$S_j(\Delta t) = e^{iH_2 b_j \Delta t} e^{iH_1 \Im(a_j) \Delta t} e^{H_1 \Re(a_j) \Delta t}$$

Specific Formulas Used

The paper utilizes two specific high-order product formulas that satisfy the constraints ($\Re(a_j) > 0, b_j > 0$):

1. Order 4: Based on Castella et al. (5 stages).
2. Order 6: Based on Bernier et al. (16 stages, symmetric-conjugate).

3. Application and Circuit Derivation

The method was demonstrated using the classical damped wave equation:
$$\frac{\partial^2 \psi}{\partial t^2} + \gamma \frac{\partial \psi}{\partial t} = c^2 \nabla^2 \psi$$

• Discretization: The equation was transformed into a first-order system of ODEs using a pseudo-spectral method (Fourier modes).
• Hamiltonian Decomposition: The generator $M$ was split into a Hermitian part (damping, $H_1$) and an anti-Hermitian part (wave propagation, $H_2$).
• Circuit Implementation:
  • Wave Propagation: Implemented using controlled $R_Y$ rotations based on the binary expansion of the wavenumber.
  • Damping: Implemented via post-selection. An ancilla qubit is rotated based on the damping coefficient; measuring the ancilla in the $|0\rangle$ state projects the system into the damped state.
  • Imaginary Time Damping: The complex coefficients introduce phase gates ($P(\theta)$) to handle the imaginary components of the damping coefficients.

4. Key Results

Statevector Simulations (Noise-Free)

• Accuracy Scaling: Simulations confirmed that the error scales as $O(\Delta t^p)$ for orders $p=1, 2, 4, 6$.
• Efficiency: While low-order methods require fewer gates for coarse accuracy, high-order methods (Order 4 and 6) become significantly more efficient for high-precision requirements.
• Gate Count: Order 6 achieved classical machine precision with approximately $10^4$ CNOT gates, whereas lower orders would require exponentially more steps to reach the same precision.

Hardware Execution (IonQ Forte 1)

The algorithm was executed on a 36-qubit trapped-ion quantum processor.

• Setup: Simulated a damped wave with $N=16$ grid points.
• Performance Comparison:
  • Order 1 & 2: Order 1 was visibly inaccurate. Order 2 showed good accuracy but was less precise than Order 4.
  • Order 4: Achieved the best balance between accuracy and noise susceptibility. It outperformed Order 2 in terms of $\ell_2$ error norm (0.1138 vs 0.1166) despite having a deeper circuit.
  • Order 6: Suffered from noise dominance. The increased circuit depth (232 CNOT gates per step) introduced too much hardware noise, resulting in lower accuracy (0.2207 error) than even Order 1.
• Post-selection Success: The success probability of the post-selection (measuring ancilla in $|0\rangle$) remained high (approx. 70–87%), validating the stability of the approach.

5. Significance and Contributions

1. Breaking the Order-2 Barrier: This work provides the first demonstration of stable high-order (Order 4 and 6) operator splitting for non-unitary, dissipative dynamics on quantum hardware.
2. Complex Coefficient Framework: It establishes a rigorous framework for using complex coefficients to convert unstable negative-time dissipative steps into stable positive-real-time and unitary imaginary-time steps.
3. Practical Viability on NISQ: The results demonstrate that Order 4 is the "sweet spot" for near-term quantum devices. It offers a significant accuracy advantage over Order 2 without the prohibitive noise sensitivity of Order 6.
4. Resource Scaling Analysis: The paper derives the gate complexity scaling as $O(2^{p/2} t^{1+1/p} \epsilon^{-1/p})$. This shows that for high-precision simulations, high-order splitting is asymptotically superior to low-order methods, reducing the dependence on error tolerance ($\epsilon$) from $\epsilon^{-1/2}$ (Order 2) to $\epsilon^{-1/6}$ (Order 6).
5. Industrial Applicability: The method is applicable to a wide range of problems involving diffusion, advection, and fluid dynamics, extending the utility of quantum computing beyond pure quantum chemistry and physics to classical computational science.

Conclusion

The paper successfully demonstrates that high-order operator splitting, previously thought impossible for dissipative systems due to instability, can be realized on quantum computers using complex-coefficient product formulas. While hardware noise currently limits the practical order to 4, the approach offers a robust path toward accurate simulation of dissipative dynamics on near-term quantum processors.