Repository: PennyLaneAI/pennylane
Pull request #3749 (2023-02-10 15:19:59)

Rydberg Hamiltonian
This PR contains the following additions/changes:
- Add `RydbergHamiltonian` class, which is an internal class (should NOT be user-facing!) that contains the representation of a Rydberg Hamiltonian together with all the information needed to execute this Hamiltonian in hardware.
- Add `rydberg_interaction` and `rydberg_transition` functions that return `RydbergHamiltonian` instances.
- Add `RydbergPulse` dataclass, which is an internal class that is used to group all the information of a single pulse. This class is not strictly needed, but it helps structure all the pulse data.
- Use `TensorLike` in `Operation` dunder methods. Add `__truediv__` dunder method.
- Add `ParametrizedHamiltonian.__radd__` dunder method to make sure coefficients are ordered correctly.

github-actions[bot]: Hello. You may have forgotten to update the changelog!
Please edit [`doc/releases/changelog-dev.md`](/PennyLaneAI/pennylane/blob/master/doc/releases/changelog-dev.md) with:
* A one-to-two sentence description of the change. You may include a small working example for new features.
* A link back to this PR.
* Your name (or GitHub username) in the contributors section.
codecov[bot]: # [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3749?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> Merging [#3749](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3749?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (d3e58e6) into [master](https://codecov.io/gh/PennyLaneAI/pennylane/commit/0f48aa6f5ac34a624e5de2b8b09d47a9a486cadd?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (0f48aa6) will **decrease** coverage by `0.04%`.
> The diff coverage is `43.75%`.

