diff --git a/src/cell.c b/src/cell.c
index 5501a259c910bf11bca4b36cfc0b21a6cfb2cffa..906efb2cc894ec7c96d7212ed5c40af49ddf0a05 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -2702,6 +2702,100 @@ void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
   } /* Otherwise, pair interation */
 }
 
+/**
+ * @brief Traverse a sub-cell task and activate the black_holes drift tasks that
+ * are required by a black_holes task
+ *
+ * @param ci The first #cell we recurse in.
+ * @param cj The second #cell we recurse in.
+ * @param s The task #scheduler.
+ */
+void cell_activate_subcell_black_holes_tasks(struct cell *ci, struct cell *cj,
+                                             struct scheduler *s) {
+  const struct engine *e = s->space->e;
+
+  /* Store the current dx_max and h_max values. */
+  ci->black_holes.dx_max_part_old = ci->black_holes.dx_max_part;
+  ci->black_holes.h_max_old = ci->black_holes.h_max;
+  ci->hydro.dx_max_part_old = ci->hydro.dx_max_part;
+  ci->hydro.h_max_old = ci->hydro.h_max;
+
+  if (cj != NULL) {
+    cj->black_holes.dx_max_part_old = cj->black_holes.dx_max_part;
+    cj->black_holes.h_max_old = cj->black_holes.h_max;
+    cj->hydro.dx_max_part_old = cj->hydro.dx_max_part;
+    cj->hydro.h_max_old = cj->hydro.h_max;
+  }
+
+  /* Self interaction? */
+  if (cj == NULL) {
+    /* Do anything? */
+    if (!cell_is_active_black_holes(ci, e) || ci->hydro.count == 0 ||
+        ci->black_holes.count == 0)
+      return;
+
+    /* Recurse? */
+    if (cell_can_recurse_in_self_black_holes_task(ci)) {
+      /* Loop over all progenies and pairs of progenies */
+      for (int j = 0; j < 8; j++) {
+        if (ci->progeny[j] != NULL) {
+          cell_activate_subcell_black_holes_tasks(ci->progeny[j], NULL, s);
+          for (int k = j + 1; k < 8; k++)
+            if (ci->progeny[k] != NULL)
+              cell_activate_subcell_black_holes_tasks(ci->progeny[j],
+                                                      ci->progeny[k], s);
+        }
+      }
+    } else {
+      /* We have reached the bottom of the tree: activate drift */
+      cell_activate_drift_bpart(ci, s);
+      cell_activate_drift_part(ci, s);
+    }
+  }
+
+  /* Otherwise, pair interation */
+  else {
+    /* Should we even bother? */
+    if (!cell_is_active_black_holes(ci, e) &&
+        !cell_is_active_black_holes(cj, e))
+      return;
+
+    /* Get the orientation of the pair. */
+    double shift[3];
+    int sid = space_getsid(s->space, &ci, &cj, shift);
+
+    /* recurse? */
+    if (cell_can_recurse_in_pair_black_holes_task(ci, cj) &&
+        cell_can_recurse_in_pair_black_holes_task(cj, ci)) {
+      struct cell_split_pair *csp = &cell_split_pairs[sid];
+      for (int k = 0; k < csp->count; k++) {
+        const int pid = csp->pairs[k].pid;
+        const int pjd = csp->pairs[k].pjd;
+        if (ci->progeny[pid] != NULL && cj->progeny[pjd] != NULL)
+          cell_activate_subcell_black_holes_tasks(ci->progeny[pid],
+                                                  cj->progeny[pjd], s);
+      }
+    }
+
+    /* Otherwise, activate the sorts and drifts. */
+    else {
+      if (cell_is_active_black_holes(ci, e)) {
+
+        /* Activate the drifts if the cells are local. */
+        if (ci->nodeID == engine_rank) cell_activate_drift_bpart(ci, s);
+        if (cj->nodeID == engine_rank) cell_activate_drift_part(cj, s);
+      }
+
+      if (cell_is_active_black_holes(cj, e)) {
+
+        /* Activate the drifts if the cells are local. */
+        if (ci->nodeID == engine_rank) cell_activate_drift_part(ci, s);
+        if (cj->nodeID == engine_rank) cell_activate_drift_bpart(cj, s);
+      }
+    }
+  } /* Otherwise, pair interation */
+}
+
 /**
  * @brief Traverse a sub-cell task and activate the gravity drift tasks that
  * are required by a self gravity task.
diff --git a/src/cell.h b/src/cell.h
index 18b4512cf1de5efed7e147a9faa52244165aa90a..79599cd40acc4c3e7be1321d18b1260e900da614 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -818,6 +818,8 @@ void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
                                       struct scheduler *s);
 void cell_activate_subcell_stars_tasks(struct cell *ci, struct cell *cj,
                                        struct scheduler *s);
+void cell_activate_subcell_black_holes_tasks(struct cell *ci, struct cell *cj,
+                                             struct scheduler *s);
 void cell_activate_subcell_external_grav_tasks(struct cell *ci,
                                                struct scheduler *s);
 void cell_activate_super_spart_drifts(struct cell *c, struct scheduler *s);
diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c
index ee5f7bfd4fcc7db8f741df2ec419eeb858fbe660..e81d47f8b2b680fbb13a576fee663ed405c8603b 100644
--- a/src/engine_maketasks.c
+++ b/src/engine_maketasks.c
@@ -1622,6 +1622,14 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
                               task_subtype_stars_feedback, flags, 0, ci, cj);
       }
 
+      /* The black hole feedback tasks */
+      if (with_black_holes) {
+        t_bh_density = scheduler_addtask(
+            sched, task_type_pair, task_subtype_bh_density, flags, 0, ci, cj);
+        t_bh_feedback = scheduler_addtask(
+            sched, task_type_pair, task_subtype_bh_feedback, flags, 0, ci, cj);
+      }
+
       engine_addlink(e, &ci->hydro.force, t_force);
       engine_addlink(e, &cj->hydro.force, t_force);
       if (with_limiter) {
@@ -1634,6 +1642,12 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
         engine_addlink(e, &ci->stars.feedback, t_star_feedback);
         engine_addlink(e, &cj->stars.feedback, t_star_feedback);
       }
