IBM circuit functions returning erroneous results on simple circuits

[nonhermitian]

Steps to reproduce the problem

Here is the circuit and operator:

op = SparsePauliOp('ZIIIII')

g = Parameter('g')
g_vals = np.linspace(0, np.pi, 11)

qc = QuantumCircuit(6)
qc.h(0)
qc.cx(0, 5)
qc.h(1)
qc.cx(1, 4)
qc.h(2)
qc.cx(2, 3)
qc.rx(-np.pi/2, 0)
qc.rx(-np.pi/2, 1)
qc.rx(-np.pi/2, 2)
qc.rz(-0.0566794, 0)
qc.rz(-0.01039906, 1)
qc.rz(-0.0632158, 2)
qc.rzz(-np.pi/2, 0, 1)
qc.rzz(-np.pi/2, 1, 2)
qc.rx(-np.pi/2, 0)
qc.rx(-np.pi/2, 1)
qc.rx(-np.pi/2, 2)
qc.rz(-0.0566794, 0)
qc.rz(-0.01039906, 1)
qc.rz(-0.0632158, 2)
qc.rzz(-np.pi/2, 0, 1)
qc.rzz(-np.pi/2, 1, 2)
qc.rx(-np.pi/2, 0)
qc.rx(-np.pi/2, 1)
qc.rx(-np.pi/2, 2)
qc.rz(-0.0566794, 0)
qc.rz(-0.01039906, 1)
qc.rz(-0.0632158, 2)
qc.rzz(-np.pi/2, 0, 1)
qc.rzz(-np.pi/2, 1, 2)
qc.reset(0)
qc.rz(0.0566794, 0)
qc.rz(0.01039906, 1)
qc.rz(0.0632158, 2)
qc.rzz(np.pi/2, 0, 1)
qc.rzz(np.pi/2, 1, 2)
qc.rx(np.pi/2, 0)
qc.rx(np.pi/2, 1)
qc.rx(np.pi/2, 2)
qc.rz(0.0566794, 0)
qc.rz(0.01039906, 1)
qc.rz(0.0632158, 2)
qc.rzz(np.pi/2, 0, 1)
qc.rzz(np.pi/2, 1, 2)
qc.rx(np.pi/2, 0)
qc.rx(np.pi/2, 1)
qc.rx(np.pi/2, 2)
qc.rz(0.0566794, 0)
qc.rz(0.01039906, 1)
qc.rz(0.0632158, 2)
qc.rzz(np.pi/2, 0, 1)
qc.rzz(np.pi/2, 1, 2)
qc.rx(np.pi/2, 0)
qc.rx(np.pi/2, 1)
qc.rx(np.pi/2, 2)
# Here are the two RZZ gates that are parameterized
qc.rzz(g, 1, 4)
qc.rzz(g, 2, 3)
# -------------------------------------------------
qc.rx(np.pi/2, 5)
qc.rx(np.pi/2, 4)
qc.rx(np.pi/2, 3)
qc.rz(0.0566794, 5)
qc.rz(0.01039906, 4)
qc.rz(0.0632158, 3)
qc.rzz(np.pi/2, 5, 4)
qc.rzz(np.pi/2, 4, 3)
qc.rx(np.pi/2, 5)
qc.rx(np.pi/2, 4)
qc.rx(np.pi/2, 3)
qc.rz(0.0566794, 5)
qc.rz(0.01039906, 4)
qc.rz(0.0632158, 3)
qc.rzz(np.pi/2, 5, 4)
qc.rzz(np.pi/2, 4, 3)
qc.rx(np.pi/2, 5)
qc.rx(np.pi/2, 4)
qc.rx(np.pi/2, 3)
qc.rz(0.0566794, 5)
qc.rz(0.01039906, 4)
qc.rz(0.0632158, 3)
qc.rzz(np.pi/2, 5, 4)
qc.rzz(np.pi/2, 4, 3)

What is the current behavior?

As shown below, the IBM circuit functions return worse results and calling the Runtime directly, irrespective of mitigation level. This is not a small difference, but is systematically worse across the board. The default circuit function should be closely matching a call to the Runtime with DD activation (here I use XY4), but this is far from the case. This same behavior has been observed on 3 different systems. As such, it strongly appears to be a bug in the circuit functions themselves.

What is the expected behavior?

The default circuit function should be nearly the same as a direct call to Runtime with DD activated.