Operator()/to_matrix() on PauliEvolutionGate silently ignores an explicitly requested approximate synthesis (e.g. LieTrotter)

[Mostafa-Atallah2020]

Environment

• Qiskit version: 2.2.3
• Python version: 3.10
• Operating system: Windows 11 (also reproduced across the range qiskit >= 2.2)

What is happening?

When a PauliEvolutionGate is constructed with an explicitly approximate product-formula synthesis, for example synthesis=LieTrotter(reps=1), taking Operator() (equivalently to_matrix()) on the un-decomposed circuit silently ignores the requested synthesis and returns the exact evolution exp(-iHt). The user explicitly asked for a first-order Trotter approximation, but the matrix that comes back has zero Trotter error, with no warning or error of any kind.

The surprising part is that an argument the user explicitly passed, synthesis=LieTrotter(...), is silently discarded by Operator()/to_matrix(). Silently ignoring a user-supplied argument and returning a result that is "more accurate" than the requested method can ever be is a hard-to-detect failure mode. It is exactly the case a user is least likely to question, and it is easy to hit when benchmarking Trotter error or comparing product formulas against other simulation methods.

How can we reproduce the issue?

import numpy as np
from scipy.linalg import expm
from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp, Operator
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.synthesis import LieTrotter

H = SparsePauliOp.from_list([("X", 1.0), ("Z", 1.0)])   # non-commuting
t = 0.5
g = PauliEvolutionGate(H, time=t, synthesis=LieTrotter(reps=1))
qc = QuantumCircuit(1); qc.append(g, [0])

U_exact = expm(-1j * H.to_matrix() * t)

print(np.linalg.norm(Operator(qc).data - U_exact, 2))             # 0.0      <- synthesis silently ignored
print(np.linalg.norm(Operator(qc.decompose()).data - U_exact, 2)) # 0.23646  <- the real first-order Trotter error
print(np.linalg.norm(Operator(g).data - Operator(g.definition).data, 2))  # 0.23646  <- same gate, two unitaries

Output:

0.0
0.23646
0.23646

The value 0.23646 is the correct first-order Trotter error for H = X + Z at t = 0.5.

What should happen?

At minimum, requesting an approximate synthesis and then evaluating the gate's matrix should not be silent. Either Operator()/to_matrix() on a gate whose synthesis is an approximate (non matrix-exponential) formula should emit a UserWarning stating that the exact evolution is being returned and that the requested approximate synthesis is ignored, with a pointer to .decompose()/transpile; or the matrix evaluation should reflect the requested synthesis so that Operator(gate) and Operator(gate.definition) agree. The current situation, where the same gate object yields two different unitaries (the exact one via to_matrix() and the approximate one via .definition) with no signal connecting them, is the surprising part.

Any suggestions?

The lightest fix is a runtime UserWarning emitted from to_matrix()/__array__ when self.synthesis is an approximate formula, for instance "returning the exact evolution; the requested synthesis (LieTrotter) is ignored. Call .decompose() or transpile to obtain the synthesized circuit." A complementary docstring note placed directly on to_matrix() would also help. The mechanism, for reference, is that to_matrix() returns exp(-itH) and does not consult self.synthesis, while _define() sets self.definition = self.synthesis.synthesize(self); Operator(gate) resolves via to_matrix() and therefore bypasses the requested synthesis.