diff --git a/.gitignore b/.gitignore
index 2a16a4a681b86d30940001fb2ce86dfb9d27a558..ba45fb06ea5e5defc3bf3592922b44d036afecd6 100644
--- a/.gitignore
+++ b/.gitignore
@@ -88,6 +88,8 @@ theory/SPH/Flavours/sph_flavours.pdf
 theory/SPH/EoS/eos.pdf
 theory/SPH/*.pdf
 theory/paper_pasc/pasc_paper.pdf
+theory/Multipoles/fmm.pdf
+theory/Multipoles/fmm_standalone.pdf
 
 m4/libtool.m4
 m4/ltoptions.m4
diff --git a/src/cell.c b/src/cell.c
index 874f03f0cad3257d866b70ec63fd8a6bbcd4b6a7..d4cea2b09ef9dc5a102e46dd216a17dc34b604c4 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1017,14 +1017,15 @@ void cell_clean_links(struct cell *c, void *data) {
 }
 
 /**
- * @brief Checks that a cell is at the current point in time
+ * @brief Checks that the particles in a cell are at the
+ * current point in time
  *
  * Calls error() if the cell is not at the current time.
  *
  * @param c Cell to act upon
  * @param data The current time on the integer time-line
  */
-void cell_check_drift_point(struct cell *c, void *data) {
+void cell_check_particle_drift_point(struct cell *c, void *data) {
 
 #ifdef SWIFT_DEBUG_CHECKS
 
@@ -1056,6 +1057,34 @@ void cell_check_drift_point(struct cell *c, void *data) {
 #endif
 }
 
+/**
+ * @brief Checks that the multipole of a cell is at the current point in time
+ *
+ * Calls error() if the cell is not at the current time.
+ *
+ * @param c Cell to act upon
+ * @param data The current time on the integer time-line
+ */
+void cell_check_multipole_drift_point(struct cell *c, void *data) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+
+  const integertime_t ti_drift = *(integertime_t *)data;
+
+  /* Only check local cells */
+  if (c->nodeID != engine_rank) return;
+
+  if (c->ti_old_multipole != ti_drift)
+    error(
+        "Cell multipole in an incorrect time-zone! c->ti_old_multipole=%lld "
+        "ti_drift=%lld (depth=%d)",
+        c->ti_old_multipole, ti_drift, c->depth);
+
+#else
+  error("Calling debugging code without debugging flag activated.");
+#endif
+}
+
 /**
  * @brief Resets all the individual cell task counters to 0.
  *
@@ -1128,7 +1157,8 @@ void cell_check_multipole(struct cell *c, void *data) {
 
     /* Now  compare the multipole expansion */
     if (!gravity_multipole_equal(&ma, c->multipole, tolerance)) {
-      message("Multipoles are not equal at depth=%d!", c->depth);
+      message("Multipoles are not equal at depth=%d! tol=%f", c->depth,
+              tolerance);
       message("Correct answer:");
       gravity_multipole_print(&ma.m_pole);
       message("Recursive multipole:");
@@ -1476,17 +1506,15 @@ void cell_drift_all_multipoles(struct cell *c, const struct engine *e) {
   /* Check that we are actually going to move forward. */
   if (ti_current < ti_old_multipole) error("Attempt to drift to the past");
 
+  /* Drift the multipole */
+  if (ti_current > ti_old_multipole) gravity_drift(c->multipole, dt);
+
   /* Are we not in a leaf ? */
   if (c->split) {
 
-    /* Loop over the progeny and drift the multipoles. */
+    /* Loop over the progeny and recurse. */
     for (int k = 0; k < 8; k++)
-      if (c->progeny[k] != NULL) cell_drift_particles(c->progeny[k], e);
-
-  } else if (ti_current > ti_old_multipole) {
-
-    /* Drift the multipole */
-    gravity_drift(c->multipole, dt);
+      if (c->progeny[k] != NULL) cell_drift_all_multipoles(c->progeny[k], e);
   }
 
   /* Update the time of the last drift */
diff --git a/src/cell.h b/src/cell.h
index b9fbb6abffe427b81daec033447fbf35804c242d..113ed28fe59bb7c6caa667b785ec72605d434fa9 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -352,7 +352,8 @@ int cell_are_neighbours(const struct cell *restrict ci,
                         const struct cell *restrict cj);
 void cell_check_multipole(struct cell *c, void *data);
 void cell_clean(struct cell *c);
-void cell_check_drift_point(struct cell *c, void *data);
+void cell_check_particle_drift_point(struct cell *c, void *data);
+void cell_check_multipole_drift_point(struct cell *c, void *data);
 void cell_reset_task_counters(struct cell *c);
 int cell_is_drift_needed(struct cell *c, const struct engine *e);
 int cell_unskip_tasks(struct cell *c, struct scheduler *s);
diff --git a/src/engine.c b/src/engine.c
index 0814f5e955772c18a65116b6ac2aba97736f7e37..c8ce267982adbf67349fb9fc81fe08989b0b12a4 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -846,7 +846,8 @@ void engine_repartition(struct engine *e) {
   fflush(stdout);
 
   /* Check that all cells have been drifted to the current time */
-  space_check_drift_point(e->s, e->ti_old);
+  space_check_drift_point(e->s, e->ti_old,
+                          e->policy & engine_policy_self_gravity);
 #endif
 
   /* Clear the repartition flag. */
@@ -2664,7 +2665,8 @@ void engine_rebuild(struct engine *e) {
   /* Check that all cells have been drifted to the current time.
    * That can include cells that have not
    * previously been active on this rank. */
-  space_check_drift_point(e->s, e->ti_old);
+  space_check_drift_point(e->s, e->ti_old,
+                          e->policy & engine_policy_self_gravity);
 #endif
 
   if (e->verbose)
@@ -2686,7 +2688,8 @@ void engine_prepare(struct engine *e) {
     /* Check that all cells have been drifted to the current time.
      * That can include cells that have not
      * previously been active on this rank. */
-    space_check_drift_point(e->s, e->ti_old);
+    space_check_drift_point(e->s, e->ti_old,
+                            e->policy & engine_policy_self_gravity);
   }
 #endif
 
@@ -2920,7 +2923,8 @@ void engine_print_stats(struct engine *e) {
   /* Check that all cells have been drifted to the current time.
    * That can include cells that have not
    * previously been active on this rank. */
-  space_check_drift_point(e->s, e->ti_current);
+  space_check_drift_point(e->s, e->ti_current,
+                          e->policy & engine_policy_self_gravity);
 
   /* Be verbose about this */
   message("Saving statistics at t=%e.", e->time);
@@ -3219,7 +3223,10 @@ void engine_step(struct engine *e) {
     for (int i = 0; i < e->s->nr_cells; ++i)
       num_gpart_mpole += e->s->cells_top[i].multipole->m_pole.num_gpart;
     if (num_gpart_mpole != e->s->nr_gparts)
-      error("Multipoles don't contain the total number of gpart");
+      error(
+          "Multipoles don't contain the total number of gpart mpoles=%zd "
+          "ngparts=%zd",
+          num_gpart_mpole, e->s->nr_gparts);
   }
 #endif
 
@@ -3228,6 +3235,9 @@ void engine_step(struct engine *e) {
   gravity_exact_force_compute(e->s, e);
 #endif
 
+  /* Do we need to drift the top-level multipoles ? */
+  if (e->policy & engine_policy_self_gravity) engine_drift_top_multipoles(e);
+
   /* Start all the tasks. */
   TIMER_TIC;
   engine_launch(e, e->nr_threads);
@@ -3238,6 +3248,10 @@ void engine_step(struct engine *e) {
   gravity_exact_force_check(e->s, e, 1e-1);
 #endif
 
+  if (!(e->policy & engine_policy_hydro) &&
+      (e->policy & engine_policy_self_gravity) && e->step % 20 == 0)
+    e->forcerebuild = 1;
+
   /* Collect the values of rebuild from all nodes and recover the (integer)
    * end of the next time-step. Do these together to reduce the collective MPI
    * calls per step, but some of the gathered information is not applied just
@@ -3317,8 +3331,8 @@ void engine_unskip(struct engine *e) {
 }
 
 /**
- * @brief Mapper function to drift ALL particle types and multipoles forward in
- * time.
+ * @brief Mapper function to drift *all* particle types and multipoles forward
+ * in time.
  *
  * @param map_data An array of #cell%s.
  * @param num_elements Chunk size.
@@ -3336,7 +3350,7 @@ void engine_do_drift_all_mapper(void *map_data, int num_elements,
       /* Drift all the particles */
       cell_drift_particles(c, e);
 
-      /* Drift the multipole */
+      /* Drift the multipoles */
       if (e->policy & engine_policy_self_gravity)
         cell_drift_all_multipoles(c, e);
     }
@@ -3344,7 +3358,8 @@ void engine_do_drift_all_mapper(void *map_data, int num_elements,
 }
 
 /**
- * @brief Drift *all* particles forward to the current time.
+ * @brief Drift *all* particles and multipoles at all levels
+ * forward to the current time.
  *
  * @param e The #engine.
  */
@@ -3361,7 +3376,54 @@ void engine_drift_all(struct engine *e) {
 
 #ifdef SWIFT_DEBUG_CHECKS
   /* Check that all cells have been drifted to the current time. */
-  space_check_drift_point(e->s, e->ti_current);
+  space_check_drift_point(e->s, e->ti_current,
+                          e->policy & engine_policy_self_gravity);
+#endif
+
+  if (e->verbose)
+    message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+            clocks_getunit());
+}
+
+/**
+ * @brief Mapper function to drift *all* top-level multipoles forward in
+ * time.
+ *
+ * @param map_data An array of #cell%s.
+ * @param num_elements Chunk size.
+ * @param extra_data Pointer to an #engine.
+ */
+void engine_do_drift_top_multipoles_mapper(void *map_data, int num_elements,
+                                           void *extra_data) {
+
+  struct engine *e = (struct engine *)extra_data;
+  struct cell *cells = (struct cell *)map_data;
+
+  for (int ind = 0; ind < num_elements; ind++) {
+    struct cell *c = &cells[ind];
+    if (c != NULL && c->nodeID == e->nodeID) {
+
+      /* Drift the multipole at this level only */
+      cell_drift_multipole(c, e);
+    }
+  }
+}
+
+/**
+ * @brief Drift *all* top-level multipoles forward to the current time.
+ *
+ * @param e The #engine.
+ */
+void engine_drift_top_multipoles(struct engine *e) {
+
+  const ticks tic = getticks();
+
+  threadpool_map(&e->threadpool, engine_do_drift_top_multipoles_mapper,
+                 e->s->cells_top, e->s->nr_cells, sizeof(struct cell), 10, e);
+
+#ifdef SWIFT_DEBUG_CHECKS
+  /* Check that all cells have been drifted to the current time. */
+  space_check_top_multipoles_drift_point(e->s, e->ti_current);
 #endif
 
   if (e->verbose)
@@ -3577,7 +3639,8 @@ void engine_dump_snapshot(struct engine *e) {
   /* Check that all cells have been drifted to the current time.
    * That can include cells that have not
    * previously been active on this rank. */
-  space_check_drift_point(e->s, e->ti_current);
+  space_check_drift_point(e->s, e->ti_current,
+                          e->policy & engine_policy_self_gravity);
 
   /* Be verbose about this */
   message("writing snapshot at t=%e.", e->time);
diff --git a/src/engine.h b/src/engine.h
index 97b978f9a22032a64f5798c8dcf14159f9f9d1af..22ee122bb082895131584942bde50509952b98ff 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -255,6 +255,7 @@ void engine_barrier(struct engine *e, int tid);
 void engine_compute_next_snapshot_time(struct engine *e);
 void engine_unskip(struct engine *e);
 void engine_drift_all(struct engine *e);
+void engine_drift_top_multipoles(struct engine *e);
 void engine_dump_snapshot(struct engine *e);
 void engine_init(struct engine *e, struct space *s,
                  const struct swift_params *params, int nr_nodes, int nodeID,
diff --git a/src/multipole.h b/src/multipole.h
index 0daef6de24039afe63f529ee53f0afda761bc6c3..47f4764daf7402cdebad2e6218fda9adfd5190fe 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -354,6 +354,11 @@ INLINE static void gravity_field_tensors_print(const struct grav_tensor *l) {
 #endif
 }
 
+INLINE static void gravity_multipole_init(struct multipole *m) {
+
+  bzero(m, sizeof(struct multipole));
+}
+
 /**
  * @brief Prints the content of a #multipole to stdout.
  *
diff --git a/src/space.c b/src/space.c
index 432873215c258b987b3c83f87486ade061ea66f0..d9ce5f5a582a47cb79057c9e60787e2a0f64714b 100644
--- a/src/space.c
+++ b/src/space.c
@@ -2150,7 +2150,16 @@ void space_split_recursive(struct space *s, struct cell *c,
     ti_beg_max = get_integer_time_begin(e->ti_current + 1, time_bin_max);
 
     /* Construct the multipole and the centre of mass*/
-    if (s->gravity) gravity_P2M(c->multipole, c->gparts, c->gcount);
+    if (s->gravity) {
+      if (gcount > 0)
+        gravity_P2M(c->multipole, c->gparts, c->gcount);
+      else {
+        gravity_multipole_init(&c->multipole->m_pole);
+        c->multipole->CoM[0] = c->loc[0] + c->width[0] / 2.;
+        c->multipole->CoM[1] = c->loc[1] + c->width[1] / 2.;
+        c->multipole->CoM[2] = c->loc[2] + c->width[2] / 2.;
+      }
+    }
   }
 
   /* Set the values for this cell. */
@@ -2819,11 +2828,26 @@ void space_link_cleanup(struct space *s) {
  *
  * @param s The #space to check.
  * @param ti_drift The (integer) time.
+ * @param multipole Are we also checking the multipoles ?
  */
-void space_check_drift_point(struct space *s, integertime_t ti_drift) {
+void space_check_drift_point(struct space *s, integertime_t ti_drift,
+                             int multipole) {
 #ifdef SWIFT_DEBUG_CHECKS
   /* Recursively check all cells */
-  space_map_cells_pre(s, 1, cell_check_drift_point, &ti_drift);
+  space_map_cells_pre(s, 1, cell_check_particle_drift_point, &ti_drift);
+  if (multipole)
+    space_map_cells_pre(s, 1, cell_check_multipole_drift_point, &ti_drift);
+#else
+  error("Calling debugging code without debugging flag activated.");
+#endif
+}
+
+void space_check_top_multipoles_drift_point(struct space *s,
+                                            integertime_t ti_drift) {
+#ifdef SWIFT_DEBUG_CHECKS
+  for (int i = 0; i < s->nr_cells; ++i) {
+    cell_check_multipole_drift_point(&s->cells_top[i], &ti_drift);
+  }
 #else
   error("Calling debugging code without debugging flag activated.");
 #endif
diff --git a/src/space.h b/src/space.h
index d2879d96b9a4ede4e96236ddd5ac19897fbd10cd..8674c39e110694ec303584627840a7ee47644552 100644
--- a/src/space.h
+++ b/src/space.h
@@ -213,7 +213,10 @@ void space_init_parts(struct space *s);
 void space_init_gparts(struct space *s);
 void space_init_sparts(struct space *s);
 void space_link_cleanup(struct space *s);
-void space_check_drift_point(struct space *s, integertime_t ti_drift);
+void space_check_drift_point(struct space *s, integertime_t ti_drift,
+                             int multipole);
+void space_check_top_multipoles_drift_point(struct space *s,
+                                            integertime_t ti_drift);
 void space_check_timesteps(struct space *s);
 void space_replicate(struct space *s, int replicate, int verbose);
 void space_reset_task_counters(struct space *s);
diff --git a/theory/Multipoles/bibliography.bib b/theory/Multipoles/bibliography.bib
new file mode 100644
index 0000000000000000000000000000000000000000..abdb5855f52aefd9b4562851946d16bd402db68a
--- /dev/null
+++ b/theory/Multipoles/bibliography.bib
@@ -0,0 +1,27 @@
+@ARTICLE{Dehnen2014,
+   author = {{Dehnen}, W.},
+    title = "{A fast multipole method for stellar dynamics}",
+  journal = {Computational Astrophysics and Cosmology},
+archivePrefix = "arXiv",
+   eprint = {1405.2255},
+ primaryClass = "astro-ph.IM",
+     year = 2014,
+    month = sep,
+   volume = 1,
+      eid = {1},
+    pages = {1},
+      doi = {10.1186/s40668-014-0001-7},
+   adsurl = {http://adsabs.harvard.edu/abs/2014ComAC...1....1D},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+@article{Hardy2006,
+  title={Combinatorics of partial derivatives},
+  author={Hardy, Michael},
+  journal={Electron. J. Combin},
+  volume={13},
+   number={1},
+   pages={13},
+   year={2006}
+ }
+	      
diff --git a/theory/Multipoles/fmm_standalone.tex b/theory/Multipoles/fmm_standalone.tex
new file mode 100644
index 0000000000000000000000000000000000000000..6948708622dbf0cfc05daa9d63917f7a3017bfc2
--- /dev/null
+++ b/theory/Multipoles/fmm_standalone.tex
@@ -0,0 +1,22 @@
+\documentclass[fleqn, usenatbib, useAMS,a4paper]{mnras}
+\usepackage{graphicx}
+\usepackage{amsmath,paralist,xcolor,xspace,amssymb}
+\usepackage{times}
+
+\newcommand{\swift}{{\sc Swift}\xspace}
+
+%opening
+\title{FMM in SWIFT}
+\author{Matthieu Schaller}
+
+\begin{document}
+
+\maketitle
+
+\input{gravity_derivatives}
+
+\bibliographystyle{mnras}
+\bibliography{./bibliography.bib}
+
+
+\end{document}
diff --git a/theory/Multipoles/gravity_derivatives.tex b/theory/Multipoles/gravity_derivatives.tex
new file mode 100644
index 0000000000000000000000000000000000000000..913bdf752c666aac5913c6e6452bdb06f3fcdc86
--- /dev/null
+++ b/theory/Multipoles/gravity_derivatives.tex
@@ -0,0 +1,54 @@
+We use the multi-index notation of \cite{Dehnen2014} to simplify expressions.
+
+\subsection{Derivatives of the gravitational potential}
+
+The calculation of all the
+$D_\mathbf{n}(x,y,z) \equiv \nabla^{\mathbf{n}}\phi(x,y,z)$ terms up
+to the relevent order can be quite tedious and it is beneficial to
+automatize the whole setup. Ideally, one would like to have an
+expression for each of this term that is only made of multiplications
+and additions of each of the coordinates and the inverse distance. We
+achieve this by writing $\phi$ as a composition of functions
+$\phi(u(x,y,z))$ and apply the \textit{Fa\`a di Bruno}
+formula \citep[i.e. the ``chain rule'' for higher order derivatives,
+e.g.][]{Hardy2006} to construct our terms:
+
+\begin{equation}
+\label{eq:faa_di_bruno}
+\frac{\partial^n}{\partial x_1 \cdots \partial x_n} \phi(u)
+= \sum_{A} \phi^{(|A|)}(u) \prod_{B \in
+A} \frac{\partial^{|B|}}{\prod_{c\in B}\partial x_c} u(x,y,z),
+\end{equation}
+where $A$ is the set of all partitions of $\lbrace1,\cdots, n\rbrace$,
+$B$ is a block of a partition $A$ and $|\cdot|$ denotes the
+cardinality of a set. For generic functions $\phi$ and $u$ this
+formula yields an untracktable number of terms; an 8th-order
+derivative will have $4140$ (!)  terms in the sum\footnote{The number
+of terms in the sum is given by the Bell number of the same order}. \\
+We choose to write
+\begin{align}
+   \phi(x,y,z) &= 1 / \sqrt{u(x,y,z)}, \\
+   u(x,y,z) &= x^2 + y^2 + z^2.
+\end{align}
+This choice allows to have derivatives of any order of $\phi(u)$ that
+only depend on powers of $u$:
+
+\begin{equation}
+f^{(n)}(u) = \frac{\Gamma(\frac{1}{2})}{\Gamma(\frac{1}{2} -
+n)}\frac{1}{u^{n+\frac{1}{2}}}.
+\end{equation}
+More importantly, this choice of decomposition allows us to have
+derivatives of $u$ only up to second order in $x$, $y$ or $z$. The
+number of non-zero terms in eq. \ref{eq:faa_di_bruno} is hence
+drastically reduced. For instance, when computing
+$D_{(4,1,3)} \equiv \frac{\partial^8}{\partial x^4 \partial y \partial
+z^3} \phi$, $4100$ of the $4140$ terms will involve at least one
+zero-valued derivative (e.g. $\partial^3/\partial x^3$ or
+$\partial^2/\partial x\partial y$) of $u$. Furthermore, among the 40
+remaining terms, many will involve the same derivatives and can be
+grouped together, leaving us with a sum of six products of $x$,$y$ and
+$z$. This is generally the case for most of the $D_\mathbf{n}$'s and
+figuring out which terms are identical in a given set of partitions of
+$\lbrace1,\cdots, n\rbrace$ is an interesting exercise in
+combinatorics left for the reader \citep[see also][]{Hardy2006}.
+
diff --git a/theory/Multipoles/run.sh b/theory/Multipoles/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..34304a37dc1cb1fad3a7442544a8df793ea94d35
--- /dev/null
+++ b/theory/Multipoles/run.sh
@@ -0,0 +1,6 @@
+#!/bin/bash
+python kernels.py
+pdflatex -jobname=fmm fmm_standalone.tex
+bibtex fmm.aux
+pdflatex -jobname=fmm fmm_standalone.tex
+pdflatex -jobname=fmm fmm_standalone.tex