Remove usage of qml.Hamiltonian from hamiltonian.py (#5429)

`hamiltonian.py` defines the legacy `Hamiltonian`. It should not be
returning `LinearCombination` in its own methods even with new op math
enabled.

Co-authored-by: qottmann <korbinian.kottmann@gmail.com>

pennylane/ops/qubit/hamiltonian.py

@@ -67,6 +67,7 @@ class Hamiltonian(Observable):
     The Hamiltonian is represented as a linear combination of other operators, e.g.,
     :math:`\sum_{k=0}^{N-1} c_k O_k`, where the :math:`c_k` are trainable parameters.
 
+
     Args:
         coeffs (tensor_like): coefficients of the Hamiltonian expression
         observables (Iterable[Observable]): observables in the Hamiltonian expression, of same length as coeffs
@@ -675,12 +676,12 @@ def __matmul__(self, H):
             coeffs = qml.math.kron(coeffs1, coeffs2)
             ops_list = itertools.product(ops1, ops2)
             terms = [qml.operation.Tensor(t[0], t[1]) for t in ops_list]
-            return qml.Hamiltonian(coeffs, terms, simplify=True)
+            return Hamiltonian(coeffs, terms, simplify=True)
 
         if isinstance(H, (Tensor, Observable)):
             terms = [op @ copy(H) for op in ops1]
 
-            return qml.Hamiltonian(coeffs1, terms, simplify=True)
+            return Hamiltonian(coeffs1, terms, simplify=True)
 
         return NotImplemented
 
@@ -697,7 +698,7 @@ def __rmatmul__(self, H):
         if isinstance(H, (Tensor, Observable)):
             terms = [copy(H) @ op for op in ops1]
 
-            return qml.Hamiltonian(coeffs1, terms, simplify=True)
+            return Hamiltonian(coeffs1, terms, simplify=True)
 
         return NotImplemented
 
@@ -712,14 +713,14 @@ def __add__(self, H):
         if isinstance(H, Hamiltonian):
             coeffs = qml.math.concatenate([self_coeffs, copy(H.coeffs)], axis=0)
             ops.extend(H.ops.copy())
-            return qml.Hamiltonian(coeffs, ops, simplify=True)
+            return Hamiltonian(coeffs, ops, simplify=True)
 
         if isinstance(H, (Tensor, Observable)):
             coeffs = qml.math.concatenate(
                 [self_coeffs, qml.math.cast_like([1.0], self_coeffs)], axis=0
             )
             ops.append(H)
-            return qml.Hamiltonian(coeffs, ops, simplify=True)
+            return Hamiltonian(coeffs, ops, simplify=True)
 
         return NotImplemented
 
@@ -730,7 +731,7 @@ def __mul__(self, a):
         if isinstance(a, (int, float)):
             self_coeffs = copy(self.coeffs)
             coeffs = qml.math.multiply(a, self_coeffs)
-            return qml.Hamiltonian(coeffs, self.ops.copy())
+            return Hamiltonian(coeffs, self.ops.copy())
 
         return NotImplemented