def hamiltonian_jac_element(prep, error, observable):
    """
    Computes the Jacobian matrix element for a Hamiltonian error: how the
    expectation value of `observable` in state `prep` changes due to
    Hamiltonian `error`.

    Parameters
    ----------
    prep : NQPauliState
        The state that is prepared.

    error : NQPauliOp
        The type (as a pauli operator) of Hamiltonian error.

    observable : NQPauliOp
        The observable whose expectation value is measured.
        Note: giving a NQPauliState will be treated as an
        N-qubit Pauli operator whose sign is the product of
        the signs of the NQPauliState's per-qubit basis signs.

    Returns
    -------
    float
    """
    # dp/deps where p = eps * i * Tr(Obs (Err*rho - rho*Err)) = eps * i * ( Tr(Obs Err rho) - Tr(Obs rho Err))
    #                 = eps * i * Tr([Obs,Err] * rho)  so dp/deps just drops eps factor
    com = error.icommutator_over_2(observable)
    return 0 if (com is None) else com.statedot(prep)