```diff
@@            Coverage Diff             @@
##           master    #3749      +/-   ##
==========================================
- Coverage   99.82%   99.79%   -0.04%     
==========================================
  Files         329      330       +1     
  Lines       28856    28871      +15     
==========================================
+ Hits        28805    28811       +6     
- Misses         51       60       +9     
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3749?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Coverage Δ | |
|---|---|---|
| [pennylane/pulse/rydberg.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3749?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3B1bHNlL3J5ZGJlcmcucHk=) | `35.71% <35.71%> (ø)` | |
| [pennylane/pulse/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3749?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3B1bHNlL19faW5pdF9fLnB5) | `100.00% <100.00%> (ø)` | |

Help us with your feedback. Take ten seconds to tell us [how you rate us](https://about.codecov.io/nps?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Have a feature suggestion? [Share it here.](https://app.codecov.io/gh/feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)

lillian542: Nice work 💯 Cool that there's already a working prototype for this! 🚀 

I had a few comments/questions, mostly about things we might consider as this gets more fleshed out, and making sure it aligns where reasonable with the QuEra hardware.
AlbertMitjans: [sc-32760]
AlbertMitjans: This class should allow the following:
1. Allow local drives to single atoms in the array.
2. Allow global drives to all the atoms.
3. Allow translation into `braket` notation:
  - The `braket` plugin loops over the operations in a tape and translates them. Thus we should queue a `RydbergMachine` to be able to extract all the relevant information of the pulses.

I believe the best way to do this is to make `RydbergMachine` a [class decorator](https://refactoring.guru/design-patterns/decorator) that wraps a `ParametrizedHamiltonian`. For each call to `drive` or `global_drive`, the class computes the term of that driving and adds it to the wrapped hamiltonian.

To be able to queue `RydbergMachine`s, we need to add support for them in `ParametrizedEvolution`. This is done by adding the `__call__`, `ops` and `coeffs` methods (that internally call the methods of the wrapped `ParametrizedHamiltonian`.

Let me know what you think!


AlbertMitjans: [sc-32760]
lillian542: Here are some things I encounter when I try to play around with it. Could definitely be some of these are user-error issues on my end though:

1. In the `local_drive`, it doesn't work if phase is a callable 
2. Possibly the intended behaviour, but if I call `local_drive` more than once, which I might be inclined to do in order to first set the drive on qubits 1, 5, and 9 that all have the same drive, and then on 2, 6, and 10, etc., it fails.
3. In the `global_drive`, it also doesn't work if `phase` is a callable
4. It's unintuitive to me that it's not possible to update the `global_drive` term or add onto the `local_drive`, term since they are just being included in the class. I think it might be better if these functions just returned a `global_drive` Hamiltonian, or a `local_drive`, or a `rydberg_interaction` based on the current positions. But this is very much a UI question and I would like to hear other people's thoughts on this. If it stays as is, I think it should be named `RydbergHamiltonian` instead of `RydbergMachine`, to make it clear what it is. But I think this could be more modular rather than a single object.
lillian542: One thing that I think is a bit difficult, and that could be easier if this were more modular, is keeping track of the order of parameters. This object is hard to examine to understand in what order things come, and the order will depend on in what order you updated `global_drive`, and potentially sequential calls to `local_drive`. 

If you could create `qubit_1_drive = rydberg_local_drive(amp, detuning, phase)`, `qubit_2_drive...` etc, and then add them together and make a params list that was `[qubit1_params, qubit2_params]`, etc. I think it would be easier to keep track of what is going on.
lillian542: Nice to see more of the parts coming together! 💯 

I think there are still some discussions to be had about UI, but this might need to wait until more people have a chance to look at it/think about it.
AlbertMitjans: > I'm not great with decorators, and I've only used them for functions. Can you give an example of how to you would use the class to decorate a ParametrizedHamiltonian? Or is that already what is being done when the class contains an attribute that is a ParametrizedHamiltonian and can update and return it?

A decorator is just a wrapper of an object. However, I misspoke before, because `RydbergMachine` is not a decorator (because it doesn't have a `ParametrizedHamiltonian` argument in the `__init__`). It is a class that contains a `ParametrizedHamiltonian`, which is a bit different (it doesn't wrap a `ParametrizedHamiltonian`, but it actually defines one).

> I.e. is the way to use this still to create a RydbergMachine instance, or do I create a ParametrizedHamiltonian and somehow decorate it with the RydbergMachine class?

No, you cannot decorate a `ParametrizedHamiltonian` with the `RydbergMachine` class. Sorry for the confusion.
AlbertMitjans: > 1. In the local_drive, it doesn't work if phase is a callable
2. In the global_drive, it also doesn't work if phase is a callable

Right now this is not supported ([see comment](https://github.com/PennyLaneAI/pennylane/pull/3749#discussion_r1108200226)).

> 3. It's unintuitive to me that it's not possible to update the global_drive term. I think it might be better if these functions just returned a global_drive Hamiltonian, or a local_drive, or a rydberg_interaction based on the current positions. But this is very much a UI question and I would like to hear other people's thoughts on this. If it stays as is, I think it should be named RydbergHamiltonian instead of RydbergMachine, to make it clear what it is. But I think this could be more modular rather than a single object.

Yes, this definitely needs discussion.
AlbertMitjans: > One thing that I think is a bit difficult, and that could be easier if this were more modular, is keeping track of the order of parameters.

Why do we want to keep track of the order of parameters?
lillian542: > Why do we want to keep track of the order of parameters?

We want **users** to keep track of the order of parameters. If I do `qml.evolve(H_drift + H1 + H2 + H3 + H_global)`, then I know to pass parameters in the order `[params_qubit_1, params_qubit2, params_qubit_3, params_global]`. If it's all 'hidden' inside the Hamiltonian object, then as a user, I don't know if it wants the global params first, or the local. Does it sort the local terms by wire order, or by the order I added them it? etc. It will require more attention to the docstrings and more of an understanding of the hidden structure than if I can just create the different terms and pass them into `evolve`. 

Basically, I'll have to learn more new things to use the RydbergHamiltonian, whereas if we go with something more modular and people already know how ParametrizedHamiltonian works, I think it will be a simpler and more intuitive switch. It will also feel like a more separate thing from the ParametrizedHamiltonian if its less modular. 

If making something that is structurally quite different from how we expect users to work with `ParametrizedHamiltonian` is what we decide we need to do for hardware, I think there should still be, as @Qottmann talked about yesterday, a 'philosophically hardware compatible' thing that can't actually upload to hardware, but describes the behaviour of hardware, and works in a way that is more in line with how `ParametrizedHamiltonian` has been made to work.
AlbertMitjans: > I think there should still be, as @Qottmann talked about yesterday, a 'philosophically hardware compatible' thing that can't actually upload to hardware

I don't like the idea of having 2 different classes/objects for executing the same thing on a simulator and on hardware. I don't think this follows PL standards.

> If it's all 'hidden' inside the Hamiltonian object, then as a user, I don't know if it wants the global params first, or the local.

It will depend on the order of the `local_drive` and the `global_drive`. If having 2 drive functions is confusing, we could always use a generic `drive` function, and use a specific string/object (such as `"global"`) to define that the pulse is global (just like `Pulser` does). Then I don't think there will be any confusion: the order of the parameters corresponds to the order of the pulses in the `drive` function.
Qottmann: > I don't like the idea of having 2 different classes/objects for executing the same thing on a simulator and on hardware. I don't think this follows PL standards.

I agree that ideally these things are all unified. I still know too little of what hardware vendors and/or aws braket require to say what is the best solution to do that. For now the purpose of this feature here is to have a convenience for the user to generate ParametrizedHamiltonians for pulses that resemble realistic interactions. Just to explain the motivation here: I was testing a lot of different things for VQE using the transmon Hamiltonian, and for every script I had to copy-paste the parameters and construction script of the pulse Hamiltonian. It would be nice to have that accessible with good coefficients in a more easy way. At the same time I would want these convenience features to be enabling and not restricting, so users can re-use or even re-purpose parts of the code.
lillian542: I also agree that ideally we should unify these with how we send signals to hardware, but if we can't, we shouldn't avoid having a more flexible option so that there is only one way - the limited, hardware compatible way - to define a RydbergHamiltonian. I don't think this interface looks or feels very much like what we are asking people to do to create a `ParametrizedHamiltonian`.

> It will depend on the order of the local_drive and the global_drive.

Of course it does - but what *is* the order of these things, for someone who doesn't know? I would honestly expect, if I called `R.add_local_drive`, then `R.add_global_drive`, then try that out, decide to modify, and call `R.add_local_drive` again to add a drive on an additional qubit (or at a later time, or even an overlapping pulse), that I would need to call `R([local_params1, local_params2, global_params]`), because I would assume these are separate categories. But that's not what it does.

This is what I feel is hidden. I had to go in, read what the class `__call__` returns, find `driving_interaction`, and check how the `local_drive` and `global_drive` functions modify `driving_interaction` to be sure I knew how this would behave. It's not that obvious.
AlbertMitjans: [sc-34211]
lillian542: Since the `rydberg_drive` and `rydberg_interaction` terms are closely linked, should be consider an example in each (maybe under Usage Details) that combines them? 

Also for `rydberg_drive` we might consider adding a Usage Details section with a more detailed example that adds a global drive, plus some local drive to a subset of atoms (this could also be the example containing the interaction term).
lillian542: We've addressed distance being in micrometers, and energy/frequency in MHz - where do we specify that times are also expected to be in μs? 

This is relevant when passing to `qml.evolve` (but can't really be documented there, because it's specific to the RydbergHamiltonian), and when defining the functions describing phase/amplitude/detuning.
Qottmann: Good point! I would say in the same line / place where we mention the microseconds and mega hertz. OR a Returns statement that mentions the callabes are expected to have the signature (p, t) and t is expected in mus 
mudit2812: I'm a bit confused about why we can't say that the unit of time is microseconds in `qml.evolve` or `qml.ParametrizedEvolution`. If `RydbergHamiltonian` is a subclass of `ParametrizedHamiltonian`, then I would imagine the units of time that are used should be the same for both. Is that not the case?
lillian542: > If RydbergHamiltonian is a subclass of ParametrizedHamiltonian, then I would imagine the units of time that are used should be the same for both. Is that not the case?

The scaling factors are different for different systems. For Rydberg, the relevant frequency scale is MHz, and the timescale um. For transmons, they will be GHz and ns. For a general ParametrizedHamiltonian, people can do whatever makes sense for the physical system they are describing wrt to 'hiding' large scaling factors in the units - they'll have to figure out what makes sense.

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 309646a1d..06ebf21a1 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -8,6 +8,18 @@
   state vector and a number of shots.
   [(#3720)](https://github.com/PennyLaneAI/pennylane/pull/3720)
 
+* Added the needed functions and classes to simulate an ensemble of Rydberg atoms:
+  * A new internal `RydbergHamiltonian` class is added, which contains the Hamiltonian of an ensemble of
+    Rydberg atoms.
+  * A new user-facing `rydberg_interaction` function is added, which returns a `RydbergHamiltonian` containing
+    the Hamiltonian of the interaction of all the Rydberg atoms.
+  * A new user-facing `rydberg_drive` function is added, which returns a `RydbergHamiltonian` containing
+    the Hamiltonian of the interaction between a driving laser field and a group of atoms.
+  [(#3749)](https://github.com/PennyLaneAI/pennylane/pull/3749)
+
+* Added `Operation.__truediv__` dunder method to be able to divide operators.
+  [(#3749)](https://github.com/PennyLaneAI/pennylane/pull/3749)
+
 * The `simulate` function added to `devices/qubit` now supports measuring expectation values of large observables such as
   `qml.Hamiltonian`, `qml.SparseHamiltonian`, `qml.Sum`.
   [(#3759)](https://github.com/PennyLaneAI/pennylane/pull/3759)
@@ -29,6 +41,9 @@
 * `AdaptiveOptimizer` is updated to use non-default user-defined qnode arguments.
   [(#3765)](https://github.com/PennyLaneAI/pennylane/pull/3765)
 
+* Use `TensorLike` type in `Operator` dunder methods.
+  [(#3749)](https://github.com/PennyLaneAI/pennylane/pull/3749)
+
 * The `apply_operation` function added to `devices/qubit` now supports broadcasting.
   [(#3852)](https://github.com/PennyLaneAI/pennylane/pull/3852)
 
@@ -42,6 +57,7 @@
   to be set as trainable), and therefore the derivatives are computed more efficiently.
   [(#3697)](https://github.com/PennyLaneAI/pennylane/pull/3697)
 
+
 <h3>Breaking changes</h3>
 
 * An operation that implements a custom `generator` method, but does not always return a valid generator, also has
@@ -70,6 +86,10 @@
 
 <h3>Bug fixes</h3>
 
+* Fixed bug where the coefficients where not ordered correctly when summing a `ParametrizedHamiltonian`
+  with other operators.
+  [(#3749)](https://github.com/PennyLaneAI/pennylane/pull/3749)
+
 * The metric tensor transform is fully compatible with Jax and therefore users can provide multiple parameters.
   [(#3847)](https://github.com/PennyLaneAI/pennylane/pull/3847)
 
@@ -95,6 +115,7 @@ Lillian M. A. Frederiksen
 Soran Jahangiri
 Christina Lee
 Vincent Michaud-Rioux
+Albert Mitjans
 Romain Moyard
 Mudit Pandey
 Matthew Silverman
diff --git a/pennylane/operation.py b/pennylane/operation.py
index 553a8983c..ae045c391 100644
--- a/pennylane/operation.py
+++ b/pennylane/operation.py
@@ -216,12 +216,10 @@ import abc
 import copy
 import functools
 import itertools
-import numbers
 import warnings
 from enum import IntEnum
 from typing import List
 
-import autoray
 import numpy as np
 from numpy.linalg import multi_dot
 from scipy.sparse import coo_matrix, eye, kron
@@ -229,6 +227,7 @@ from scipy.sparse import coo_matrix, eye, kron
 import pennylane as qml
 from pennylane.math import expand_matrix
 from pennylane.queuing import QueuingManager
+from pennylane.typing import TensorLike
 from pennylane.wires import Wires
 
 from .utils import pauli_eigs
@@ -1411,14 +1410,9 @@ class Operator(abc.ABC):
         """The addition operation of Operator-Operator objects and Operator-scalar."""
         if isinstance(other, Operator):
             return qml.sum(self, other)
-        if other == 0:
-            return self
-        if isinstance(other, qml.pulse.ParametrizedHamiltonian):
-            return other.__add__(self)
-        backend = autoray.infer_backend(other)
-        if (
-            backend == "builtins" and isinstance(other, numbers.Number)
-        ) or backend in SUPPORTED_INTERFACES:
+        if isinstance(other, TensorLike):
+            if qml.math.allequal(other, 0):
+                return self
             return qml.sum(self, qml.s_prod(scalar=other, operator=qml.Identity(self.wires)))
         return NotImplemented
 
@@ -1428,13 +1422,16 @@ class Operator(abc.ABC):
         """The scalar multiplication between scalars and Operators."""
         if callable(other):
             return qml.pulse.ParametrizedHamiltonian([other], [self])
-        backend = autoray.infer_backend(other)
-        if (
-            backend == "builtins" and isinstance(other, numbers.Number)
-        ) or backend in SUPPORTED_INTERFACES:
+        if isinstance(other, TensorLike):
             return qml.s_prod(scalar=other, operator=self)
         return NotImplemented
 
+    def __truediv__(self, other):
+        """The division between an Operator and a number."""
+        if isinstance(other, TensorLike):
+            return self.__mul__(1 / other)
+        return NotImplemented
+
     __rmul__ = __mul__
 
     def __matmul__(self, other):
@@ -1443,8 +1440,8 @@ class Operator(abc.ABC):
 
     def __sub__(self, other):
         """The subtraction operation of Operator-Operator objects and Operator-scalar."""
-        if isinstance(other, (Operator, numbers.Number)):
-            return self + (-other)
+        if isinstance(other, (Operator, TensorLike)):
+            return self + (qml.math.multiply(-1, other))
         return NotImplemented
 
     def __rsub__(self, other):
@@ -1457,10 +1454,7 @@ class Operator(abc.ABC):
 
     def __pow__(self, other):
         r"""The power operation of an Operator object."""
-        backend = autoray.infer_backend(other)
-        if (
-            backend == "builtins" and isinstance(other, numbers.Number)
-        ) or backend in SUPPORTED_INTERFACES:
+        if isinstance(other, TensorLike):
             return qml.pow(self, z=other)
         return NotImplemented
 
diff --git a/pennylane/pulse/__init__.py b/pennylane/pulse/__init__.py
index c63dcce2a..62c7dadf2 100644
--- a/pennylane/pulse/__init__.py
+++ b/pennylane/pulse/__init__.py
@@ -59,6 +59,17 @@ Convenience Functions
     ~pwc_from_function
     ~rect
 
+Hardware Compatible Hamiltonians
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. currentmodule:: pennylane.pulse
+
+.. autosummary::
+    :toctree: api
+
+    ~rydberg_interaction
+    ~rydberg_drive
+
 
 Creating a parametrized Hamiltonian
 -----------------------------------
@@ -245,6 +256,7 @@ JIT-compiling is optional, and one can remove the decorator when only single exe
     See Usage Details of :class:`~.ParametrizedEvolution` for a detailed example.
 """
 
-from .convenience_functions import constant, rect, pwc, pwc_from_function
+from .convenience_functions import constant, pwc, pwc_from_function, rect
 from .parametrized_evolution import ParametrizedEvolution
 from .parametrized_hamiltonian import ParametrizedHamiltonian
+from .rydberg_hamiltonian import RydbergHamiltonian, rydberg_drive, rydberg_interaction
diff --git a/pennylane/pulse/parametrized_hamiltonian.py b/pennylane/pulse/parametrized_hamiltonian.py
index ef72ab4b9..2ffeeacf9 100644
--- a/pennylane/pulse/parametrized_hamiltonian.py
+++ b/pennylane/pulse/parametrized_hamiltonian.py
@@ -18,6 +18,7 @@ This submodule contains the ParametrizedHamiltonian class
 import pennylane as qml
 from pennylane.operation import Operator
 from pennylane.ops.qubit.hamiltonian import Hamiltonian
+from pennylane.typing import TensorLike
 from pennylane.wires import Wires
 
 
@@ -289,20 +290,57 @@ class ParametrizedHamiltonian:
         coeffs = self.coeffs.copy()
 
         if isinstance(H, (Hamiltonian, ParametrizedHamiltonian)):
-            coeffs.extend(H.coeffs.copy())
-            ops.extend(H.ops.copy())
-            return ParametrizedHamiltonian(coeffs, ops)
+            new_coeffs = coeffs + H.coeffs.copy()
+            new_ops = ops + H.ops.copy()
+            return ParametrizedHamiltonian(new_coeffs, new_ops)
 
         if isinstance(H, qml.ops.SProd):  # pylint: disable=no-member
-            coeffs.append(H.scalar)
-            ops.append(H.base)
-            return ParametrizedHamiltonian(coeffs, ops)
+            new_coeffs = coeffs + [H.scalar]
+            new_ops = ops + [H.base]
+            return ParametrizedHamiltonian(new_coeffs, new_ops)
 
         if isinstance(H, Operator):
-            coeffs.append(1)
-            ops.append(H)
-            return ParametrizedHamiltonian(coeffs, ops)
+            new_coeffs = coeffs + [1]
+            new_ops = ops + [H]
 
+            return ParametrizedHamiltonian(new_coeffs, new_ops)
+
+        return NotImplemented
+
+    def __radd__(self, H):
+        r"""The addition operation between a ``ParametrizedHamiltonian`` and an ``Operator``
+        or ``ParametrizedHamiltonian``."""
+        ops = self.ops.copy()
+        coeffs = self.coeffs.copy()
+
+        if isinstance(H, (Hamiltonian, ParametrizedHamiltonian)):
+            new_coeffs = H.coeffs.copy() + coeffs
+            new_ops = H.ops.copy() + ops
+            return ParametrizedHamiltonian(new_coeffs, new_ops)
+
+        if isinstance(H, qml.ops.SProd):  # pylint: disable=no-member
+            new_coeffs = [H.scalar] + coeffs
+            new_ops = [H.base] + ops
+            return ParametrizedHamiltonian(new_coeffs, new_ops)
+
+        if isinstance(H, Operator):
+            new_coeffs = [1] + coeffs
+            new_ops = [H] + ops
+
+            return ParametrizedHamiltonian(new_coeffs, new_ops)
+
+        return NotImplemented
+
+    def __mul__(self, other):
+        ops = self.ops.copy()
+        coeffs_fixed = self.coeffs_fixed.copy()
+        coeffs_parametrized = self.coeffs_parametrized.copy()
+        if isinstance(other, TensorLike) and qml.math.ndim(other) == 0:
+            coeffs_fixed = [other * c for c in coeffs_fixed]
+            coeffs_parametrized = [
+                lambda p, t, new_c=c: other * new_c(p, t) for c in coeffs_parametrized
+            ]
+            return ParametrizedHamiltonian(coeffs_fixed + coeffs_parametrized, ops)
         return NotImplemented
 
-    __radd__ = __add__
+    __rmul__ = __mul__
diff --git a/pennylane/pulse/rydberg_hamiltonian.py b/pennylane/pulse/rydberg_hamiltonian.py
new file mode 100644
index 000000000..144ecf531
--- /dev/null
+++ b/pennylane/pulse/rydberg_hamiltonian.py
@@ -0,0 +1,340 @@
+# Copyright 2018-2023 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""This module contains the classes/functions needed to simulate the evolution of ensembles of
+individual (trapped) rydberg atoms under the excitation of several laser fields."""
+import warnings
+from dataclasses import dataclass
+from typing import Callable, List, Union
+
+import numpy as np
+
+import pennylane as qml
+from pennylane.wires import Wires
+
+from .parametrized_hamiltonian import ParametrizedHamiltonian
+
+
+def rydberg_interaction(register: list, wires=None, interaction_coeff: float = 862690):
+    r"""Returns a :class:`ParametrizedHamiltonian` representing the interaction of an ensemble of
+    Rydberg atoms due to the Rydberg blockade
+
+    .. math::
+
+        \sum_{i<j} V_{ij} n_i n_j
+
+    where :math:`n_i` corresponds to the projector on the Rydberg state of the atom :math:`i`, and
+    :math:`V_{ij}` is the van der Waals potential:
+
+    .. math::
+
+        V_{ij} = \frac{C_6}{R_{ij}^6}
+
+    where :math:`R_{ij}` is the distance between the atoms :math:`i` and :math:`j`, and :math:`C_6`
+    is the Rydberg interaction constant, which defaults to :math:`862690 \text{MHz} \times \mu \text{m}^6`.
+    The unit of time for the evolution of this Rydberg interaction term is in :math:`\mu \text{s}`.
+    This interaction term can be combined with laser drive terms (:func:`~.rydberg_drive`) to create
+    a Hamiltonian describing a driven Rydberg atom system.
+
+    .. seealso::
+
+        :func:`~.rydberg_drive`
+
+    Args:
+        register (list): list of coordinates of the Rydberg atoms (in micrometers)
+        wires (list): List of wires containing the wire values for all the atoms. This list should
+            have the same length as ``register``. If ``None``, each atom's wire value will
+            correspond to its index in the ``register`` list.
+        interaction_coeff (float): Rydberg interaction constant in units: :math:`\text{MHz} \times \mu \text{m}^6`.
+            Defaults to :math:`862690 \text{ MHz} \times \mu \text{m}^6`. This value is based on an assumption that
+            frequencies and energies in the Hamiltonian are provided in units of MHz.
+
+    Returns:
+        RydbergHamiltonian: a :class:`~.ParametrizedHamiltonian` representing the atom interaction
+
+    **Example**
+
+    We create a Hamiltonian describing the van der Waals interaction among 9 Rydberg atoms in a square lattice:
+
+    .. code-block:: python
+
+        atom_coordinates = [[0, 0], [0, 5], [0, 10], [5, 0], [5, 5], [5, 10], [10, 0], [10, 5], [10, 10]]
+        wires = [1, 5, 0, 2, 4, 3, 8, 6, 7]
+        H_i = qml.pulse.rydberg_interaction(atom_coordinates, wires=wires)
+
+    >>> H_i
+    ParametrizedHamiltonian: terms=36
+
+    As expected, we have :math:`\frac{N(N-1)}{2} = 36` terms for N=6 atoms.
+
+    The interaction term is dependent only on the number and positions of the Rydberg atoms. We can execute this
+    pulse program, which corresponds to all driving laser fields being turned off and therefore has no trainable
+    parameters. To add a driving laser field, see :func:`~.rydberg_drive`.
+
+    .. code-block:: python
+
+        dev = qml.device("default.qubit.jax", wires=9)
+
+        @qml.qnode(dev, interface="jax")
+        def circuit():
+            qml.evolve(H_i)([], t=[0, 10])
+            return qml.expval(qml.PauliZ(0))
+
+    >>> circuit()
+    Array(1., dtype=float32)
+    """
+    if wires is None:
+        wires = list(range(len(register)))
+    elif len(wires) != len(register):
+        raise ValueError("The length of the wires and the register must match.")
+
+    coeffs = []
+    observables = []
+    for idx, (pos1, wire1) in enumerate(zip(register[:-1], wires[:-1])):
+        for pos2, wire2 in zip(register[(idx + 1) :], wires[(idx + 1) :]):
+            atom_distance = np.linalg.norm(qml.math.array(pos2) - pos1)
+            Vij = interaction_coeff / (abs(atom_distance) ** 6)  # van der Waals potential
+            coeffs.append(Vij)
+            observables.append(qml.prod(qml.Projector([1], wire1), qml.Projector([1], wire2)))
+            # Rydberg projectors
+
+    return RydbergHamiltonian(
+        coeffs, observables, register=register, interaction_coeff=interaction_coeff
+    )
+
+
+def rydberg_drive(amplitude, detuning, phase, wires):
+    r"""Returns a :class:`ParametrizedHamiltonian` representing the action of a driving laser
+    field with the given rabi frequency, detuning and phase acting on the given wires
+
+    .. math::
+
+        \frac{1}{2} \Omega(t) \sum_{i \in \text{wires}} (\cos(\phi)\sigma_i^x - \sin(\phi)\sigma_i^y) -
+        \frac{1}{2} \delta(t) \sum_{i \in \text{wires}} \sigma_i^z
+
+    where :math:`\Omega`, :math:`\delta` and :math:`\phi` correspond to the rabi frequency, detuning
+    and phase of the laser, :math:`i` correspond to the wire index, and :math:`\sigma^\alpha` for
+    :math:`\alpha = x,y,z` are the Pauli matrices. The unit of time for the  evolution of this Rydberg
+    drive term is :math:`\mu \text{s}`. This driving term can be combined with an interaction term to
+    create a Hamiltonian describing a driven Rydberg atom system. Multiple driving terms can be combined
+    by summing them (see example).
+
+    Args:
+        amplitude (Union[float, Callable]): float or callable returning the amplitude (in MHz) of a
+            laser field
+        detuning (Union[float, Callable]): float or callable returning the detuning (in MHz) of a
+            laser field
+        phase (float): float containing the phase (in radians) of the laser field
+        wires (Union[int, List[int]]): integer or list containing wire values for the Rydberg atoms that
+            the laser field acts on
+
+    Returns:
+        RydbergHamiltonian: a :class:`~.ParametrizedHamiltonian` representing the action of the laser field
+        on the Rydberg atoms.
+
+    .. seealso::
+
+        :func:`~.rydberg_interaction`, :class:`~.ParametrizedHamiltonian`, :class:`~.ParametrizedEvolution`
+        and :func:`~.evolve`
+
+    **Example**
+
+    We create a Hamiltonian describing a laser acting on 4 wires (Rydberg atoms) with a fixed detuning and
+    phase, and a parametrized, time-dependent amplitude. The Hamiltonian includes an interaction term for
+    inter-atom interactions due to van der Waals forces, as well as the driving term for the laser driving
+    the atoms:
+
+    .. code-block:: python
+
+        atom_coordinates = [[0, 0], [0, 4], [4, 0], [4, 4]]
+        H_i = qml.pulse.rydberg_interaction(atom_coordinates, wires)
+
+        amplitude = lambda p, t: p * jnp.sin(jnp.pi * t)
+        detuning = 3 * jnp.pi / 4
+        phase = jnp.pi / 2
+        wires = [0, 1, 2, 3]
+        H_d = qml.pulse.rydberg_drive(amplitude, detuning, phase, wires)
+
+    >>> H_i
+    ParametrizedHamiltonian: terms=6
+    >>> H_d
+    ParametrizedHamiltonian: terms=2
+
+    The two terms of the drive field correspond to the first and second sum, corresponding to the drive and the shift term.
+    This drive term corresponds to a global drive that acts on all wires of the device.
+
+    .. code-block:: python
+
+        dev = qml.device("default.qubit.jax", wires=wires)
+
+        @qml.qnode(dev, interface="jax")
+        def circuit(params):
+            qml.evolve(H_i + H_d)(params, t=[0, 10])
+            return qml.expval(qml.PauliZ(0))
+
+    >>> params = [2.4]
+    >>> circuit(params)
+    Array(0.97137696, dtype=float32)
+    >>> jax.grad(circuit)(params)
+    [Array(0.10493923, dtype=float32)]
+
+    We can also create a Hamiltonian with multiple local drives. The following circuit corresponds to the
+    evolution where an additional local drive acting on wires ``[0, 1]`` is added to the Hamiltonian:
+
+    .. code-block:: python
+
+        amplitude_local = lambda p, t: p[0] * jnp.sin(2 * jnp.pi * t) + p[1]
+        detuning_local = lambda p, t: p * jnp.exp(-0.25 * t)
+        phase_local = jnp.pi / 4
+        H_local = qml.pulse.rydberg_drive(amplitude_local, detuning_local, phase_local, [0, 1])
+
+        H = H_i + H_d + H_local
+
+        @jax.jit
+        @qml.qnode(dev, interface="jax")
+        def circuit_local(params):
+            qml.evolve(H)(params, t=[0, 10])
+            return qml.expval(qml.PauliZ(0))
+
+    >>> params = [2.4, [1.3, -2.0]]
+    >>> circuit_local(params)
+    Array(0.45782223, dtype=float64)
+    >>> jax.grad(circuit_local_drives)(params)
+    [Array(-0.33522988, dtype=float64),
+     [Array(0.40320718, dtype=float64, weak_type=True),
+      Array(-0.12003976, dtype=float64, weak_type=True)]]
+    """
+    if isinstance(wires, int):
+        wires = [wires]
+
+    # We compute the `coeffs` and `observables` of the laser field
+    coeffs = [amplitude, detuning]
+    amplitude_observable = sum(
+        qml.math.cos(phase) * qml.PauliX(wire) - qml.math.sin(phase) * qml.PauliY(wire)
+        for wire in wires
+    )
+    detuning_observable = sum(qml.PauliZ(wire) for wire in wires)
+    observables = [amplitude_observable, detuning_observable]
+
+    # We convert the pulse data into a list of ``RydbergPulse`` objects
+    pulses = [RydbergPulse(amplitude, detuning, phase, wires)]
+    return RydbergHamiltonian(coeffs, observables, pulses=pulses)
+
+
+class RydbergHamiltonian(ParametrizedHamiltonian):
+    r"""Internal class used to keep track of the required information to translate a ``ParametrizedHamiltonian``
+    into hardware.
+
+    This class contains the ``coeffs`` and the ``observables`` that represent one or more
+    terms of the Hamiltonian of an ensemble of Rydberg atoms under the action of local and global
+    laser fields:
+
+    .. math::
+
+        H = \frac{1}{2} \sum_i  \Omega_i(t) (\cos(\phi_i)\sigma_i^x - \sin(\phi_i)\sigma_i^y) -
+        \frac{1}{2} \sum_i \delta_i(t) \sigma_i^z + \sum_{i<j} V_{ij} n_i n_j
+
+    Additionally, it also contains two more attributes (``register`` and ``pulses``) that contain
+    the information that the hardware needs to execute this Hamiltonian.
+
+    .. warning::
+
+        This class should NEVER be initialized directly! Please use the functions
+        :func:`rydberg_interaction` and :func:`rydberg_drive` instead.
+
+    .. seealso:: :func:`rydberg_interaction`, :func:`rydberg_drive`, :class:`ParametrizedHamiltonian`
+
+    Args:
+        coeffs (Union[float, callable]): coefficients of the Hamiltonian expression, which may be
+            constants or parametrized functions. All functions passed as ``coeffs`` must have two
+            arguments, the first one being the trainable parameters and the second one being time.
+        observables (Iterable[Observable]): observables in the Hamiltonian expression, of same
+            length as ``coeffs``
+
+    Keyword Args:
+        register (list): list of coordinates (in micrometers) of each atom in the ensemble
+        pulses (list): list of ``RydbergPulse`` classes containing the information about the
+            amplitude, phase, detuning and wires of each pulse
+        interaction_coeff (float): Rydberg interaction constant in units: :math:`\text{MHz} \times \mu m^6`.
+            Defaults to :math:`862690 \text{MHz} \times \mu m^6`.
+
+    Returns:
+        RydbergHamiltonian: class representing the Hamiltonian of an ensemble of Rydberg atoms
+    """
+
+    # pylint: disable=too-many-arguments
+    def __init__(
+        self,
+        coeffs,
+        observables,
+        register: list = None,
+        pulses: List["RydbergPulse"] = None,
+        interaction_coeff: float = 862690,
+    ):
+        self.register = register
+        self.pulses = [] if pulses is None else pulses
+        self.interaction_coeff = interaction_coeff
+        super().__init__(coeffs, observables)
+
+    def __add__(self, other):
+        if not isinstance(other, RydbergHamiltonian):
+            return super().__add__(other)
+
+        # Update coeffs, obs and hardware attributes
+        if self.register is not None:
+            if other.register is not None:
+                raise ValueError("We cannot add two Hamiltonians with an interaction term!")
+            if not self.wires.contains_wires(other.wires):
+                warnings.warn(
+                    "The wires of the laser fields are not present in the Rydberg ensemble."
+                )
+        elif other.register is not None and not other.wires.contains_wires(self.wires):
+            warnings.warn("The wires of the laser fields are not present in the Rydberg ensemble.")
+
+        new_register = self.register or other.register
+        new_pulses = self.pulses + other.pulses
+        new_ops = self.ops + other.ops
+        new_coeffs = self.coeffs + other.coeffs
+        return RydbergHamiltonian(new_coeffs, new_ops, register=new_register, pulses=new_pulses)
+
+
+@dataclass
+class RydbergPulse:
+    """Dataclass that contains the information of a single Rydberg pulse. This class is used
+    internally in PL to group into a single object all the data related to a single laser field.
+
+    Args:
+        amplitude (Union[float, Callable]): float or callable returning the amplitude (in MHz) of a laser
+            field
+        detuning (Union[float, Callable]): float or callable returning the detuning (in MHz) of a
+            laser field
+        phase (float): float containing the phase (in radians) of the laser field
+        wires (Union[int, List[int]]): integer or list containing wire values that the laser field
+            acts on
+    """
+
+    amplitude: Union[float, Callable]
+    detuning: Union[float, Callable]
+    phase: float
+    wires: List[Wires]
+
+    def __post_init__(self):
+        self.wires = Wires(self.wires)
+
+    def __eq__(self, other):
+        return (
+            self.amplitude == other.amplitude
+            and self.phase == other.phase
+            and self.detuning == other.detuning
+            and self.wires == other.wires
+        )


Test patch:
diff --git a/tests/pulse/test_parametrized_hamiltonian.py b/tests/pulse/test_parametrized_hamiltonian.py
index 000a2e4ed..52a4daa29 100644
--- a/tests/pulse/test_parametrized_hamiltonian.py
+++ b/tests/pulse/test_parametrized_hamiltonian.py
@@ -233,38 +233,48 @@ class TestInteractionWithOperators:
         """Test that a Hamiltonian and SProd can be added to a ParametrizedHamiltonian, and
         will be incorporated in the H_fixed term, with their coefficients included in H_coeffs_fixed
         """
-        pH = ParametrizedHamiltonian([f1, f2], [qml.PauliX(0), qml.PauliY(1)])
+        pH = ParametrizedHamiltonian([f1, f2, 2], [qml.PauliX(0), qml.PauliY(1), qml.PauliZ(2)])
+        pH_fixed = qml.s_prod(2, qml.PauliZ(2))
         params = [1, 2]
         # Adding on the right
         new_pH = pH + H
-        assert pH.H_fixed() == 0
-        assert qml.equal(new_pH.H_fixed(), qml.s_prod(coeff, qml.PauliZ(0)))
-        assert new_pH.coeffs_fixed[0] == coeff
-        assert qml.math.allequal(new_pH(params, t=0.5).matrix(), qml.matrix(pH(params, t=0.5) + H))
+        assert qml.equal(pH.H_fixed(), pH_fixed)
+        assert qml.equal(new_pH.H_fixed(), sum((pH_fixed, qml.s_prod(coeff, qml.PauliZ(0)))))
+        assert new_pH.coeffs_fixed == [2, coeff]
+        assert qml.math.allclose(
+            new_pH(params, t=0.5).matrix(),
+            qml.matrix(pH(params, t=0.5) + H, wire_order=new_pH.wires),
+        )
         # Adding on the left
         new_pH = H + pH
-        assert pH.H_fixed() == 0
-        assert qml.equal(new_pH.H_fixed(), qml.s_prod(coeff, qml.PauliZ(0)))
-        assert new_pH.coeffs_fixed[0] == coeff
-        assert qml.math.allequal(new_pH(params, t=0.5).matrix(), qml.matrix(pH(params, t=0.5) + H))
+        assert qml.equal(pH.H_fixed(), pH_fixed)
+        assert qml.equal(new_pH.H_fixed(), sum((qml.s_prod(coeff, qml.PauliZ(0)), pH_fixed)))
+        assert new_pH.coeffs_fixed == [coeff, 2]
+        assert qml.math.allclose(
+            new_pH(params, t=0.5).matrix(),
+            qml.matrix(pH(params, t=0.5) + H, wire_order=new_pH.wires),
+        )
 
     @pytest.mark.parametrize("op", ops)
     def test_add_other_operators(self, op):
         """Test that a Hamiltonian, SProd, Tensor or Operator can be added to a
         ParametrizedHamiltonian, and will be incorporated in the H_fixed term"""
-        pH = ParametrizedHamiltonian([f1, f2], [qml.PauliX(0), qml.PauliY(1)])
+        pH = ParametrizedHamiltonian([f1, f2, 2], [qml.PauliX(0), qml.PauliY(1), qml.PauliZ(2)])
+        pH_fixed = qml.s_prod(2, qml.PauliZ(2))
 
         # Adding on the right
         new_pH = pH + op
-        assert pH.H_fixed() == 0
-        assert qml.equal(new_pH.H_fixed(), qml.s_prod(1, op))
+        assert qml.equal(pH.H_fixed(), pH_fixed)
+        assert qml.equal(new_pH.H_fixed(), sum((pH_fixed, qml.s_prod(1, op))))
 
         # Adding on the left
         new_pH = op + pH
-        assert pH.H_fixed() == 0
-        assert qml.equal(new_pH.H_fixed(), qml.s_prod(1, op))
+        assert qml.equal(pH.H_fixed(), pH_fixed)
+        assert qml.equal(new_pH.H_fixed(), sum((qml.s_prod(1, op), pH_fixed)))
 
     def test_add_invalid_object_raises_error(self):
+        """Test that an error is raised when adding a ``ParametrizedHamiltonian`` with an invalid
+        object."""
         H = ParametrizedHamiltonian([f1, f2], [qml.PauliX(0), qml.PauliY(1)])
 
         class DummyObject:  # pylint: disable=too-few-public-methods
@@ -273,6 +283,33 @@ class TestInteractionWithOperators:
         with pytest.raises(TypeError, match="unsupported operand type"):
             _ = H + DummyObject()
 
+        with pytest.raises(TypeError, match="unsupported operand type"):
+            _ = DummyObject() + H
+
+    def test_multiply_with_scalar(self):
+        """Test the __mul__ dunder method with a scalar."""
+        H = ParametrizedHamiltonian([f1, f2], [qml.PauliX(0), qml.PauliY(1)])
+        new_H = 3 * H
+        expected_H = ParametrizedHamiltonian(
+            [lambda p, t: 3 * f1(p, t), lambda p, t: 3 * f2(p, t)], [qml.PauliX(0), qml.PauliY(1)]
+        )
+        assert new_H.H_fixed() == expected_H.H_fixed() == 0
+        assert qml.equal(new_H.H_parametrized([1, 2], t=4), expected_H.H_parametrized([1, 2], t=4))
+
+    def test_multiply_invalid_object_raises_error(self):
+        """Test that an error is raised when multiplying a ``ParametrizedHamiltonian`` with an invalid
+        object."""
+        H = ParametrizedHamiltonian([f1, f2], [qml.PauliX(0), qml.PauliY(1)])
+
+        class DummyObject:  # pylint: disable=too-few-public-methods
+            pass
+
+        with pytest.raises(TypeError, match="unsupported operand type"):
+            _ = H * DummyObject()
+
+        with pytest.raises(TypeError, match="unsupported operand type"):
+            _ = DummyObject() * H
+
     def test_adding_two_parametrized_hamiltonians(self):
         """Test that two ParametrizedHamiltonians can be added together and
         the H_fixed and H_parametrized terms are correctly combined."""
diff --git a/tests/pulse/test_rydberg_hamiltonian.py b/tests/pulse/test_rydberg_hamiltonian.py
new file mode 100644
index 000000000..ef3d8fc4f
--- /dev/null
+++ b/tests/pulse/test_rydberg_hamiltonian.py
@@ -0,0 +1,385 @@
+# Copyright 2018-2023 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+Unit tests for the RydbergHamiltonian class.
+"""
+import warnings
+
+# pylint: disable=too-few-public-methods
+import numpy as np
+import pytest
+
+import pennylane as qml
+from pennylane.pulse import RydbergHamiltonian, rydberg_drive, rydberg_interaction
+from pennylane.pulse.rydberg_hamiltonian import RydbergPulse
+from pennylane.wires import Wires
+
+atom_coordinates = [[0, 0], [0, 5], [5, 0], [10, 5], [5, 10], [10, 10]]
+wires = [1, 6, 0, 2, 4, 3]
+
+
+class TestRydbergHamiltonian:
+    """Unit tests for the properties of the RydbergHamiltonian class."""
+
+    # pylint: disable=protected-access
+    def test_initialization(self):
+        """Test the RydbergHamiltonian class is initialized correctly."""
+        rm = RydbergHamiltonian(coeffs=[], observables=[], register=atom_coordinates)
+
+        assert qml.math.allequal(rm.register, atom_coordinates)
+        assert rm.pulses == []
+        assert rm.wires == Wires([])
+        assert rm.interaction_coeff == 862690
+
+    def test_add(self):
+        """Test that the __add__ dunder method works correctly."""
+        rm1 = RydbergHamiltonian(
+            coeffs=[1, 2],
+            observables=[qml.PauliX(4), qml.PauliZ(8)],
+            register=atom_coordinates,
+            pulses=[RydbergPulse(1, 2, 3, [4, 8])],
+        )
+        rm2 = RydbergHamiltonian(
+            coeffs=[2],
+            observables=[qml.PauliY(8)],
+            pulses=[RydbergPulse(5, 6, 7, 8)],
+        )
+        with warnings.catch_warnings():
+            # We make sure that no warning is raised
+            warnings.simplefilter("error")
+            sum_rm = rm1 + rm2
+        assert isinstance(sum_rm, RydbergHamiltonian)
+        assert qml.math.allequal(sum_rm.coeffs, [1, 2, 2])
+        assert all(
+            qml.equal(op1, op2)
+            for op1, op2 in zip(sum_rm.ops, [qml.PauliX(4), qml.PauliZ(8), qml.PauliY(8)])
+        )
+        assert qml.math.allequal(sum_rm.register, atom_coordinates)
+        assert sum_rm.pulses == [RydbergPulse(1, 2, 3, [4, 8]), RydbergPulse(5, 6, 7, 8)]
+
+    def test_add_parametrized_hamiltonian(self):
+        """Tests that adding a `RydbergHamiltonian` and `ParametrizedHamiltonian` works as
+        expected."""
+        coeffs = [2, 3]
+        ops = [qml.PauliZ(0), qml.PauliX(2)]
+        h_wires = [0, 2]
+
+        rh = RydbergHamiltonian(
+            coeffs=[coeffs[0]],
+            observables=[ops[0]],
+            pulses=[RydbergPulse(5, 6, 7, 8)],
+        )
+        ph = qml.pulse.ParametrizedHamiltonian(coeffs=[coeffs[1]], observables=[ops[1]])
+
+        res1 = rh + ph
+        res2 = ph + rh
+
+        assert res1.coeffs_fixed == coeffs
+        assert res1.coeffs_parametrized == []
+        assert all(qml.equal(op1, op2) for op1, op2 in zip(res1.ops_fixed, ops))
+        assert res1.ops_parametrized == []
+        assert res1.wires == qml.wires.Wires(h_wires)
+
+        coeffs.reverse()
+        ops.reverse()
+        h_wires.reverse()
+
+        assert res2.coeffs_fixed == coeffs
+        assert res2.coeffs_parametrized == []
+        assert all(qml.equal(op1, op2) for op1, op2 in zip(res2.ops_fixed, ops))
+        assert res2.ops_parametrized == []
+        assert res2.wires == qml.wires.Wires(h_wires)
+
+    def test_add_raises_error(self):
+        """Test that an error is raised if two RydbergHamiltonians with registers are added."""
+        rm1 = RydbergHamiltonian(
+            coeffs=[1],
+            observables=[qml.PauliX(0)],
+            register=atom_coordinates,
+            pulses=[RydbergPulse(1, 2, 3, 4)],
+        )
+        with pytest.raises(
+            ValueError, match="We cannot add two Hamiltonians with an interaction term"
+        ):
+            _ = rm1 + rm1
+
+    def test_add_raises_warning(self):
+        """Test that an error is raised when adding two RydbergHamiltonians where one Hamiltonian
+        contains pulses on wires that are not present in the register."""
+        coords = [[0, 0], [0, 5], [5, 0]]
+
+        Hd = rydberg_interaction(register=coords, wires=[0, 1, 2])
+        Ht = rydberg_drive(2, 3, 4, wires=3)
+
+        with pytest.warns(
+            UserWarning,
+            match="The wires of the laser fields are not present in the Rydberg ensemble",
+        ):
+            _ = Hd + Ht
+
+        with pytest.warns(
+            UserWarning,
+            match="The wires of the laser fields are not present in the Rydberg ensemble",
+        ):
+            _ = Ht + Hd
+
+
+def f1(p, t):
+    return p * np.sin(t) * (t - 1)
+
+
+def f2(p, t):
+    return p * np.cos(t**2)
+
+
+param = [1.2, 2.3]
+
+
+class TestInteractionWithOperators:
+    """Test that the interaction between a ``RydbergHamiltonian`` and other operators work as
+    expected."""
+
+    ops_with_coeffs = (
+        (qml.Hamiltonian([2], [qml.PauliZ(0)]), 2),
+        (qml.Hamiltonian([1.7], [qml.PauliZ(0)]), 1.7),
+        (qml.ops.SProd(3, qml.PauliZ(0)), 3),
+    )
+    ops = (
+        qml.PauliX(2),
+        qml.PauliX(2) @ qml.PauliX(3),
+        qml.CNOT([0, 1]),
+    )  # ToDo: maybe add more operators to test here?
+
+    @pytest.mark.parametrize("H, coeff", ops_with_coeffs)
+    def test_add_special_operators(self, H, coeff):
+        """Test that a Hamiltonian and SProd can be added to a RydbergHamiltonian, and
+        will be incorporated in the H_fixed term, with their coefficients included in H_coeffs_fixed.
+        """
+        R = rydberg_drive(amplitude=f1, detuning=f2, phase=0, wires=[0, 1])
+        params = [1, 2]
+        # Adding on the right
+        new_pH = R + H
+        assert R.H_fixed() == 0
+        assert qml.equal(new_pH.H_fixed(), qml.s_prod(coeff, qml.PauliZ(0)))
+        assert new_pH.coeffs_fixed[0] == coeff
+        assert qml.math.allequal(new_pH(params, t=0.5).matrix(), qml.matrix(R(params, t=0.5) + H))
+        # Adding on the left
+        new_pH = H + R
+        assert R.H_fixed() == 0
+        assert qml.equal(new_pH.H_fixed(), qml.s_prod(coeff, qml.PauliZ(0)))
+        assert new_pH.coeffs_fixed[0] == coeff
+        assert qml.math.allequal(new_pH(params, t=0.5).matrix(), qml.matrix(R(params, t=0.5) + H))
+
+    @pytest.mark.parametrize("op", ops)
+    def test_add_other_operators(self, op):
+        """Test that a Hamiltonian, SProd, Tensor or Operator can be added to a
+        ParametrizedHamiltonian, and will be incorporated in the H_fixed term"""
+        R = rydberg_drive(amplitude=f1, detuning=f2, phase=0, wires=[0, 1])
+
+        # Adding on the right
+        new_pH = R + op
+        assert R.H_fixed() == 0
+        assert qml.equal(new_pH.H_fixed(), qml.s_prod(1, op))
+
+        # Adding on the left
+        new_pH = op + R
+        assert R.H_fixed() == 0
+        assert qml.equal(new_pH.H_fixed(), qml.s_prod(1, op))
+
+
+class TestRydbergInteraction:
+    """Unit tests for the ``rydberg_interaction`` function."""
+
+    def test_attributes_and_number_of_terms(self):
+        """Test that the attributes and the number of terms of the ``ParametrizedHamiltonian`` returned by
+        ``rydberg_interaction`` are correct."""
+        Hd = rydberg_interaction(register=atom_coordinates, wires=wires, interaction_coeff=1)
+
+        assert isinstance(Hd, RydbergHamiltonian)
+        assert Hd.interaction_coeff == 1
+        assert Hd.wires == Wires(wires)
+        assert qml.math.allequal(Hd.register, atom_coordinates)
+        N = len(wires)
+        num_combinations = N * (N - 1) / 2  # number of terms on the rydberg_interaction hamiltonian
+        assert len(Hd.ops) == num_combinations
+        assert Hd.pulses == []
+
+    def test_wires_is_none(self):
+        """Test that when wires is None the wires correspond to an increasing list of values with
+        the same length as the atom coordinates."""
+        Hd = rydberg_interaction(register=atom_coordinates)
+
+        assert Hd.wires == Wires(list(range(len(atom_coordinates))))
+
+    def test_coeffs(self):
+        """Test that the generated coefficients are correct."""
+        coords = [[0, 0], [0, 1], [1, 0]]
+        Hd = rydberg_interaction(coords, interaction_coeff=1)
+        assert Hd.coeffs == [1, 1, 1 / np.sqrt(2) ** 6]
+
+    def test_different_lengths_raises_error(self):
+        """Test that using different lengths for the wires and the register raises an error."""
+        with pytest.raises(ValueError, match="The length of the wires and the register must match"):
+            _ = rydberg_interaction(register=atom_coordinates, wires=[0])
+
+
+class TestRydbergDrive:
+    """Unit tests for the ``rydberg_drive`` function."""
+
+    def test_attributes_and_number_of_terms(self):
+        """Test that the attributes and the number of terms of the ``ParametrizedHamiltonian`` returned by
+        ``rydberg_drive`` are correct."""
+        Hd = rydberg_drive(amplitude=1, detuning=2, phase=3, wires=[1, 2])
+
+        assert isinstance(Hd, RydbergHamiltonian)
+        assert Hd.interaction_coeff == 862690
+        assert Hd.wires == Wires([1, 2])
+        assert Hd.register is None
+        assert len(Hd.ops) == 2  # amplitude and detuning terms of the Hamiltonian
+        assert Hd.pulses == [RydbergPulse(1, 2, 3, [1, 2])]
+
+    def test_multiple_local_drives(self):
+        """Test that adding multiple drive terms behaves as expected"""
+
+        def fa(p, t):
+            return np.sin(p * t)
+
+        def fb(p, t):
+            return np.cos(p * t)
+
+        H1 = rydberg_drive(amplitude=fa, detuning=3, phase=1, wires=[0, 3])
+        H2 = rydberg_drive(amplitude=1, detuning=fb, phase=3, wires=[1, 2])
+        Hd = H1 + H2
+
+        ops_expected = [
+            qml.Hamiltonian([1, 1], [qml.PauliZ(0), qml.PauliZ(3)]),
+            qml.Hamiltonian(
+                [np.cos(3), -np.sin(3), np.cos(3), -np.sin(3)],
+                [qml.PauliX(1), qml.PauliY(1), qml.PauliX(2), qml.PauliY(2)],
+            ),
+            qml.Hamiltonian(
+                [np.cos(1), -np.sin(1), np.cos(1), -np.sin(1)],
+                [qml.PauliX(0), qml.PauliY(0), qml.PauliX(3), qml.PauliY(3)],
+            ),
+            qml.Hamiltonian([1, 1], [qml.PauliZ(1), qml.PauliZ(2)]),
+        ]
+        coeffs_expected = [3, 1, fa, fb]
+        H_expected = qml.dot(coeffs_expected, ops_expected)
+
+        assert isinstance(Hd, RydbergHamiltonian)
+        assert Hd.interaction_coeff == 862690
+        assert Hd.wires == Wires([0, 3, 1, 2])
+        assert Hd.register is None
+        assert len(Hd.ops) == 4  # amplitude and detuning terms of the Hamiltonian
+        # Callable coefficients are shifted to the end of the list.
+        assert Hd.coeffs[0:2] == [3, 1]
+        assert Hd.coeffs[2] is fa
+        assert Hd.coeffs[3] is fb
+        assert len(Hd.pulses) == 2
+        assert Hd.pulses == H1.pulses + H2.pulses
+        assert qml.equal(Hd([0.5, -0.5], t=5), H_expected([0.5, -0.5], t=5))
+
+
+class TestRydbergPulse:
+    """Unit tests for the ``RydbergPulse`` class."""
+
+    def test_init(self):
+        """Test the initialization of the ``RydbergPulse`` class."""
+        p = RydbergPulse(amplitude=4, phase=8, detuning=9, wires=[0, 4, 7])
+        assert p.amplitude == 4
+        assert p.phase == 8
+        assert p.detuning == 9
+        assert p.wires == Wires([0, 4, 7])
+
+    def test_equal(self):
+        """Test the ``__eq__`` method of the ``RydbergPulse`` class."""
+        p1 = RydbergPulse(1, 2, 3, [0, 1])
+        p2 = RydbergPulse(1, 2, 3, 0)
+        p3 = RydbergPulse(1, 2, 3, [0, 1])
+        assert p1 != p2
+        assert p2 != p3
+        assert p1 == p3
+
+
+class TestIntegration:
+    """Integration tests for the ``RydbergHamiltonian`` class."""
+
+    @pytest.mark.jax
+    def test_jitted_qnode(self):
+        """Test that a ``RydbergHamiltonian`` class can be executed within a jitted qnode."""
+        import jax
+        import jax.numpy as jnp
+
+        Hd = rydberg_interaction(register=atom_coordinates, wires=wires)
+
+        def fa(p, t):
+            return jnp.polyval(p, t)
+
+        def fb(p, t):
+            return p[0] * jnp.sin(p[1] * t)
+
+        Ht = rydberg_drive(amplitude=fa, detuning=fb, phase=0, wires=1)
+
+        dev = qml.device("default.qubit", wires=wires)
+
+        ts = jnp.array([0.0, 3.0])
+        H_obj = sum(qml.PauliZ(i) for i in range(2))
+
+        @jax.jit
+        @qml.qnode(dev, interface="jax")
+        def qnode(params):
+            qml.evolve(Hd + Ht)(params, ts)
+            return qml.expval(H_obj)
+
+        params = (jnp.ones(5), jnp.array([1.0, jnp.pi]))
+        res = qnode(params)
+
+        assert isinstance(res, jax.Array)
+
+    @pytest.mark.jax
+    def test_jitted_qnode_multidrive(self):
+        """Test that a ``RydbergHamiltonian`` class with multiple drive terms can be
+        executed within a jitted qnode."""
+        import jax
+        import jax.numpy as jnp
+
+        Hd = rydberg_interaction(register=atom_coordinates, wires=wires)
+
+        def fa(p, t):
+            return jnp.polyval(p, t)
+
+        def fb(p, t):
+            return p[0] * jnp.sin(p[1] * t)
+
+        def fc(p, t):
+            return p[0] * jnp.sin(t) + jnp.cos(p[1] * t)
+
+        H1 = rydberg_drive(amplitude=fa, detuning=fb, phase=0, wires=1)
+        H2 = rydberg_drive(amplitude=fc, detuning=jnp.pi / 4, phase=3 * jnp.pi, wires=4)
+
+        dev = qml.device("default.qubit", wires=wires)
+
+        ts = jnp.array([0.0, 3.0])
+        H_obj = sum(qml.PauliZ(i) for i in range(2))
+
+        @jax.jit
+        @qml.qnode(dev, interface="jax")
+        def qnode(params):
+            qml.evolve(Hd + H1 + H2)(params, ts)
+            return qml.expval(H_obj)
+
+        params = (jnp.ones(5), jnp.array([1.0, jnp.pi]), jnp.array([jnp.pi / 2, 0.5]))
+        res = qnode(params)
+
+        assert isinstance(res, jax.Array)
diff --git a/tests/test_operation.py b/tests/test_operation.py
index 2b2613079..13d986d0c 100644
--- a/tests/test_operation.py
+++ b/tests/test_operation.py
@@ -1033,6 +1033,33 @@ class TestOperatorIntegration:
         assert isinstance(prod_op, SProd)
         assert prod_op.scalar is scalar
 
+    def test_divide_with_scalar(self):
+        """Test the __truediv__ dunder method with a scalar value."""
+        sprod_op = qml.RX(1, 0) / 4
+        final_op = qml.s_prod(scalar=1 / 4, operator=qml.RX(1, 0))
+        assert isinstance(sprod_op, qml.ops.SProd)
+        assert sprod_op.name == final_op.name
+        assert sprod_op.wires == final_op.wires
+        assert sprod_op.data == final_op.data
+        assert np.allclose(sprod_op.matrix(), final_op.matrix(), rtol=0)
+
+    def test_divide_with_scalar_tensor(self):
+        """Test the __truediv__ dunder method with a scalar tensor."""
+        scalar = pnp.array(5)
+        prod_op = qml.RX(1.23, 0) / scalar
+        assert isinstance(prod_op, SProd)
+        assert prod_op.scalar == 1 / scalar
+
+    def test_divide_not_supported(self):
+        """Test that the division of an operator with an unknown object is not supported."""
+        obs = qml.PauliX(0)
+
+        class UnknownObject:
+            pass
+
+        with pytest.raises(TypeError, match="unsupported operand type"):
+            _ = obs / UnknownObject()
+
     def test_dunder_method_with_new_class(self):
         """Test that when calling any Operator dunder method with a non-supported class that
         has its right dunder method defined, the class' right dunder method is called."""
@@ -1895,24 +1922,6 @@ class TestTensorObservableOperations:
         """Tests the tensor product between Observables"""
         assert res.compare(obs1 @ obs2)
 
-    def test_arithmetic_errors(self):
-        """Tests that the arithmetic operations throw the correct errors"""
-        obs = qml.PauliZ(0)
-        tensor = qml.PauliZ(0) @ qml.PauliX(1)
-        A = [[1, 0], [0, -1]]
-        with pytest.raises(TypeError, match="unsupported operand type"):
-            _ = obs + A
-        with pytest.raises(TypeError, match="unsupported operand type"):
-            _ = tensor + A
-        with pytest.raises(TypeError, match="can't multiply sequence by non-int of type 'PauliZ'"):
-            _ = obs * A
-        with pytest.raises(TypeError, match="can't multiply sequence by non-int of type 'Tensor'"):
-            _ = A * tensor
-        with pytest.raises(TypeError, match="unsupported operand type"):
-            _ = obs - A
-        with pytest.raises(TypeError, match="unsupported operand type"):
-            _ = tensor - A
-
 
 # Dummy class inheriting from Operator
 class MyOp(Operator):