diff --git a/src/hydro/AnarchyDU/hydro_iact.h b/src/hydro/AnarchyDU/hydro_iact.h
index 27353e78bac4f848a5d7d1c53ff62648613f1b2d..cba945ecae8cbf2971b3d0d810cd565cb8ecc8eb 100644
--- a/src/hydro/AnarchyDU/hydro_iact.h
+++ b/src/hydro/AnarchyDU/hydro_iact.h
@@ -397,12 +397,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   const float visc_du_term = 0.5f * visc_acc_term * dvdr_Hubble;
 
   /* Diffusion term */
-  const float v_diff =
-      max(pi->force.soundspeed + pj->force.soundspeed + dvdr_Hubble, 0.f);
   const float alpha_diff = 0.5f * (pi->diffusion.alpha + pj->diffusion.alpha);
+  const float v_diff =
+      alpha_diff * sqrtf(0.5f * fabsf(pressurei - pressurej) / rho_ij) +
+      fabsf(fac_mu * r_inv * dvdr_Hubble);
   /* wi_dx + wj_dx / 2 is F_ij */
   const float diff_du_term =
-      alpha_diff * fac_mu * v_diff * (pi->u - pj->u) * (wi_dr + wj_dr) / rho_ij;
+      v_diff * (pi->u - pj->u) * (wi_dr / rhoi + wj_dr / rhoj);
 
   /* Assemble the energy equation term */
   const float du_dt_i = sph_du_term_i + visc_du_term + diff_du_term;
@@ -518,12 +519,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   const float visc_du_term = 0.5f * visc_acc_term * dvdr_Hubble;
 
   /* Diffusion term */
-  const float v_diff =
-      max(pi->force.soundspeed + pj->force.soundspeed + dvdr_Hubble, 0.f);
   const float alpha_diff = 0.5f * (pi->diffusion.alpha + pj->diffusion.alpha);
+  const float v_diff =
+      alpha_diff * sqrtf(0.5f * fabsf(pressurei - pressurej) / rho_ij) +
+      fabsf(fac_mu * r_inv * dvdr_Hubble);
   /* wi_dx + wj_dx / 2 is F_ij */
   const float diff_du_term =
-      alpha_diff * fac_mu * v_diff * (pi->u - pj->u) * (wi_dr + wj_dr) / rho_ij;
+      v_diff * (pi->u - pj->u) * (wi_dr / rhoi + wj_dr / rhoj);
 
   /* Assemble the energy equation term */
   const float du_dt_i = sph_du_term_i + visc_du_term + diff_du_term;
diff --git a/src/hydro/AnarchyPU/hydro_iact.h b/src/hydro/AnarchyPU/hydro_iact.h
index 106fbdbcd42b4b7f2f837d542496551445c4b752..bb53221908d71535d090e8d2581e466fcf1009f9 100644
--- a/src/hydro/AnarchyPU/hydro_iact.h
+++ b/src/hydro/AnarchyPU/hydro_iact.h
@@ -414,7 +414,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
 
   /* Diffusion term */
   const float v_diff =
-      max(pi->force.soundspeed + pj->force.soundspeed + dvdr_Hubble, 0.f);
+      max(pi->force.soundspeed + pj->force.soundspeed + mu_ij, 0.f);
   const float alpha_diff = 0.5f * (pi->diffusion.alpha + pj->diffusion.alpha);
   /* wi_dx + wj_dx / 2 is F_ij */
   const float diff_du_term =
@@ -538,7 +538,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
 
   /* Diffusion term */
   const float v_diff =
-      max(pi->force.soundspeed + pj->force.soundspeed + dvdr_Hubble, 0.f);
+      max(pi->force.soundspeed + pj->force.soundspeed + mu_ij, 0.f);
   const float alpha_diff = 0.5f * (pi->diffusion.alpha + pj->diffusion.alpha);
   /* wi_dx + wj_dx / 2 is F_ij */
   const float diff_du_term =