diff --git a/pennylane/measurements/probs.py b/pennylane/measurements/probs.py
index a9f4b5deb59..3fdbb778e8d 100644
--- a/pennylane/measurements/probs.py
+++ b/pennylane/measurements/probs.py
@@ -266,9 +266,12 @@ def process_density_matrix(self, density_matrix: TensorLike, wire_order: Wires):
         # Since we only care about the probabilities, we can simplify the task here by creating a 'pseudo-state' to carry the diagonal elements and reuse the process_state method
         # Introduce a small epsilon to avoid sqrt(0) to ensure absolute positivity.
         prob = qml.math.real(prob)
-        prob = qml.math.where(prob < 0, -prob, prob)
-        prob = qml.math.where(prob == 0, prob + np.finfo(float).eps, prob)
         p_state = qml.math.sqrt(prob)
+        precision_issue_indices = qml.math.isnan(p_state) & list(
+            qml.math.allclose(p, 0) for p in prob
+        )
+        p_state = qml.math.where(precision_issue_indices, 0, p_state)
+
         return self.process_state(p_state, wire_order)
 
     @staticmethod