+      if (with_black_holes) {
+        engine_addlink(e, &ci->black_holes.density, t_bh_density);
+        engine_addlink(e, &cj->black_holes.density, t_bh_density);
+        engine_addlink(e, &ci->black_holes.feedback, t_bh_feedback);
+        engine_addlink(e, &cj->black_holes.feedback, t_bh_feedback);
+      }
 
 #ifdef EXTRA_HYDRO_LOOP
 
@@ -1702,6 +1716,22 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
                               ci->hydro.super->stars.stars_out);
         }
 
+        if (with_black_holes) {
+
+          scheduler_addunlock(sched, ci->hydro.super->black_holes.drift,
+                              t_bh_density);
+          scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
+                              t_bh_density);
+          scheduler_addunlock(
+              sched, ci->hydro.super->black_holes.black_holes_in, t_bh_density);
+          scheduler_addunlock(sched, t_bh_density,
+                              ci->hydro.super->black_holes.ghost);
+          scheduler_addunlock(sched, ci->hydro.super->black_holes.ghost,
+                              t_bh_feedback);
+          scheduler_addunlock(sched, t_bh_feedback,
+                              ci->hydro.super->black_holes.black_holes_out);
+        }
+
         if (with_limiter) {
           scheduler_addunlock(sched, ci->super->kick2, t_limiter);
           scheduler_addunlock(sched, t_limiter, ci->super->timestep);
@@ -1738,6 +1768,23 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
                                 cj->hydro.super->stars.stars_out);
           }
 
