Repository: PennyLaneAI/pennylane

Issue:
[BUG] `sum` not working with hardware Hamiltonians (transmons and rydberg)
### Expected behavior

Working with hardware Hamiltonians, we add them together with similar syntax to `H = H_interaction + H_drive`. Often, we repeatedly add terms manually, e.g., `H_drive += qml.pulse.transmon_drive(amplitude, phase, freq, wires)`.

I would expect to be able to use `sum()` on a list of terms, such as:
`H_drive = sum([qml.pulse.transmon_drive(amplitude, phase, freq, q) for q in wires])`
Similar to how we can add other operators in PennyLane `sum([qml.PauliX(i) for i in range(wires)])`.

### Actual behavior

We can't do it.

### Additional information

I have tested this with a few `HardwareHamiltonian` instances using `qml.pulse.transmon_interaction`, `qml.pulse.transmon_drive`, `qml.pulse.rydberg_interaction` and `qml.pulse.rydberg_drive`.

We can manipulate all of them with stuff such as:
```python
H = qml.pulse.transmon_interaction(qubit_freqs, connection, couplings, wires)
for q in wires:
    H += qml.pulse.transmon_drive(amplitude, phase, freq, q)
```
Or even:
```python
H_int = qml.pulse.transmon_interaction(qubit_freqs, connection, couplings, wires)
drive_terms = [qml.pulse.transmon_drive(amplitude, phase, freq, q) for q in range(2)]
H_drive = drive_terms[0] + drive_terms[1]
H = H_int + H_drive
```

But we wouldn't be able to use `sum()` on any iterable containing `HardwareHamiltonian`s. I also tried changing the input types, e.g., setting `amplitude`, `phase`, and `freq` to be all callables, floats, or ints, or a combination of those.

### Source code

```shell
import pennylane as qml
sum([qml.pulse.transmon_drive(0., 0., 0., q) for q in range(2)])
```


### Tracebacks

```shell
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[1], line 2
      1 import pennylane as qml
----> 2 sum([qml.pulse.transmon_drive(0., 0., 0., q) for q in range(2)])

TypeError: unsupported operand type(s) for +: 'int' and 'HardwareHamiltonian'
```


### System information

```shell
Latest master

Name: PennyLane
Version: 0.32.0.dev0
Summary: PennyLane is a Python quantum machine learning library by Xanadu Inc.
Home-page: https://github.com/PennyLaneAI/pennylane
Author: 
Author-email: 
License: Apache License 2.0
Location: /home/borja.requena/.conda/envs/pl_dev/lib/python3.10/site-packages
Editable project location: /home/borja.requena/GitHub/pennylane
Requires: appdirs, autograd, autoray, cachetools, networkx, numpy, pennylane-lightning, requests, rustworkx, scipy, semantic-version, toml
Required-by: PennyLane-Lightning

Platform info:           Linux-5.10.102.1-microsoft-standard-WSL2-x86_64-with-glibc2.31
Python version:          3.10.11
Numpy version:           1.23.5
Scipy version:           1.10.0
Installed devices:
- default.gaussian (PennyLane-0.32.0.dev0)
- default.mixed (PennyLane-0.32.0.dev0)
- default.qubit (PennyLane-0.32.0.dev0)
- default.qubit.autograd (PennyLane-0.32.0.dev0)
- default.qubit.jax (PennyLane-0.32.0.dev0)
- default.qubit.tf (PennyLane-0.32.0.dev0)
- default.qubit.torch (PennyLane-0.32.0.dev0)
- default.qutrit (PennyLane-0.32.0.dev0)
- null.qubit (PennyLane-0.32.0.dev0)
- lightning.qubit (PennyLane-Lightning-0.31.0)
```


### Existing GitHub issues

- [X] I have searched existing GitHub issues to make sure the issue does not already exist.

Discussion:
dwierichs: There was a mistake in timing the `qml.matrix` evaluation: I did not pass parameters to that call, so that it did not really compute the matrices. I updated the times.

