diff --git a/src/gravity_cache.h b/src/gravity_cache.h
index 6d073db60171a9d30d6f397213849e4c3d5314a1..f6a848f25173ffacab29b0b589914efcece65900 100644
--- a/src/gravity_cache.h
+++ b/src/gravity_cache.h
@@ -253,9 +253,11 @@ __attribute__((always_inline)) INLINE static void gravity_cache_populate(
       dz = nearestf(dz, dim[2]);
     }
     const float r2 = dx * dx + dy * dy + dz * dz;
+    const float epsilon2 = epsilon[i] * epsilon[i];
 
     /* Check whether we can use the multipole instead of P-P */
-    use_mpole[i] = allow_mpole && gravity_M2P_accept(r_max2, theta_crit2, r2);
+    use_mpole[i] =
+        allow_mpole && gravity_M2P_accept(r_max2, theta_crit2, r2, epsilon2);
   }
 
 #ifdef SWIFT_DEBUG_CHECKS
@@ -435,8 +437,9 @@ gravity_cache_populate_all_mpole(const timebin_t max_active_bin,
       dz = nearestf(dz, dim[2]);
     }
     const float r2 = dx * dx + dy * dy + dz * dz;
+    const float epsilon2 = epsilon[i] * epsilon[i];
 
-    if (!gravity_M2P_accept(r_max2, theta_crit2, r2))
+    if (!gravity_M2P_accept(r_max2, theta_crit2, r2, epsilon2))
       error("Using m-pole where the test fails");
 #endif
   }
diff --git a/src/multipole.h b/src/multipole.h
index fe17a624491ac88a8d9061cba99b47c69d48bc74..041be93330b3c77676e6204d67a9afd4ebbd9563 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -2601,18 +2601,24 @@ __attribute__((always_inline, const)) INLINE static int gravity_M2L_accept(
  * We use the multipole acceptance criterion of Dehnen, 2002, JCoPh, Volume 179,
  * Issue 1, pp.27-42, equation 10.
  *
+ * We also additionally check that the distance between the particle and the
+ * multipole is larger than then softening length (here the distance at which
+ * the gravity becomes Newtonian again, not the Plummer-equivalent quantity).
+ *
  * @param r_max2 The square of the size of the multipole.
  * @param theta_crit2 The square of the critical opening angle.
  * @param r2 Square of the distance (periodically wrapped) between the
  * particle and the multipole.
+ * @param epsilon2 The square of the softening length of the particle.
  */
 __attribute__((always_inline, const)) INLINE static int gravity_M2P_accept(
-    const float r_max2, const float theta_crit2, const float r2) {
+    const float r_max2, const float theta_crit2, const float r2,
+    const float epsilon2) {
 
   // MATTHIEU: Make this mass-dependent ?
 
   /* Multipole acceptance criterion (Dehnen 2002, eq.10) */
-  return (r2 * theta_crit2 > r_max2);
+  return (r2 * theta_crit2 > r_max2) && (r2 > epsilon2);
 }
 
 #endif /* SWIFT_MULTIPOLE_H */
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 72d00f7fd9db0385856f7daae08f51e87bddfffc..27e811b0044886e305913b19db0e9d48585ba4c7 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -519,7 +519,7 @@ static INLINE void runner_dopair_grav_pm_full(
     const float theta_crit2 = e->gravity_properties->theta_crit2;
 
     /* Note: 1.1 to avoid FP rounding false-positives */
-    if (!gravity_M2P_accept(r_max2, theta_crit2 * 1.1, r2))
+    if (!gravity_M2P_accept(r_max2, theta_crit2 * 1.1, r2, h_i * h_i))
       error(
           "use_mpole[i] set when M2P accept fails CoM=[%e %e %e] pos=[%e %e "
           "%e], rmax=%e",
@@ -650,7 +650,7 @@ static INLINE void runner_dopair_grav_pm_truncated(
     const float theta_crit2 = e->gravity_properties->theta_crit2;
 
     /* 1.1 to avoid FP rounding false-positives */
-    if (!gravity_M2P_accept(r_max2, theta_crit2 * 1.1, r2))
+    if (!gravity_M2P_accept(r_max2, theta_crit2 * 1.1, r2, h_i * h_i))
       error(
           "use_mpole[i] set when M2P accept fails CoM=[%e %e %e] pos=[%e %e "
           "%e], rmax=%e",