diff --git a/docs/source/index.rst b/docs/source/index.rst
index 5d9b5dcb..63eadfc9 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -31,6 +31,7 @@ See the `contributing page <https://github.com/ICHEC/qse/blob/main/CONTRIBUTIONS
 
    tutorials/creating_and_manipulating_qbits
    tutorials/generating_lattices
+   tutorials/working_with_operators
    tutorials/pulser-calc-example
    tutorials/ssh_model
 
diff --git a/docs/source/tutorials/working_with_operators.ipynb b/docs/source/tutorials/working_with_operators.ipynb
new file mode 100644
index 00000000..bcab74d0
--- /dev/null
+++ b/docs/source/tutorials/working_with_operators.ipynb
@@ -0,0 +1,274 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Working with operators\n",
+    "`qse` provides an `Operator` class and an `Operators` class to easily manipulate quantum operators. The `Operators` class can be easily created from a `Qbits` object, making creating interaction hamiltonians very convienient."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import qse"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The Operator class\n",
+    "The `Operator` class represents a single quantum operator."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# To create 0.3 ZII\n",
+    "qse.Operator(operator=\"Z\", qubits=0, nqbits=3, coef=0.3)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# To create -2. IIXIY\n",
+    "qse.Operator(operator=[\"X\", \"Y\"], qubits=[2, 4], nqbits=5, coef=-2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Operators can be multiplied by a float, which multiplies the coefficient."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "op = qse.Operator(\"Z\", 0, 3)\n",
+    "print(op)\n",
+    "print(op * 4.4)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "An operator can be converted into a QuTiP object by the `to_qutip` method."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create ZII\n",
+    "# Note that `coef` defaults to 1\n",
+    "op = qse.Operator(operator=\"Z\", qubits=0, nqbits=3)\n",
+    "op.to_qutip()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# The Operators class\n",
+    "For dealing with multiple `Operator` classes we have the `Operators` class.\n",
+    "\n",
+    "We can fill the `Operators` class either by passing `Operator` classes as a list or by addition."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Creating Operators by passing a list\n",
+    "nqbits = 3\n",
+    "qse.Operators(\n",
+    "    [qse.Operator(operator=\"Z\", qubits=i, nqbits=nqbits) for i in range(nqbits)]\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Creating Operators by addition\n",
+    "hamiltonian = qse.Operators()\n",
+    "\n",
+    "for i in range(nqbits):\n",
+    "    hamiltonian += qse.Operator(operator=\"Z\", qubits=i, nqbits=nqbits)\n",
+    "\n",
+    "hamiltonian"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can use both methods to create complex models, for example to create a Heisenberg Model\n",
+    "$$ H = -J \\sum_{i=1}^{N-1} (X_i X_{i+1} + Y_i Y_{i+1} + Z_i Z_{i+1})$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# For example\n",
+    "nqbits = 3\n",
+    "coupling = -1.0\n",
+    "ham_xx = qse.Operators(\n",
+    "    [qse.Operator(\"X\", [i, i + 1], nqbits, coupling) for i in range(nqbits - 1)]\n",
+    ")\n",
+    "ham_yy = qse.Operators(\n",
+    "    [qse.Operator(\"Y\", [i, i + 1], nqbits, coupling) for i in range(nqbits - 1)]\n",
+    ")\n",
+    "ham_zz = qse.Operators(\n",
+    "    [qse.Operator(\"Z\", [i, i + 1], nqbits, coupling) for i in range(nqbits - 1)]\n",
+    ")\n",
+    "ham = ham_xx + ham_yy + ham_zz\n",
+    "ham"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The `Operators` class also supports conversion to `QuTiP`"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ham.to_qutip()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Creating an interaction Hamiltonian from a Qbits object\n",
+    "`Qbits` objects have a method `compute_interaction_hamiltonian` which allows one to easily create an interaction Hamiltonian from a `Qbits` geometry."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create a nearest-neighbor Ising model on a hexagonal lattice\n",
+    "qbits = qse.lattices.hexagonal(1.0, 2, 3)\n",
+    "del qbits[[0, 11]]\n",
+    "qbits.labels = range(qbits.nqbits)  # relabel\n",
+    "qbits.draw(show_labels=True, radius=\"nearest\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "coupling = 1.0\n",
+    "\n",
+    "\n",
+    "# You can either explicity create a callable to pass to the function\n",
+    "def distance_func(d):\n",
+    "    return coupling * (d < 1.01)  # will be non-zero for nearest-neighbours only\n",
+    "\n",
+    "\n",
+    "ham_int = qbits.compute_interaction_hamiltonian(distance_func, interaction=\"Z\")\n",
+    "\n",
+    "# Or use a lambda function\n",
+    "# ham_int = qbits.compute_interaction_hamiltonian(lambda d: coupling * (d < 1.01), interaction=\"Z\")\n",
+    "\n",
+    "ham_int"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Now add an external field to create the full hamiltonian\n",
+    "coupling_ext = 0.1\n",
+    "ham_ext = qse.Operators(\n",
+    "    [qse.Operator(\"X\", i, qbits.nqbits, coupling_ext) for i in range(qbits.nqbits)]\n",
+    ")\n",
+    "ham_ext"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ham = ham_int + ham_ext\n",
+    "ham"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Version"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "qse.utils.print_environment()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "qse_testbed",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/qse/operator/operator.py b/qse/operator/operator.py
index 509fafaf..756af3b9 100644
--- a/qse/operator/operator.py
+++ b/qse/operator/operator.py
@@ -47,6 +47,11 @@ def __init__(self, operator, qubits, nqbits, coef=1.0):
             raise Exception(
                 "The number of passed qubits must equal the number of passed operators."
             )
+        if max(qubits) >= nqbits:
+            raise Exception(
+                "One or more of the qubits indicies is out of range. "
+                "They must each be less than nqbits."
+            )
 
         self.operator = operator
         self.qubits = qubits
@@ -90,6 +95,21 @@ def __repr__(self):
             [f"{op}{q}" for op, q in zip(self.operator, self.qubits)]
         )
 
+    def copy(self):
+        """
+        Return a copy.
+        """
+        return self.__class__(self.operator, self.qubits, self.nqbits, self.coef)
+
+    def __mul__(self, other):
+        op = self.copy()
+        op *= other
+        return op
+
+    def __imul__(self, other):
+        self.coef *= float(other)
+        return self
+
 
 def _check_operator(op_list):
     for op in op_list:
diff --git a/qse/operator/operators.py b/qse/operator/operators.py
index 0c5a5184..8bf9a5ab 100644
--- a/qse/operator/operators.py
+++ b/qse/operator/operators.py
@@ -69,7 +69,8 @@ def extend(self, other):
         other: Operators or Operator
             The operators to be added to the current object.
         """
-        if other.nqbits != self.nqbits:
+        # We only run this check when nterms > 0
+        if (other.nqbits != self.nqbits) and self.nterms:
             raise Exception("other must have the same number of qubits.")
 
         if isinstance(other, Operator):
@@ -79,6 +80,12 @@ def extend(self, other):
         else:
             raise Exception("other must be Operators or Operator.")
 
+    def copy(self):
+        """
+        Return a copy.
+        """
+        return self.__class__(self.operator_list)
+
     def __add__(self, other):
         op = self.copy()
         op += other
@@ -100,6 +107,8 @@ def __getitem__(self, i):
 
     @property
     def nqbits(self):
+        if self.nterms == 0:
+            return 0
         return self[0].nqbits
 
     @property