diff --git a/.jules/bolt.md b/.jules/bolt.md
index 96776c7..63e8323 100644
--- a/.jules/bolt.md
+++ b/.jules/bolt.md
@@ -5,3 +5,7 @@
 ## 2025-11-01 - [Avoid Temporary Array Allocations in Reductions]
 **Learning:** In loops over parameters like temperature (`thermal_scan`), doing `np.sum(arr1 * arr2)` where array sizes match the system Hilbert space creates huge temporary arrays per loop iteration. Memory allocation dominates execution time.
 **Action:** Use `np.dot(arr1, arr2)` instead of `np.sum(arr1 * arr2)` for 1D arrays to evaluate reductions in C-level BLAS/NumPy functions, bypassing Python/NumPy array allocations and boosting performance significantly with less peak memory.
+
+## 2025-11-01 - [O(N^3) Matrix Multiplication Avoidance in Susceptibility]
+**Learning:** In computing `np.sum(rho * np.real(np.diag(A @ A)))` for static susceptibilities where `A` is Hermitian, performing the full matrix multiplication `A @ A` costs $O(N^3)$ and dominates execution time for large N.
+**Action:** Since `A` is Hermitian, `diag(A@A)_i = sum_j |A_ij|^2`. This allows replacing the $O(N^3)$ matrix multiplication with an $O(N^2)$ operation: `np.sum(np.abs(A)**2, axis=1)`. Combining this with `np.dot` instead of `np.sum(rho * arr)` avoids temporary allocations and massive speedups for large matrices.
diff --git a/physics/response/susceptibility.py b/physics/response/susceptibility.py
index de715ad..acad0ba 100644
--- a/physics/response/susceptibility.py
+++ b/physics/response/susceptibility.py
@@ -250,11 +250,12 @@ def static_susceptibility(
         rho     /= Z
 
         # <A> = \sum _n rho _n <n|A|n>
-        A_avg   = np.sum(rho * np.real(np.diag(A_q_eigen)))
+        A_avg = np.dot(rho, np.real(np.diag(A_q_eigen)))
 
         # <A^2> = \sum _n rho _n <n|A^2|n>
-        A2_eigen    = A_q_eigen @ A_q_eigen
-        A2_avg      = np.sum(rho * np.real(np.diag(A2_eigen)))
+        # Optimized: A is Hermitian, so diag(A@A)_i = sum_j |A_ij|^2. Replaces O(N^3) with O(N^2).
+        A2_diag = np.sum(np.abs(A_q_eigen)**2, axis=1)
+        A2_avg = np.dot(rho, A2_diag)
 
         # chi = \beta (<A^2> - <A>^2)
         chi_static  = beta * (A2_avg - A_avg**2)