diff --git a/pulser-simulation/pulser_simulation/simresults.py b/pulser-simulation/pulser_simulation/simresults.py
index 387916554..8e954c543 100644
--- a/pulser-simulation/pulser_simulation/simresults.py
+++ b/pulser-simulation/pulser_simulation/simresults.py
@@ -521,28 +521,34 @@ def sample_state(
             quantum states at time t.
         """
         sampled_state = super().sample_state(t, n_samples, t_tol)
-        if self._meas_errors is None:
+        if self._meas_errors is None or (
+            self._meas_errors["epsilon"] == 0.0
+            and self._meas_errors["epsilon_prime"] == 0
+        ):
             return sampled_state
 
-        detected_sample_dict: Counter = Counter()
-        for shot, n_detects in sampled_state.items():
-            eps = self._meas_errors["epsilon"]
-            eps_p = self._meas_errors["epsilon_prime"]
-            # Shot as an array of 1s and 0s
-            shot_arr = np.array(list(shot), dtype=int)
-            # Probability of flipping each bit
-            flip_probs = np.array([eps_p if x == "1" else eps for x in shot])
-            # 1 if it flips, 0 if it stays the same
-            flips = (
-                np.random.uniform(size=(n_detects, len(flip_probs)))
-                < flip_probs
-            ).astype(int)
-            # XOR betwen the original array and the flips
-            # Gives an array of n_detects individual shots
-            new_shots = shot_arr ^ flips
-            # Count all the new_shots
-            detected_sample_dict += Counter(
-                "".join(map(str, measured)) for measured in new_shots
-            )
-
-        return detected_sample_dict
+        eps = self._meas_errors["epsilon"]
+        eps_p = self._meas_errors["epsilon_prime"]
+        shots = list(sampled_state.keys())
+        n_detects_list = list(sampled_state.values())
+
+        # Convert shots to a 2D array
+        shot_arr = np.array([list(shot) for shot in shots], dtype=int)
+        # Compute flip probabilities
+        flip_probs = np.where(shot_arr == 1, eps_p, eps)
+        # Repeat flip_probs based on n_detects_list
+        flip_probs_repeated = np.repeat(flip_probs, n_detects_list, axis=0)
+        # Generate random matrix of shape (sum(n_detects_list), len(shot))
+        random_matrix = np.random.uniform(
+            size=(np.sum(n_detects_list), len(shot_arr[0]))
+        )
+        # Compare random matrix with flip probabilities
+        flips = random_matrix < flip_probs_repeated
+        # Perform XOR between original array and flips
+        new_shots = shot_arr.repeat(n_detects_list, axis=0) ^ flips
+        # Count all the new_shots
+        # We are not converting to str before because tuple indexing is faster
+        detected_sample_dict: Counter = Counter(map(tuple, new_shots))
+        return Counter(
+            {"".join(map(str, k)): v for k, v in detected_sample_dict.items()}
+        )
diff --git a/tests/test_simulation.py b/tests/test_simulation.py
index f08950943..3959376a3 100644
--- a/tests/test_simulation.py
+++ b/tests/test_simulation.py
@@ -643,6 +643,26 @@ def test_noise(seq, matrices):
                 assert np.all(sim2.samples["Local"][basis][t][qty] == 0.0)
 
 
+def test_noise_with_zero_epsilons(seq, matrices):
+    np.random.seed(3)
+    sim = Simulation(seq, sampling_rate=0.01)
+
+    sim2 = Simulation(
+        seq,
+        sampling_rate=0.01,
+        config=SimConfig(
+            noise=("SPAM"), eta=0.0, epsilon=0.0, epsilon_prime=0.0
+        ),
+    )
+    assert sim2.config.spam_dict == {
+        "eta": 0,
+        "epsilon": 0.0,
+        "epsilon_prime": 0.0,
+    }
+
+    assert sim.run().sample_final_state() == sim2.run().sample_final_state()
+
+
 def test_dephasing():
     np.random.seed(123)
     reg = Register.from_coordinates([(0, 0)], prefix="q")