Repository: PennyLaneAI/pennylane (python)
Issue: [BUG] fourier.qnode_spectrum() returns different results depending on how the generator of a gate is defined
### Expected behavior

The output of `fourier.qnode_spectrum()` should be the same for two QNodes that differ only in that the generator of one of the gates is defined in a different way, which however represents the same operator.

### Actual behavior

Strangely enough I get diffrent spectra depending on whether `generator(self)` returns a `pennylane.ops.op_math.linear_combination.LinearCombination` or a `pennylane.ops.op_math.sum.Sum` representing the same operator.

See below for a minimal non-working example, which for me produces the output:

```
orig_generator    <class 'pennylane.ops.op_math.linear_combination.LinearCombination'> 0.25 * (X(0) @ Y(1)) + -0.25 * (Y(0) @ X(1))
[[0.+0.j  0.+0.j  0.+0.j  0.+0.j ]
 [0.+0.j  0.+0.j  0.+0.5j 0.+0.j ]
 [0.+0.j  0.-0.5j 0.+0.j  0.+0.j ]
 [0.+0.j  0.+0.j  0.+0.j  0.+0.j ]]
patched_generator <class 'pennylane.ops.op_math.sum.Sum'> (-0.25 * Y(0)) @ X(1) + (0.25 * X(0)) @ Y(1)
[[0.+0.j  0.+0.j  0.+0.j  0.+0.j ]
 [0.+0.j  0.+0.j  0.+0.5j 0.+0.j ]
 [0.+0.j  0.-0.5j 0.+0.j  0.+0.j ]
 [0.+0.j  0.+0.j  0.+0.j  0.+0.j ]]
<class 'pennylane.ops.qubit.qchem_ops.SingleExcitation'>
{'params': {(0,): [-1.0, -0.5, 0.0, 0.5, 1.0]}}
<class '__main__.MySingleExcitation'>
{'params': {(0,): [-1.0, 0.0, 1.0]}}
```

### Additional information

This used to work with v0.35.1.

### Source code

```shell
from pennylane import numpy as np
import pennylane as qml

dev = qml.device("default.qubit", wires=2, shots=None)

class MySingleExcitation(qml.SingleExcitation):
    # the definition from qml.SingleExcitation is
    # def generator2(self):
    #     w1, w2 = self.wires
    #     return qml.Hamiltonian([0.25, -0.25], [qml.X(w1) @ qml.Y(w2), qml.Y(w1) @ qml.X(w2)])
    # which should be exactly the same as:
    def generator(self):
        w1, w2 = self.wires
        return -0.25 * qml.PauliY(w1) @ qml.PauliX(w2) + 0.25 * qml.PauliX(w1) @ qml.PauliY(w2)

orig_generator = qml.SingleExcitation(0.1, wires=range(2)).generator()
patched_generator = MySingleExcitation(0.1, wires=range(2)).generator()

print("orig_generator   ", type(orig_generator), orig_generator)
print(qml.matrix(orig_generator, wire_order=range(2)))
print("patched_generator", type(patched_generator), patched_generator)
print(qml.matrix(patched_generator, wire_order=range(2)))

# check that they are the same generator
assert np.allclose(
    qml.matrix(orig_generator, wire_order=range(2)),
    qml.matrix(patched_generator, wire_order=range(2)),
)

for gate in [qml.SingleExcitation, MySingleExcitation]:

    print(gate)

    @qml.qnode(dev)
    def qnode(params):

        qml.PauliX(wires=0)
        qml.PauliX(wires=1)
        qml.OrbitalRotation(0.4, wires=range(4), id="OrbitalRotation")
        gate(params[0], wires=range(2), id="ReFermionEx")

        return qml.expval(qml.PauliZ(wires=0))

    params = np.array([0.1], requires_grad=True)
    spectrum_fn = qml.fourier.qnode_spectrum(qnode)
    spectra = spectrum_fn(params)
    print(spectra)
```

### Tracebacks

```shell

```

### System information

```shell
code was run with v0.40.0 I will provide this later.
```

### Existing GitHub issues

- [x] I have searched existing GitHub issues to make sure the issue does not already exist.
Pull request description: Fix gen_is_multi_term_hamiltonian
**Context:**

`qml.operation.gen_is_multiterm_hamiltonian` was never properly updated for new operator arithmetic.

**Description of the Change:**

Allow is to work with general generators.

**Benefits:**

The `fourier.qnode_spectrum` works for general generators.

**Possible Drawbacks:**

**Related GitHub Issues:**
Fixes #7116 

**Related Shortcut Stories:**
[sc-86870]


Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 96b025e93..ad96e177d 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -714,6 +714,12 @@
 
 <h3>Bug fixes 🐛</h3>
 
+* Fixes a bug that caused the output of `qml.fourier.qnode_spectrum()` to
+  differ depending if equivalent gate generators are defined using
+  different PennyLane operators. This was resolved by updating
+  `qml.operation.gen_is_multi_term_hamiltonian` to work with more complicated generators.
+  [(#7121)])(https://github.com/PennyLaneAI/pennylane/pull/7121)
+
 * Modulo operator calls on MCMs now correctly offload to the autoray-backed `qml.math.mod` dispatch.
   [(#7085)](https://github.com/PennyLaneAI/pennylane/pull/7085)
 
diff --git a/pennylane/operation.py b/pennylane/operation.py
index 36f5448b0..f2703631f 100644
--- a/pennylane/operation.py
+++ b/pennylane/operation.py
@@ -2416,14 +2416,15 @@ def defines_diagonalizing_gates(obj):
 def gen_is_multi_term_hamiltonian(obj):
     """Returns ``True`` if an operator has a generator defined and it is a Hamiltonian
     with more than one term."""
-
+    if not isinstance(obj, Operator) or not obj.has_generator:
+        return False
     try:
-        o = obj.generator()
-    except (AttributeError, OperatorPropertyUndefined, GeneratorUndefinedError):
+        generator = obj.generator()
+        _, ops = generator.terms()  # len(coeffs) can be weird sometimes
+        return len(ops) > 1
+    except TermsUndefinedError:
         return False
 
-    return isinstance(o, qml.ops.LinearCombination) and len(o.coeffs) > 1
-
 
 def __getattr__(name):
     """To facilitate StatePrep rename"""


Test patch:
diff --git a/tests/test_operation.py b/tests/test_operation.py
index 3d464d519..138e70b8b 100644
--- a/tests/test_operation.py
+++ b/tests/test_operation.py
@@ -1611,6 +1611,20 @@ class TestCriteria:
         assert not qml.operation.gen_is_multi_term_hamiltonian(self.rot)
         assert not qml.operation.gen_is_multi_term_hamiltonian(self.exp)
 
+        class SProdGen(Operator):
+
+            def generator(self):
+                return 2.0 * (qml.X(0) + qml.Y(0))
+
+        assert qml.operation.gen_is_multi_term_hamiltonian(SProdGen(wires=0))
+
+        class SumGen(Operator):
+
+            def generator(self):
+                return qml.X(0) + qml.Y(1)
+
+        assert qml.operation.gen_is_multi_term_hamiltonian(SumGen(wires=0))
+
     def test_has_multipar(self):
         """Test has_multipar criterion."""
         assert not qml.operation.has_multipar(self.rx)