+          if (with_black_holes) {
+
+            scheduler_addunlock(sched, cj->hydro.super->black_holes.drift,
+                                t_bh_density);
+            scheduler_addunlock(sched, cj->hydro.super->hydro.drift,
+                                t_bh_density);
+            scheduler_addunlock(sched,
+                                cj->hydro.super->black_holes.black_holes_in,
+                                t_bh_density);
+            scheduler_addunlock(sched, t_bh_density,
+                                cj->hydro.super->black_holes.ghost);
+            scheduler_addunlock(sched, cj->hydro.super->black_holes.ghost,
+                                t_bh_feedback);
+            scheduler_addunlock(sched, t_bh_feedback,
+                                cj->hydro.super->black_holes.black_holes_out);
+          }
+
           if (with_limiter) {
             scheduler_addunlock(sched, cj->super->kick2, t_limiter);
             scheduler_addunlock(sched, t_limiter, cj->super->timestep);
@@ -1780,6 +1827,16 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
                               task_subtype_stars_feedback, flags, 0, ci, NULL);
       }
 
+      /* The black hole feedback tasks */
+      if (with_black_holes) {
+        t_bh_density =
+            scheduler_addtask(sched, task_type_sub_self,
+                              task_subtype_bh_density, flags, 0, ci, NULL);
+        t_bh_feedback =
+            scheduler_addtask(sched, task_type_sub_self,
+                              task_subtype_bh_feedback, flags, 0, ci, NULL);
+      }
+
       /* Add the link between the new loop and the cell */
       engine_addlink(e, &ci->hydro.force, t_force);
       if (with_limiter) {
@@ -1789,6 +1846,10 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
         engine_addlink(e, &ci->stars.density, t_star_density);
         engine_addlink(e, &ci->stars.feedback, t_star_feedback);
       }
+      if (with_black_holes) {
+        engine_addlink(e, &ci->black_holes.density, t_bh_density);
+        engine_addlink(e, &ci->black_holes.feedback, t_bh_feedback);
+      }
 
 #ifdef EXTRA_HYDRO_LOOP
 
@@ -1835,6 +1896,22 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
                             ci->hydro.super->stars.stars_out);
       }
 
+      if (with_black_holes) {
+
+        scheduler_addunlock(sched, ci->hydro.super->black_holes.drift,
+                            t_bh_density);
+        scheduler_addunlock(sched, ci->hydro.super->hydro.drift, t_bh_density);
+        scheduler_addunlock(sched, ci->hydro.super->hydro.sorts, t_bh_density);
+        scheduler_addunlock(sched, ci->hydro.super->black_holes.black_holes_in,
+                            t_bh_density);
+        scheduler_addunlock(sched, t_bh_density,
+                            ci->hydro.super->black_holes.ghost);
+        scheduler_addunlock(sched, ci->hydro.super->black_holes.ghost,
+                            t_bh_feedback);
+        scheduler_addunlock(sched, t_bh_feedback,
+                            ci->hydro.super->black_holes.black_holes_out);
+      }
+
       if (with_limiter) {
         scheduler_addunlock(sched, ci->super->kick2, t_limiter);
         scheduler_addunlock(sched, t_limiter, ci->super->timestep);
@@ -1881,6 +1958,16 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
                               task_subtype_stars_feedback, flags, 0, ci, cj);
       }
 
+      /* The black hole feedback tasks */
+      if (with_feedback) {
+        t_bh_density =
+            scheduler_addtask(sched, task_type_sub_pair,
+                              task_subtype_bh_density, flags, 0, ci, cj);
+        t_bh_feedback =
+            scheduler_addtask(sched, task_type_sub_pair,
+                              task_subtype_bh_feedback, flags, 0, ci, cj);
+      }
+
       engine_addlink(e, &ci->hydro.force, t_force);
       engine_addlink(e, &cj->hydro.force, t_force);
       if (with_limiter) {
@@ -1893,6 +1980,12 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
         engine_addlink(e, &ci->stars.feedback, t_star_feedback);
         engine_addlink(e, &cj->stars.feedback, t_star_feedback);
       }
