diff --git a/notebooks/Chemistry usage.ipynb b/notebooks/Chemistry usage.ipynb
index fd6ebed..2a9d9f5 100644
--- a/notebooks/Chemistry usage.ipynb	
+++ b/notebooks/Chemistry usage.ipynb	
@@ -6276,7 +6276,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3.8.10 64-bit",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -6290,7 +6290,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.10"
+   "version": "3.9.7"
   },
   "vscode": {
    "interpreter": {
diff --git a/symmer/symplectic/base.py b/symmer/symplectic/base.py
index 13ec579..0f933c8 100644
--- a/symmer/symplectic/base.py
+++ b/symmer/symplectic/base.py
@@ -11,6 +11,7 @@
 from openfermion import QubitOperator
 from qiskit.quantum_info import Pauli
 from qiskit.opflow import PauliOp, PauliSumOp
+from scipy.stats import unitary_group
 import warnings
 warnings.simplefilter('always', UserWarning)
 
@@ -63,6 +64,22 @@ def random(cls,
             coeff_vec += 1j * np.random.randn(n_terms)
         return cls(symp_matrix, coeff_vec)
 
+    @classmethod
+    def haar_random(cls,
+            n_qubits: int,
+        ) -> "PauliwordOp":
+        """ Generate a Haar random U(N) matrix (N^n_qubits) as a linear combination of Pauli operators.
+        aka generate a uniform random unitary from a Hilbert space.
+
+        Args:
+            n_qubits: number of qubits
+        Returns:
+            p_random (PauliwordOp): Haar random matrix in Pauli basis
+        """
+        haar_matrix = unitary_group.rvs(2**n_qubits)
+        p_random = cls.from_matrix(haar_matrix)
+        return p_random
+
     @classmethod
     def from_list(cls, 
             pauli_terms :List[str], 
@@ -846,6 +863,31 @@ def copy(self) -> "QuantumState":
         """
         return deepcopy(self)
 
+    @classmethod
+    def haar_random(cls,
+                    n_qubits: int,
+                    vec_type: str) -> "QuantumState":
+        """
+        Generate a Haar random quantum state - (uniform random quantum state).
+
+        Args:
+            n_qubits: number of qubits
+            vec_type (str): bra or ket
+
+        Returns:
+            qstate_random (QuantumState): Haar random quantum state
+
+        """
+        if vec_type=='ket':
+            haar_vec = (unitary_group.rvs(2**n_qubits)[:,0]).reshape([-1, 1])
+        elif vec_type == 'bra':
+            haar_vec = (unitary_group.rvs(2**n_qubits)[0,:]).reshape([1, -1])
+        else:
+            raise ValueError(f'vector type: {vec_type} unkown')
+
+        qstate_random = cls.from_array(haar_vec)
+        return qstate_random
+
     def __str__(self) -> str:
         """ 
         Defines the print behaviour of QuantumState - differs depending on vec_type