diff --git a/benchmarks/scripts/expval_benchmark.py b/benchmarks/scripts/expval_benchmark.py
index 2c60948..97d5753 100644
--- a/benchmarks/scripts/expval_benchmark.py
+++ b/benchmarks/scripts/expval_benchmark.py
@@ -7,7 +7,7 @@
 import pytket
 import qiskit
 from qiskit import qasm2
-from qiskit.quantum_info import Operator, Statevector
+from qiskit.quantum_info import Operator, Statevector, SparsePauliOp
 from qiskit_aer import AerSimulator
 import numpy as np
 
@@ -207,6 +207,57 @@ def estimate_heavy_output_prob(
     return hop
 
 
+def generate_qaoa_observable(num_qubits):
+    """Generates the problem Hamiltonian as the observable for the QAOA
+    benchmarking circuits, based on the binary encoding described in
+    Franz G. Fuchs, Herman Øie Kolden, Niels Henrik Aase, and Giorgio
+    Sartor "Efficient encoding of the weighted MAX k-CUT on a quantum computer
+    using QAOA". (2020) arXiv 2009.01095 (https://arxiv.org/abs/2009.01095).
+    The weights of the edges between vertices and of the resulting unitary
+    evolution come from the 10-vertex Barabasi-Albert graph in Fig 4(c)
+    of the paper.
+    """
+    pauli_strings = []
+    # Weights of edges between vertices and of the resulting unitary evolution
+    weighted_edges = [
+        (0, 1, 6.720),
+        (0, 2, 3.246),
+        (1, 2, 6.462),
+        (1, 3, 3.386),
+        (1, 5, 5.014),
+        (1, 6, 6.596),
+        (2, 3, 8.579),
+        (2, 4, 0.62),
+        (2, 5, 0.708),
+        (2, 6, 2.275),
+        (2, 7, 5.0),
+        (2, 8, 4.034),
+        (3, 4, 4.987),
+        (3, 6, 1.089),
+        (3, 9, 2.961),
+        (4, 5, 1.134),
+        (4, 6, 6.865),
+        (5, 6, 8.184),
+        (5, 7, 9.459),
+        (6, 7, 2.268),
+        (6, 8, 8.197),
+        (6, 9, 1.212),
+        (7, 9, 4.265),
+        (8, 9, 1.690),
+    ]
+    for i, j, _ in weighted_edges:
+        # Start with identity string
+        pauli_string = ["I"] * num_qubits
+        # Place Z operators on the chosen qubits
+        pauli_string[i] = "Z"
+        pauli_string[j] = "Z"
+        # Convert to PauliSumOp
+        pauli_strings.append("".join(pauli_string))
+    coeffs = [weight for _, _, weight in weighted_edges]
+    observable = SparsePauliOp(pauli_strings, coeffs)
+    return observable
+
+
 def simulate_expvals(
     uncompiled_circuit: qiskit.QuantumCircuit,
     compiled_circuit: qiskit.QuantumCircuit,
@@ -235,8 +286,15 @@ def simulate_expvals(
 
     else:
         density_matrix = simulate_density_matrix(compiled_circuit)
-        obs_str = "Z" * compiled.num_qubits
-        observable = Operator.from_label(obs_str)
+        if circuit_name == "qaoa_barabasi_albert":
+            # observable is the problem Hamiltonian
+            num_qubits = compiled_circuit.num_qubits
+            observable = generate_qaoa_observable(num_qubits)
+            obs_str = "".join(("H_p = ", str(observable.to_sparse_list())))
+
+        else:
+            obs_str = "Z" * compiled_circuit.num_qubits
+            observable = Operator.from_label(obs_str)
 
         compiled_ev = np.real(density_matrix.expectation_value(observable))