diff --git a/pennylane_rigetti/numpy_wavefunction.py b/pennylane_rigetti/numpy_wavefunction.py
index 38eb5976..05bc4a39 100644
--- a/pennylane_rigetti/numpy_wavefunction.py
+++ b/pennylane_rigetti/numpy_wavefunction.py
@@ -41,7 +41,7 @@ class NumpyWavefunctionDevice(RigettiDevice):
     name = "pyQVM NumpyWavefunction Simulator Device"
     short_name = "rigetti.numpy_wavefunction"
 
-    observables = {"PauliX", "PauliY", "PauliZ", "Hadamard", "Hermitian", "Identity"}
+    observables = {"PauliX", "PauliY", "PauliZ", "Hadamard", "Hermitian", "Identity", "Prod"}
 
     def __init__(self, wires, *, shots=None):
         super().__init__(wires, shots)
diff --git a/pennylane_rigetti/qpu.py b/pennylane_rigetti/qpu.py
index 6993175b..12553c23 100644
--- a/pennylane_rigetti/qpu.py
+++ b/pennylane_rigetti/qpu.py
@@ -25,6 +25,7 @@
 import numpy as np
 from pennylane.measurements import Expectation
 from pennylane.operation import Tensor
+from pennylane.ops import Prod
 from pennylane.tape import QuantumTape
 from pyquil import get_qc
 from pyquil.api import QuantumComputer
@@ -141,7 +142,7 @@ def expval(self, observable, shot_range=None, bin_size=None):
                 pauli_obs = sZ(wire)
 
             # Multi-qubit observable
-            elif len(device_wires) > 1 and isinstance(observable, Tensor):
+            elif len(device_wires) > 1 and isinstance(observable, (Tensor, Prod)):
                 # All observables are rotated to be measured in the Z-basis, so we just need to
                 # check which wires exist in the observable, map them to physical qubits, and measure
                 # the product of PauliZ operators on those qubits
diff --git a/pennylane_rigetti/wavefunction.py b/pennylane_rigetti/wavefunction.py
index 4d402642..2d1dd3d9 100644
--- a/pennylane_rigetti/wavefunction.py
+++ b/pennylane_rigetti/wavefunction.py
@@ -60,7 +60,7 @@ class WavefunctionDevice(RigettiDevice):
     name = "Rigetti Wavefunction Simulator Device"
     short_name = "rigetti.wavefunction"
 
-    observables = {"PauliX", "PauliY", "PauliZ", "Hadamard", "Hermitian", "Identity"}
+    observables = {"PauliX", "PauliY", "PauliZ", "Hadamard", "Hermitian", "Identity", "Prod"}
 
     def __init__(self, wires, *, shots=None):
         super().__init__(wires, shots)
diff --git a/tests/test_converter.py b/tests/test_converter.py
index 9bef99cd..949a534f 100644
--- a/tests/test_converter.py
+++ b/tests/test_converter.py
@@ -1148,9 +1148,7 @@ def test_load_quil_file_via_entry_point(self):
         cur_dir = os.path.dirname(os.path.abspath(__file__))
 
         with OperationRecorder() as rec:
-            qml.from_quil_file(os.path.join(cur_dir, "simple_program.quil"))(
-                wires=range(5)
-            )
+            qml.from_quil_file(os.path.join(cur_dir, "simple_program.quil"))(wires=range(5))
 
         # The wires should be assigned as
         # 0  1  2  3  7