diff --git a/src/qibo/hamiltonians/hamiltonians.py b/src/qibo/hamiltonians/hamiltonians.py
index 28a78d68fe..0aa68d08a6 100644
--- a/src/qibo/hamiltonians/hamiltonians.py
+++ b/src/qibo/hamiltonians/hamiltonians.py
@@ -163,6 +163,27 @@ 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.
+
+        Return:
+            Energy fluctuation value (float).
+        """
+        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:
diff --git a/tests/test_hamiltonians.py b/tests/test_hamiltonians.py
index 960dbc797d..19a22437d8 100644
--- a/tests/test_hamiltonians.py
+++ b/tests/test_hamiltonians.py
@@ -423,3 +423,18 @@ def construct_hamiltonian():
 
     backend.assert_allclose(H.exp(0.5), target_matrix)
     backend.assert_allclose(H1.exp(0.5), target_matrix)
+
+
+def test_hamiltonian_energy_fluctuation(backend):
+    """Test energy fluctuation."""
+    # define hamiltonian
+    ham = hamiltonians.XXZ(nqubits=2, backend=backend)
+    # take ground state and zero state
+    ground_state = ham.ground_state()
+    zero_state = np.ones(2**2) / np.sqrt(2**2)
+    # collect energy fluctuations
+    gs_energy_fluctuation = ham.energy_fluctuation(ground_state)
+    zs_energy_fluctuation = ham.energy_fluctuation(zero_state)
+
+    assert np.isclose(gs_energy_fluctuation, 0, atol=1e-5)
+    assert gs_energy_fluctuation < zs_energy_fluctuation