diff --git a/src/qibo/hamiltonians/hamiltonians.py b/src/qibo/hamiltonians/hamiltonians.py
index 28a78d68fe..2f4ab2c0d6 100644
--- a/src/qibo/hamiltonians/hamiltonians.py
+++ b/src/qibo/hamiltonians/hamiltonians.py
@@ -163,6 +163,24 @@ def eye(self, dim: Optional[int] = None):
             dim = int(self.matrix.shape[0])
         return self.backend.cast(self.backend.matrices.I(dim), dtype=self.matrix.dtype)
 
+    def energy_fluctuation(self, state):
+        """
+        Evaluate energy fluctuation:
+
+        .. math::
+            \\Xi_{k}(\\mu) = \\sqrt{\\langle\\mu|\\hat{H}^2|\\mu\\rangle - \\langle\\mu|\\hat{H}|\\mu\\rangle^2} \\,
+
+        for a given state :math:`|\\mu\\rangle`.
+
+        Args:
+            state (np.ndarray): quantum state to be used to compute the energy fluctuation.
+        """
+        energy = self.expectation(state)
+        h = self.matrix
+        h2 = Hamiltonian(nqubits=self.nqubits, matrix=h @ h, backend=self.backend)
+        average_h2 = self.backend.calculate_expectation_state(h2, state, normalize=True)
+        return np.sqrt(np.abs(average_h2 - energy**2))
+
     def __add__(self, o):
         if isinstance(o, self.__class__):
             if self.nqubits != o.nqubits:
diff --git a/src/qibo/models/dbi/double_bracket.py b/src/qibo/models/dbi/double_bracket.py
index 9b7beedbba..976fcbc021 100644
--- a/src/qibo/models/dbi/double_bracket.py
+++ b/src/qibo/models/dbi/double_bracket.py
@@ -180,9 +180,15 @@ def loss(self, step: float, look_ahead: int = 1):
         return loss
 
     def energy_fluctuation(self, state):
-        """Evaluate energy fluctuations"""
-        energy = self.h.expectation(state)
-        h = self.h.matrix
-        h2 = Hamiltonian(nqubits=self.h.nqubits, matrix=h @ h, backend=self.backend)
-        average_h2 = self.backend.calculate_expectation_state(h2, state, normalize=True)
-        return np.sqrt(average_h2 - energy**2)
+        """
+        Evaluate energy fluctuation
+
+        .. math::
+            \\Xi_{k}(\\mu) = \\sqrt{\\langle\\mu|\\hat{H}^2|\\mu\\rangle - \\langle\\mu|\\hat{H}|\\mu\\rangle^2} \\,
+
+        for a given state :math:`|\\mu\\rangle`.
+
+        Args:
+            state (np.ndarray): quantum state to be used to compute the energy fluctuation with H.
+        """
+        return self.h.energy_fluctuation(state)
diff --git a/src/qibo/models/variational.py b/src/qibo/models/variational.py
index d031e58ae0..950b8f3a5c 100644
--- a/src/qibo/models/variational.py
+++ b/src/qibo/models/variational.py
@@ -132,13 +132,7 @@ def energy_fluctuation(self, state):
         Args:
             state (np.ndarray): quantum state to be used to compute the energy fluctuation with H.
         """
-        energy = self.hamiltonian.expectation(state)
-        h = self.hamiltonian.matrix
-        h2 = Hamiltonian(
-            nqubits=self.hamiltonian.nqubits, matrix=h @ h, backend=self.backend
-        )
-        average_h2 = self.backend.calculate_expectation_state(h2, state, normalize=True)
-        return np.sqrt(average_h2 - energy**2)
+        return self.hamiltonian.energy_fluctuation(state)
 
 
 class AAVQE: