diff --git a/src/cell.c b/src/cell.c
index 5128ef727488fb4fca330f79e26eff849ffbf5b3..f1f29027cd09cec08bfc52d286cbcfd320d96dc6 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -2707,7 +2707,7 @@ void cell_drift_all_multipoles(struct cell *c, const struct engine *e) {
 
   /* Drift the multipole */
   if (ti_current > ti_old_multipole)
-    gravity_drift(c->multipole, dt_drift, 0.f /*c->dx_max_gpart*/);
+    gravity_drift(c->multipole, dt_drift);
 
   /* Are we not in a leaf ? */
   if (c->split) {
@@ -2749,7 +2749,7 @@ void cell_drift_multipole(struct cell *c, const struct engine *e) {
     dt_drift = (ti_current - ti_old_multipole) * e->time_base;
 
   if (ti_current > ti_old_multipole)
-    gravity_drift(c->multipole, dt_drift, 0.f /*c->dx_max_gpart*/);
+    gravity_drift(c->multipole, dt_drift);
 
   /* Update the time of the last drift */
   c->ti_old_multipole = ti_current;
diff --git a/src/multipole.h b/src/multipole.h
index 51fa99dde3c54f8c3190bf564e934e2de1b3d1e4..080cf669a11b9959ee83a392e024f7f88bfeba3d 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -108,10 +108,16 @@ struct grav_tensor {
 
 struct multipole {
 
-  /* Bulk velocity */
+  /*! Bulk velocity */
   float vel[3];
 
-  /* 0th order terms */
+  /*! Maximal velocity along each axis of all #gpart */
+  float max_delta_vel[3];
+
+  /*! Minimal velocity along each axis of all #gpart */
+  float min_delta_vel[3];
+
+  /* 0th order term */
   float M_000;
 
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 0
@@ -215,14 +221,15 @@ INLINE static void gravity_reset(struct gravity_tensors *m) {
 /**
  * @brief Drifts a #multipole forward in time.
  *
+ * This uses a first-order approximation in time. We only move the CoM
+ * using the bulk velocity measured at the last rebuild.
+ *
  * @param m The #multipole.
  * @param dt The drift time-step.
- * @param x_diff The maximal distance moved by any particle since the last
- * rebuild.
  */
-INLINE static void gravity_drift(struct gravity_tensors *m, double dt,
-                                 float x_diff) {
+INLINE static void gravity_drift(struct gravity_tensors *m, double dt) {
 
+  /* Motion of the centre of mass */
   const double dx = m->m_pole.vel[0] * dt;
   const double dy = m->m_pole.vel[1] * dt;
   const double dz = m->m_pole.vel[2] * dt;
@@ -232,8 +239,35 @@ INLINE static void gravity_drift(struct gravity_tensors *m, double dt,
   m->CoM[1] += dy;
   m->CoM[2] += dz;
 
+#ifdef SWIFT_DEBUG_CHECKS
+  if(m->m_pole.vel[0] > m->m_pole.max_delta_vel[0]) 
+    error("Invalid maximal velocity");
+  if(m->m_pole.vel[0] < m->m_pole.min_delta_vel[0]) 
+    error("Invalid minimal velocity");
+  if(m->m_pole.vel[1] > m->m_pole.max_delta_vel[1]) 
+    error("Invalid maximal velocity");
+  if(m->m_pole.vel[1] < m->m_pole.min_delta_vel[1]) 
+    error("Invalid minimal velocity");
+  if(m->m_pole.vel[2] > m->m_pole.max_delta_vel[2]) 
+    error("Invalid maximal velocity");
+  if(m->m_pole.vel[2] < m->m_pole.min_delta_vel[2]) 
+    error("Invalid minimal velocity");
+#endif
+
+  /* Maximal distance covered by any particle */
+  float dv[3];
+  dv[0] = max(m->m_pole.max_delta_vel[0] - m->m_pole.vel[0],
+	      m->m_pole.vel[0] - m->m_pole.min_delta_vel[0]);
+  dv[1] = max(m->m_pole.max_delta_vel[1] - m->m_pole.vel[1],
+	      m->m_pole.vel[1] - m->m_pole.min_delta_vel[1]);
+  dv[2] = max(m->m_pole.max_delta_vel[2] - m->m_pole.vel[2],
+	      m->m_pole.vel[2] - m->m_pole.min_delta_vel[2]);
+
+  const float max_delta_vel = sqrt(dv[0] * dv[0] + dv[1] * dv[1] + dv[2] * dv[2]);
+  const float x_diff = max_delta_vel * dt;
+
   /* Conservative change in maximal radius containing all gpart */
-  m->r_max = m->r_max_rebuild + x_diff;
+  m->r_max += x_diff;
 }
 
 /**
@@ -470,16 +504,8 @@ INLINE static void gravity_multipole_print(const struct multipole *m) {
 INLINE static void gravity_multipole_add(struct multipole *ma,
                                          const struct multipole *mb) {
 
-  const float M_000 = ma->M_000 + mb->M_000;
-  const float inv_M_000 = 1.f / M_000;
-
-  /* Add the bulk velocities */
-  ma->vel[0] = (ma->vel[0] * ma->M_000 + mb->vel[0] * mb->M_000) * inv_M_000;
-  ma->vel[1] = (ma->vel[1] * ma->M_000 + mb->vel[1] * mb->M_000) * inv_M_000;
-  ma->vel[2] = (ma->vel[2] * ma->M_000 + mb->vel[2] * mb->M_000) * inv_M_000;
-
-  /* Add 0th order terms */
-  ma->M_000 = M_000;
+  /* Add 0th order term */
+  ma->M_000 += mb->M_000;
 
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 0
   /* Add 1st order terms */
@@ -555,11 +581,6 @@ INLINE static void gravity_multipole_add(struct multipole *ma,
 #error "Missing implementation for order >5"
 #endif
 
-  // MATTHIEU
-  ma->M_100 = 0.f;
-  ma->M_010 = 0.f;
-  ma->M_001 = 0.f;
-
 #ifdef SWIFT_DEBUG_CHECKS
   ma->num_gpart += mb->num_gpart;
 #endif
@@ -1026,6 +1047,9 @@ INLINE static void gravity_P2M(struct gravity_tensors *multi,
 
   /* Prepare some local counters */
   double r_max2 = 0.;
+  float max_delta_vel[3] = {0., 0., 0.};
+  float min_delta_vel[3] = {0., 0., 0.};
+
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 0
   double M_100 = 0., M_010 = 0., M_001 = 0.;
 #endif
@@ -1068,6 +1092,16 @@ INLINE static void gravity_P2M(struct gravity_tensors *multi,
     /* Maximal distance CoM<->gpart */
     r_max2 = max(r_max2, dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
 
+    /* Store the vector of the maximal vel difference */
+    max_delta_vel[0] = max(gparts[k].v_full[0], max_delta_vel[0]);
+    max_delta_vel[1] = max(gparts[k].v_full[1], max_delta_vel[1]);
+    max_delta_vel[2] = max(gparts[k].v_full[2], max_delta_vel[2]);
+
+    /* Store the vector of the minimal vel difference */
+    min_delta_vel[0] = min(gparts[k].v_full[0], min_delta_vel[0]);
+    min_delta_vel[1] = min(gparts[k].v_full[1], min_delta_vel[1]);
+    min_delta_vel[2] = min(gparts[k].v_full[2], min_delta_vel[2]);
+
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 0
     const double m = gparts[k].mass;
 
@@ -1150,7 +1184,9 @@ INLINE static void gravity_P2M(struct gravity_tensors *multi,
   }
 
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 0
-  M_100 = M_010 = M_001 = 0.f; /* Matthieu */
+
+  /* We know the first-order multipole (dipole) is 0. */
+  M_100 = M_010 = M_001 = 0.f;
 #endif
 
   /* Store the data on the multipole. */
@@ -1162,6 +1198,12 @@ INLINE static void gravity_P2M(struct gravity_tensors *multi,
   multi->m_pole.vel[0] = vel[0];
   multi->m_pole.vel[1] = vel[1];
   multi->m_pole.vel[2] = vel[2];
+  multi->m_pole.max_delta_vel[0] = max_delta_vel[0];
+  multi->m_pole.max_delta_vel[1] = max_delta_vel[1];
+  multi->m_pole.max_delta_vel[2] = max_delta_vel[2];
+  multi->m_pole.min_delta_vel[0] = min_delta_vel[0];
+  multi->m_pole.min_delta_vel[1] = min_delta_vel[1];
+  multi->m_pole.min_delta_vel[2] = min_delta_vel[2];
 
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 0
 
@@ -1260,10 +1302,6 @@ INLINE static void gravity_P2M(struct gravity_tensors *multi,
 INLINE static void gravity_M2M(struct multipole *m_a,
                                const struct multipole *m_b,
                                const double pos_a[3], const double pos_b[3]) {
-  /* Shift bulk velocity */
-  m_a->vel[0] = m_b->vel[0];
-  m_a->vel[1] = m_b->vel[1];
-  m_a->vel[2] = m_b->vel[2];
 
   /* Shift 0th order term */
   m_a->M_000 = m_b->M_000;
diff --git a/src/space.c b/src/space.c
index a4c0d9b595e979381ea02de5a4ea3e466f5e4163..fa4dbba7082d3a3be603522a44b7fcc22523cd0e 100644
--- a/src/space.c
+++ b/src/space.c
@@ -1886,22 +1886,53 @@ void space_split_recursive(struct space *s, struct cell *c,
       /* Reset everything */
       gravity_reset(c->multipole);
 
-      /* Compute CoM of all progenies */
+      /* Compute CoM and bulk velocity from all progenies */
       double CoM[3] = {0., 0., 0.};
+      double vel[3] = {0., 0., 0.};
+      float max_delta_vel[3] = {0.f, 0.f, 0.f};
+      float min_delta_vel[3] = {0.f, 0.f, 0.f};
       double mass = 0.;
 
       for (int k = 0; k < 8; ++k) {
         if (c->progeny[k] != NULL) {
           const struct gravity_tensors *m = c->progeny[k]->multipole;
+
+          mass += m->m_pole.M_000;
+
           CoM[0] += m->CoM[0] * m->m_pole.M_000;
           CoM[1] += m->CoM[1] * m->m_pole.M_000;
           CoM[2] += m->CoM[2] * m->m_pole.M_000;
-          mass += m->m_pole.M_000;
+
+	  vel[0] += m->m_pole.vel[0] * m->m_pole.M_000;
+	  vel[1] += m->m_pole.vel[1] * m->m_pole.M_000;	  
+	  vel[2] += m->m_pole.vel[2] * m->m_pole.M_000;
+
+	  max_delta_vel[0] = max(m->m_pole.max_delta_vel[0], max_delta_vel[0]);
+	  max_delta_vel[1] = max(m->m_pole.max_delta_vel[1], max_delta_vel[1]);
+	  max_delta_vel[2] = max(m->m_pole.max_delta_vel[2], max_delta_vel[2]);
+
+	  min_delta_vel[0] = min(m->m_pole.min_delta_vel[0], min_delta_vel[0]);
+	  min_delta_vel[1] = min(m->m_pole.min_delta_vel[1], min_delta_vel[1]);
+	  min_delta_vel[2] = min(m->m_pole.min_delta_vel[2], min_delta_vel[2]);
         }
       }
-      c->multipole->CoM[0] = CoM[0] / mass;
-      c->multipole->CoM[1] = CoM[1] / mass;
-      c->multipole->CoM[2] = CoM[2] / mass;
+
+      /* Final operation on the CoM and bulk velocity */
+      const double inv_mass = 1. / mass;
+      c->multipole->CoM[0] = CoM[0] * inv_mass;
+      c->multipole->CoM[1] = CoM[1] * inv_mass;
+      c->multipole->CoM[2] = CoM[2] * inv_mass;
+      c->multipole->m_pole.vel[0] = vel[0] * inv_mass;
+      c->multipole->m_pole.vel[1] = vel[1] * inv_mass;
+      c->multipole->m_pole.vel[2] = vel[2] * inv_mass;
+
+      /* Min max velocity along each axis */
+      c->multipole->m_pole.max_delta_vel[0] = max_delta_vel[0];
+      c->multipole->m_pole.max_delta_vel[1] = max_delta_vel[1];
+      c->multipole->m_pole.max_delta_vel[2] = max_delta_vel[2];
+      c->multipole->m_pole.min_delta_vel[0] = min_delta_vel[0];
+      c->multipole->m_pole.min_delta_vel[1] = min_delta_vel[1];
+      c->multipole->m_pole.min_delta_vel[2] = min_delta_vel[2];
 
       /* Now shift progeny multipoles and add them up */
       struct multipole temp;
@@ -1923,6 +1954,7 @@ void space_split_recursive(struct space *s, struct cell *c,
           r_max = max(r_max, cp->multipole->r_max + sqrt(r2));
         }
       }
+
       /* Alternative upper limit of max CoM<->gpart distance */
       const double dx = c->multipole->CoM[0] > c->loc[0] + c->width[0] / 2.
                             ? c->multipole->CoM[0] - c->loc[0]
@@ -1943,6 +1975,11 @@ void space_split_recursive(struct space *s, struct cell *c,
       c->multipole->CoM_rebuild[1] = c->multipole->CoM[1];
       c->multipole->CoM_rebuild[2] = c->multipole->CoM[2];
 
+      /* We know the first-order multipole (dipole) is 0. */
+      c->multipole->m_pole.M_100 = 0.f;
+      c->multipole->m_pole.M_010 = 0.f;
+      c->multipole->m_pole.M_001 = 0.f;
+
     } /* Deal with gravity */
   }   /* Split or let it be? */