dwierichs: After rerunning timings with different methods (`timeit`, custom timer using `time.process_time`), I come to the conclusion that the fluctuations in the timings are too large to judge the impact of this PR properly :/ 
BorjaRequena: Nice! Then it is just as we expected with the constant term overhead only being noticeable for very small systems 😄 Also, it is only when we use `sum` or explicitly add a constant term so, if our goal was performance, we could always avoid this construction.

We could check if `other == 0` and return self to remove the term coming from `sum` (if it does not break jitability), as we discussed offline.
dwierichs: >We could check if other == 0 and return self to remove the term coming from sum (if it does not break jitability), as we discussed offline.

This actually seems to work fine, because the implementation is in particular for `int` and `float`, so that tracers are not "accidentally" supported anyways. :) 

I added tests.
dwierichs: Yep, this is equivalent to manually initializing the first term and then adding the others, so it never hits the newly added logic :)

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 93c39b88d..ef846119c 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -13,6 +13,10 @@
 
 <h3>Improvements 🛠</h3>
 
+* `HardwareHamiltonian`s can now be summed with `int` or `float`.
+  A sequence of `HardwareHamiltonian`s can now be summed via the builtin `sum`.
+  [(#4343)](https://github.com/PennyLaneAI/pennylane/pull/4343)
+
 * All `Operator` objects now define `Operator._flatten` and `Operator._unflatten` methods that separate
   trainable from untrainable components. These methods will be used in serialization and pytree registration.
   Custom operations may need an update to ensure compatibility with new PennyLane features.
diff --git a/pennylane/pulse/hardware_hamiltonian.py b/pennylane/pulse/hardware_hamiltonian.py
index eb20eb782..f3c54d391 100644
--- a/pennylane/pulse/hardware_hamiltonian.py
+++ b/pennylane/pulse/hardware_hamiltonian.py
@@ -26,9 +26,6 @@ from pennylane.ops.qubit.hamiltonian import Hamiltonian
 
 from .parametrized_hamiltonian import ParametrizedHamiltonian
 
-# from .transmon import TransmonSettings
-# from .rydberg import RydbergSettings
-
 
 def drive(amplitude, phase, wires):
     r"""Returns a :class:`ParametrizedHamiltonian` representing the action of a driving electromagnetic
@@ -325,11 +322,14 @@ class HardwareHamiltonian(ParametrizedHamiltonian):
     def __repr__(self):
         return f"HardwareHamiltonian: terms={qml.math.shape(self.coeffs)[0]}"
 
-    def __add__(self, other):
+    def __add__(self, other):  # pylint: disable=too-many-return-statements
         if isinstance(other, HardwareHamiltonian):
             if not self.reorder_fn == other.reorder_fn:
                 raise ValueError(
-                    f"Cannot add two HardwareHamiltonians with different reorder functions. Received reorder_fns {self.reorder_fn} and {other.reorder_fn}. This is likely due to an attempt to add hardware compatible Hamiltonians for different target systems."
+                    "Cannot add two HardwareHamiltonians with different reorder functions. "
+                    f"Received reorder_fns {self.reorder_fn} and {other.reorder_fn}. This is "
+                    "likely due to an attempt to add hardware compatible Hamiltonians for "
+                    "different target systems."
                 )
             if self.settings is None and other.settings is None:
                 new_settings = None
@@ -374,6 +374,19 @@ class HardwareHamiltonian(ParametrizedHamiltonian):
                 new_coeffs, new_ops, reorder_fn=self.reorder_fn, settings=settings, pulses=pulses
             )
 
+        if isinstance(other, (int, float)):
+            if other in (0, 0.0):
+                return HardwareHamiltonian(
+                    coeffs, ops, reorder_fn=self.reorder_fn, settings=settings, pulses=pulses
+                )
+            new_coeffs = coeffs + [other]
+            with qml.queuing.QueuingManager.stop_recording():
+                new_ops = ops + [qml.Identity(self.wires[0])]
+
+            return HardwareHamiltonian(
+                new_coeffs, new_ops, reorder_fn=self.reorder_fn, settings=settings, pulses=pulses
+            )
+
         return NotImplemented
 
     def __radd__(self, other):


Test patch:
diff --git a/tests/pulse/test_hardware_hamiltonian.py b/tests/pulse/test_hardware_hamiltonian.py
index 315aedf92..566c3c22e 100644
--- a/tests/pulse/test_hardware_hamiltonian.py
+++ b/tests/pulse/test_hardware_hamiltonian.py
@@ -166,6 +166,55 @@ class TestHardwareHamiltonian:
         assert res2.ops_parametrized == []
         assert res2.wires == qml.wires.Wires(h_wires)
 
+    @pytest.mark.parametrize("scalars", [(0.2, 1), (0.9, 0.2), (1, 5), (3, 0.5)])
+    def test_add_scalar(self, scalars):
+        """Test that adding a scalar/number to a `HardwareHamiltonian` works as expected."""
+        coeffs = [2, 3]
+        ops = [qml.PauliZ(2), qml.PauliX(0)]
+
+        H = HardwareHamiltonian(
+            coeffs=coeffs,
+            observables=ops,
+            pulses=[HardwarePulse(qml.pulse.constant, 6, 7, 8)],
+        )
+        orig_matrix = qml.matrix(H([0.3], 0.1))
+        H1 = H + scalars[0]
+        assert len(H1.ops) == len(H1.coeffs) == 3
+        assert qml.equal(H1.ops[-1], qml.Identity(2))
+        assert qml.math.allclose(qml.matrix(H1([0.3], 0.1)), orig_matrix + np.eye(4) * scalars[0])
+        H2 = scalars[1] + H1
+        assert len(H2.ops) == len(H2.coeffs) == 4
+        assert qml.equal(H2.ops[-2], qml.Identity(2))
+        assert qml.equal(H2.ops[-1], qml.Identity(2))
+        assert qml.math.allclose(
+            qml.matrix(H2([0.3], 0.1)), orig_matrix + np.eye(4) * (scalars[0] + scalars[1])
+        )
+        H += scalars[0]
+        assert len(H.ops) == len(H.coeffs) == 3
+        assert qml.equal(H.ops[-1], qml.Identity(2))
+        assert qml.math.allclose(qml.matrix(H([0.3], 0.1)), orig_matrix + np.eye(4) * scalars[0])
+
+    def test_add_zero(self):
+        """Test that adding an int or a float that is zero to a `HardwareHamiltonian`
+        returns an unchanged copy."""
+        coeffs = [2, 3]
+        ops = [qml.PauliZ(2), qml.PauliX(0)]
+
+        H = HardwareHamiltonian(
+            coeffs=coeffs,
+            observables=ops,
+            pulses=[HardwarePulse(qml.pulse.constant, 6, 7, 8)],
+        )
+        orig_matrix = qml.matrix(H([0.3], 0.1))
+        H1 = H + 0
+        assert len(H1.ops) == len(H1.coeffs) == 2
+        assert qml.math.allclose(qml.matrix(H1([0.3], 0.1)), orig_matrix)
+        assert H1 is not H
+        H2 = 0.0 + H
+        assert len(H2.ops) == len(H2.coeffs) == 2
+        assert qml.math.allclose(qml.matrix(H2([0.3], 0.1)), orig_matrix)
+        assert H2 is not H
+
     @pytest.mark.xfail
     def test_add_raises_warning(self):
         """Test that an error is raised when adding two HardwareHamiltonians where one Hamiltonian
@@ -296,11 +345,18 @@ class TestInteractionWithOperators:
         assert qml.equal(new_pH.H_fixed(), qml.s_prod(1, op))
         new_pH(params, 2)  # confirm calling does not raise error
 
+    def test_adding_scalar_does_not_queue_id(self):
+        """Test that no additional Identity operation is queued when adding a scalar."""
+        R = drive(amplitude=f1, phase=0, wires=[0, 1])
+        with qml.queuing.AnnotatedQueue() as q:
+            R += 3
+        assert len(q) == 0
+
     def test_unknown_type_raises_error(self):
         """Test that adding an invalid object to a HardwareHamiltonian raises an error."""
         R = drive(amplitude=f1, phase=0, wires=[0, 1])
         with pytest.raises(TypeError, match="unsupported operand type"):
-            R += 3
+            R += 3j
 
 
 class TestDrive: