diff --git a/src/qibotn/eval_qu.py b/src/qibotn/eval_qu.py
index 302a62c0..457a6824 100644
--- a/src/qibotn/eval_qu.py
+++ b/src/qibotn/eval_qu.py
@@ -65,8 +65,7 @@ def expectation_qu(
     # use cotengra package for tensor contractions
     import cotengra as ctg
 
-    # used to add one more qubit to the circuit to make the observable local
-    qasm_mod, nqubits = modify_qasm(qasm)
+    nqubits = int(np.log2(len(initial_state)))
 
     if initial_state is not None:
         initial_state = init_state_tn(nqubits, initial_state)
@@ -77,7 +76,7 @@ def expectation_qu(
     )
 
     # generates the global observable
-    obs = pauli_string_gen(nqubits - 1, pauli_string_pattern)
+    obs = pauli_string_gen(nqubits, pauli_string_pattern)
 
     # parameters to find the contraction path using cotengra
     opt = ctg.ReusableHyperOptimizer(
@@ -99,41 +98,12 @@ def expectation_qu(
 
     # expectation value
     expectation = circ_quimb.local_expectation(
-        obs, where=list(range(nqubits - 1)), optimize=opt, simplify_sequence="DRC"
+        obs, where=list(range(nqubits)), optimize=opt, simplify_sequence="DRC"
     )
 
     return expectation
 
 
-def modify_qasm(qasm_circ):
-    """Generate a modified qasm string.
-
-    Args:
-        qasm (str): QASM program.
-
-    Returns:
-        string: QASM program with an additional auxillary qubit for the calculation of expectation
-    """
-
-    import re
-
-    lines = qasm_circ.split("\n")
-
-    qasm_circ_mod = []
-    while lines:
-        line = lines.pop(0).strip()
-        sta = re.compile(r"qreg\s+(\w+)\s*\[(\d+)\];")
-        match = sta.match(line)
-        if match:
-            name, nqubits = match.groups()
-            qasm_circ_mod.append(f"qreg q[{int(nqubits)+1}];")
-        else:
-            qasm_circ_mod.append(line)
-    qasm_circ_mod = "\n".join(qasm_circ_mod)
-
-    return qasm_circ_mod, int(nqubits) + 1
-
-
 def pauli_string_gen(nqubits, pauli_string_pattern):
     """Used internally to generate the string based on given pattern and number
     of qubit.
diff --git a/tests/test_quimb_backend.py b/tests/test_quimb_backend.py
index e32aefe5..d72d5647 100644
--- a/tests/test_quimb_backend.py
+++ b/tests/test_quimb_backend.py
@@ -64,3 +64,71 @@ def test_eval(nqubits: int, tolerance: float, is_mps: bool):
     assert np.allclose(
         result_sv, result_tn, atol=tolerance
     ), "Resulting dense vectors do not match"
+
+
+@pytest.mark.parametrize(
+    "nqubits, tolerance, is_mps",
+    [(1, 1e-6, True), (2, 1e-6, False), (5, 1e-3, True), (10, 1e-3, False)],
+)
+def test_eval_expectation(nqubits: int, tolerance: float, is_mps: bool):
+    """Evaluate circuit with Quimb backend and calculate the expecatation value
+    of the given Pauli string.
+
+    Args:
+        nqubits (int): Total number of qubits in the system.
+        tolerance (float): Maximum limit allowed for difference in results
+        is_mps (bool): True if state is MPS and False for tensor network structure
+    """
+    # hack quimb to use the correct number of processes
+    # TODO: remove completely, or at least delegate to the backend
+    # implementation
+    os.environ["QUIMB_NUM_PROCS"] = str(os.cpu_count())
+    import qibotn.eval_qu
+
+    init_state = create_init_state(nqubits=nqubits)
+    init_state_tn = copy.deepcopy(init_state)
+
+    # Test qibo
+    from qibo.hamiltonians import SymbolicHamiltonian
+    from qibo.symbols import X, Y
+
+    qibo.set_backend(backend=config.qibo.backend, platform=config.qibo.platform)
+    qibo_circ, state_vec = qibo_qft(nqubits, init_state, swaps=True)
+
+    # creating the Pauli string observable
+
+    list_of_objects = []
+    for i in range(nqubits):
+        if i % 4 == 0:
+            list_of_objects.append(X(i))
+        elif i % 4 == 1:
+            list_of_objects.append(X(i))
+        elif i % 4 == 2:
+            list_of_objects.append(X(i))
+        else:
+            list_of_objects.append(Y(i))
+    obs = np.prod(list_of_objects)
+    obs = SymbolicHamiltonian(obs)
+
+    # Noise-free expected value
+    qibo_exp = obs.expectation(state_vec)
+
+    # Convert to qasm for other backends
+    qasm_circ = qibo_circ.to_qasm()
+
+    # Test quimb
+    if is_mps:
+        gate_opt = {}
+        gate_opt["method"] = "svd"
+        gate_opt["cutoff"] = 1e-6
+        gate_opt["cutoff_mode"] = "abs"
+    else:
+        gate_opt = None
+
+    qibotn_exp = qibotn.eval_qu.expectation_qu(
+        qasm_circ, "XXXY", init_state_tn, gate_opt, backend=config.quimb.backend
+    )
+
+    assert np.allclose(
+        qibo_exp, qibotn_exp, atol=tolerance
+    ), "Resulting dense vectors do not match"