diff --git a/docs/_toc.yml b/docs/_toc.yml
index ff2049b..92825e6 100644
--- a/docs/_toc.yml
+++ b/docs/_toc.yml
@@ -11,6 +11,7 @@ parts:
       - file: tutorials/generating_lattices
       - file: tutorials/pulser-calc-example
       - file: tutorials/finite-geometries
+      - file: tutorials/operators
       - file: tutorials/vqe_qiskit
   - caption: Use-Cases
     chapters:
diff --git a/docs/tutorials/operators.ipynb b/docs/tutorials/operators.ipynb
new file mode 100644
index 0000000..658b927
--- /dev/null
+++ b/docs/tutorials/operators.ipynb
@@ -0,0 +1,193 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "f2f2f320",
+   "metadata": {},
+   "source": [
+    "# The Operator and Operators classes\n",
+    "\n",
+    "The Operator and Operators classes provide a way to interact with Operators and Hamiltonians."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3b08d35b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import qse"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a042dded",
+   "metadata": {},
+   "source": [
+    "## The Operator class\n",
+    "\n",
+    "The Operator class represents a single operator acting on a multiple qubit space.\n",
+    "\n",
+    "For example:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "82b62140",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create an operator acting with a Pauli X on the\n",
+    "# third qubit of a 4-qubit system with coefficient 1.0\n",
+    "qse.Operator(\"X\", qubits=2, nqbits=4, coef=1.0)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cc3a8639",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create an operator acting with a Pauli Z on the\n",
+    "# first qubit and a Pauli Y on the second qubit of a 3-qubit system with coefficient 4.0\n",
+    "qse.Operator([\"Z\", \"Y\"], qubits=[0, 1], nqbits=3, coef=4.0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5ae6da3e",
+   "metadata": {},
+   "source": [
+    "The Operator class supports multiplication, where the coefficient will be multiplied by a scalar."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "90458598",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "op = qse.Operator(\"X\", qubits=0, nqbits=2)\n",
+    "print(op)\n",
+    "\n",
+    "op *= 3.14\n",
+    "print(op)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "566f7566",
+   "metadata": {},
+   "source": [
+    "## The Operators class\n",
+    "\n",
+    "The Operators class is for representing a sum of Operator classes.\n",
+    "\n",
+    "For Example:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "960dd86d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "op1 = qse.Operator(\"X\", qubits=0, nqbits=2)\n",
+    "op2 = qse.Operator(\"Y\", qubits=1, nqbits=2)\n",
+    "ops = qse.Operators([op1, op2])\n",
+    "ops"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "96fe72f8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ops = qse.Operators([qse.Operator(\"X\", [i, i + 1], nqbits=4) for i in range(3)])\n",
+    "ops"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cad1861a",
+   "metadata": {},
+   "source": [
+    "## Creating Operators from Qbits\n",
+    "Operators can be created from the Qbits class, given an interaction function."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7442e079",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "spacing = 1.0\n",
+    "qbits = qse.lattices.ring(spacing, 12)\n",
+    "\n",
+    "\n",
+    "def distance_func(d):\n",
+    "    # Nearest neighbours interaction\n",
+    "    return -1.0 * (d < spacing * 1.01)\n",
+    "\n",
+    "\n",
+    "ham_int = qbits.compute_interaction_hamiltonian(\n",
+    "    distance_func=distance_func, interaction=\"Z\"\n",
+    ")\n",
+    "ham_int"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "10fb5e24",
+   "metadata": {},
+   "source": [
+    "We can then combine this with other terms to make complex Hamiltonians. For example, the transverse Ising model:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cb9ac4d4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ham_ext = qse.Operators(\n",
+    "    [qse.Operator(\"X\", i, nqbits=qbits.nqbits) for i in range(qbits.nqbits)]\n",
+    ")\n",
+    "\n",
+    "ham = ham_int + ham_ext\n",
+    "ham"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "qse",
+   "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.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/qse/operator/operator.py b/qse/operator/operator.py
index b68df54..f384230 100644
--- a/qse/operator/operator.py
+++ b/qse/operator/operator.py
@@ -43,6 +43,12 @@ def __init__(self, operator, qubits, nqbits, coef=1.0):
 
         _check_operator(operator)
 
+        for q in qubits:
+            if q < 0 or q >= nqbits:
+                raise Exception(
+                    "Qubit indices must be between 0 and nqbits-1, inclusive."
+                )
+
         if len(qubits) != len(operator):
             raise Exception(
                 "The number of passed qubits must equal the number of passed operators."
@@ -85,11 +91,36 @@ def to_qutip(self):
             op[qi] = _qutip_converter(op_str)
         return self.coef * qp.tensor(op)
 
+    def copy(self):
+        """
+        Creates a copy of the operator.
+
+        Returns
+        -------
+        Operator
+            A new instance of Operator with the same attributes.
+        """
+        return Operator(self.operator, self.qubits, self.nqbits, self.coef)
+
     def __repr__(self):
         return f"{self.coef:.2f} " + " ".join(
             [f"{op}{q}" for op, q in zip(self.operator, self.qubits)]
         )
 
+    def __mul__(self, other):
+        if isinstance(other, (int, float)):
+            return Operator(self.operator, self.qubits, self.nqbits, self.coef * other)
+        else:
+            raise NotImplementedError("Multiplication only supported for scalars.")
+
+    def __imul__(self, other):
+        if isinstance(other, (int, float)):
+            self.coef *= other
+            return self
+
+    def __rmul__(self, other):
+        return self.__mul__(other)
+
 
 def _check_operator(op_list):
     for op in op_list:
diff --git a/qse/operator/operators.py b/qse/operator/operators.py
index 0c5a518..ce1129a 100644
--- a/qse/operator/operators.py
+++ b/qse/operator/operators.py
@@ -79,6 +79,17 @@ def extend(self, other):
         else:
             raise Exception("other must be Operators or Operator.")
 
+    def copy(self):
+        """
+        Creates a copy of the Operators object.
+
+        Returns
+        -------
+        Operators
+            A copy of the current Operators object.
+        """
+        return Operators(self.operator_list.copy())
+
     def __add__(self, other):
         op = self.copy()
         op += other
@@ -98,6 +109,22 @@ def __repr__(self):
     def __getitem__(self, i):
         return self.operator_list[i]
 
+    def __mul__(self, other):
+        if isinstance(other, (int, float)):
+            return Operators([op * other for op in self.operator_list])
+        else:
+            raise NotImplementedError("Multiplication only supported for scalars.")
+
+    def __imul__(self, other):
+        if isinstance(other, (int, float)):
+            self.operator_list = [op * other for op in self.operator_list]
+            return self
+        else:
+            raise NotImplementedError("Multiplication only supported for scalars.")
+
+    def __rmul__(self, other):
+        return self.__mul__(other)
+
     @property
     def nqbits(self):
         return self[0].nqbits
diff --git a/tests/qse/operator_test.py b/tests/qse/operator_test.py
index b4592e6..7f43c75 100644
--- a/tests/qse/operator_test.py
+++ b/tests/qse/operator_test.py
@@ -86,6 +86,27 @@ def test_operator_fail(qubits, operator):
         qse.Operator(operator, qubits, nqubits=4)
 
 
+@pytest.mark.parametrize("qubits", [-1, 4, [-1, 2], [3, 4]])
+def test_operator_fail_qubit_index(qubits):
+    """Test the class raises an error for a qubit outside the index range."""
+    with pytest.raises(Exception):
+        qse.Operator("Z", qubits, nqubits=4)
+
+
+def test_operator_mul():
+    """Test multiplying an Operator by a scalar."""
+    op = qse.Operator("X", 0, nqbits=2, coef=1.0)
+
+    op_scaled = op * 3.5
+
+    assert isinstance(op_scaled, qse.Operator)
+    assert op_scaled.to_str() == "XI"
+    assert np.isclose(op_scaled.coef, 3.5)
+
+    op *= 2.0
+    assert np.isclose(op.coef, 2.0)
+
+
 def test_operators_init():
     """Test initializing empty Operators."""
     ops = qse.Operators()
@@ -123,3 +144,18 @@ def test_operators_qutip():
 
     op_sum = op1.to_qutip() + op2.to_qutip()
     assert np.allclose(op_sum.full(), ops.to_qutip().full())
+
+
+def test_operators_mul():
+    """Test multiplying Operators by a scalar."""
+    nqbits = 4
+    ops = qse.Operators(
+        [qse.Operator("X", i, nqbits=nqbits + 5) for i in range(nqbits)]
+    )
+    assert all(np.isclose(op.coef, 1.0) for op in ops)
+
+    ops_scaled = ops * 2.5
+    assert all(np.isclose(op.coef, 2.5) for op in ops_scaled)
+
+    ops *= 3.5
+    assert all(np.isclose(op.coef, 3.5) for op in ops)
diff --git a/uv.lock b/uv.lock
index 8d23dc2..4fe71df 100644
--- a/uv.lock
+++ b/uv.lock
@@ -4005,7 +4005,7 @@ wheels = [
 
 [[package]]
 name = "qse"
-version = "1.1.14"
+version = "1.1.15"
 source = { editable = "." }
 dependencies = [
     { name = "matplotlib" },