+      if (with_black_holes) {
+        engine_addlink(e, &ci->black_holes.density, t_bh_density);
+        engine_addlink(e, &cj->black_holes.density, t_bh_density);
+        engine_addlink(e, &ci->black_holes.feedback, t_bh_feedback);
+        engine_addlink(e, &cj->black_holes.feedback, t_bh_feedback);
+      }
 
 #ifdef EXTRA_HYDRO_LOOP
 
@@ -1960,6 +2053,22 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
                               ci->hydro.super->stars.stars_out);
         }
 
+        if (with_black_holes) {
+
+          scheduler_addunlock(sched, ci->hydro.super->black_holes.drift,
+                              t_bh_density);
+          scheduler_addunlock(sched, ci->hydro.super->hydro.drift,
+                              t_bh_density);
+          scheduler_addunlock(
+              sched, ci->hydro.super->black_holes.black_holes_in, t_bh_density);
+          scheduler_addunlock(sched, t_bh_density,
+                              ci->hydro.super->black_holes.ghost);
+          scheduler_addunlock(sched, ci->hydro.super->black_holes.ghost,
+                              t_bh_feedback);
+          scheduler_addunlock(sched, t_bh_feedback,
+                              ci->hydro.super->black_holes.black_holes_out);
+        }
+
         if (with_limiter) {
           scheduler_addunlock(sched, ci->super->kick2, t_limiter);
           scheduler_addunlock(sched, t_limiter, ci->super->timestep);
@@ -1998,6 +2107,23 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
                                 cj->hydro.super->stars.stars_out);
           }
 
+          if (with_feedback) {
+
+            scheduler_addunlock(sched, cj->hydro.super->black_holes.drift,
+                                t_bh_density);
+            scheduler_addunlock(sched, cj->hydro.super->hydro.drift,
+                                t_bh_density);
+            scheduler_addunlock(sched,
+                                cj->hydro.super->black_holes.black_holes_in,
+                                t_bh_density);
+            scheduler_addunlock(sched, t_bh_density,
+                                cj->hydro.super->black_holes.ghost);
+            scheduler_addunlock(sched, cj->hydro.super->black_holes.ghost,
+                                t_bh_feedback);
+            scheduler_addunlock(sched, t_bh_feedback,
+                                cj->hydro.super->black_holes.black_holes_out);
+          }
+
           if (with_limiter) {
             scheduler_addunlock(sched, cj->super->kick2, t_limiter);
             scheduler_addunlock(sched, t_limiter, cj->super->timestep);
diff --git a/src/engine_marktasks.c b/src/engine_marktasks.c
index ba69ab029cdcb29c7da0ed5b61e8f2e23187ede8..11e0bbe10e2531b8b91a7f30c8a3c4d996ebda17 100644
--- a/src/engine_marktasks.c
+++ b/src/engine_marktasks.c
@@ -179,8 +179,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
                t_subtype == task_subtype_bh_density) {
         if (cell_is_active_black_holes(ci, e)) {
           scheduler_activate(s, t);
-          error("TODO!");
-          // cell_activate_subcell_bh_tasks(ci, NULL, s);
+          cell_activate_subcell_black_holes_tasks(ci, NULL, s);
         }
       }
 
@@ -699,6 +698,16 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
       if (cell_is_active_stars(t->ci, e)) scheduler_activate(s, t);
     }
 
+    /* Black hole ghost tasks ? */
+    else if (t_type == task_type_bh_ghost) {
+      if (cell_is_active_black_holes(t->ci, e)) scheduler_activate(s, t);
+    }
+
+    /* Black holes implicit tasks? */
+    else if (t_type == task_type_bh_in || t_type == task_type_bh_out) {
+      if (cell_is_active_black_holes(t->ci, e)) scheduler_activate(s, t);
+    }
+
     /* Time-step? */
     else if (t_type == task_type_timestep) {
       t->ci->hydro.updated = 0;