diff --git a/src/black_holes/EAGLE/black_holes.h b/src/black_holes/EAGLE/black_holes.h
index e9c873f69b40cb2bebfe02355a50a99ab37d3503..ed096303b4bf95b0e033ce4b93ce6e56cc687b53 100644
--- a/src/black_holes/EAGLE/black_holes.h
+++ b/src/black_holes/EAGLE/black_holes.h
@@ -54,8 +54,9 @@ __attribute__((always_inline)) INLINE static void black_holes_first_init_bpart(
 
   bp->time_bin = 0;
   bp->subgrid_mass = bp->mass;
-  bp->energy_reservoir = 0.;
-  bp->formation_time = -1.;
+  bp->energy_reservoir = 0.f;
+  bp->accretion_rate = 0.f;
+  bp->formation_time = -1.f;
 }
 
 /**
@@ -164,7 +165,14 @@ black_holes_bpart_has_no_neighbours(struct bpart* restrict bp,
 }
 
 /**
+ * @brief Compute the accretion rate of the black hole and all the quantites
+ * required for the feedback loop.
  *
+ * @param bp The black hole particle.
+ * @param props The properties of the black hole scheme.
+ * @param constants The physical constants (in internal units).
+ * @param cosmo The cosmological model.
+ * @param dt The time-step size (in physical internal units).
  */
 __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
     struct bpart* restrict bp, const struct black_holes_props* props,
@@ -219,6 +227,7 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
 
   /* Limit the accretion rate to the Eddington fraction */
   const double accr_rate = max(Bondi_rate, f_Edd * Eddington_rate);
+  bp->accretion_rate = accr_rate;
 
   /* Factor in the radiative efficiency */
   const double mass_rate = (1. - epsilon_r) * accr_rate;
diff --git a/src/black_holes/EAGLE/black_holes_io.h b/src/black_holes/EAGLE/black_holes_io.h
index f794d7a298838eee8d6da6fd6ca75d59c8110677..869927527252e368a1d72df0bdf6c1013e007dd6 100644
--- a/src/black_holes/EAGLE/black_holes_io.h
+++ b/src/black_holes/EAGLE/black_holes_io.h
@@ -61,7 +61,7 @@ INLINE static void black_holes_write_particles(const struct bpart *bparts,
                                                int *num_fields) {
 
   /* Say how much we want to write */
-  *num_fields = 10;
+  *num_fields = 11;
 
   /* List what we want to write */
   list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
@@ -84,6 +84,9 @@ INLINE static void black_holes_write_particles(const struct bpart *bparts,
                                  bparts, sound_speed_gas);
   list[9] = io_make_output_field("EnergyReservoir", FLOAT, 1, UNIT_CONV_ENERGY,
                                  bparts, energy_reservoir);
+  list[10] =
+      io_make_output_field("AccretionRate", FLOAT, 1,
+                           UNIT_CONV_MASS_PER_UNIT_TIME, bparts, accretion_rate);
 
 #ifdef DEBUG_INTERACTIONS_BLACK_HOLES
 
diff --git a/src/black_holes/EAGLE/black_holes_part.h b/src/black_holes/EAGLE/black_holes_part.h
index bbb19936b6d33c2b42279f93d53e5fc31d60491d..7d892e91ec5903c3c01a9fde837eb6677ae65bad 100644
--- a/src/black_holes/EAGLE/black_holes_part.h
+++ b/src/black_holes/EAGLE/black_holes_part.h
@@ -76,6 +76,9 @@ struct bpart {
   /*! Energy reservoir for feedback */
   float energy_reservoir;
 
+  /*! Instantaneous accretion rate */
+  float accretion_rate;
+
   /*! Density of the gas surrounding the black hole. */
   float rho_gas;
 
diff --git a/src/parser.h b/src/parser.h
index 9bbc1abcebd052a1160d4fdb7bbef7b11fd2bba5..87e618503dca567f7cf381475fcc91f1575b3913 100644
--- a/src/parser.h
+++ b/src/parser.h
@@ -32,7 +32,7 @@
 
 /* Some constants. */
 #define PARSER_MAX_LINE_SIZE 256
-#define PARSER_MAX_NO_OF_PARAMS 256
+#define PARSER_MAX_NO_OF_PARAMS 512
 #define PARSER_MAX_NO_OF_SECTIONS 64
 
 /* A parameter in the input file */
diff --git a/src/runner.c b/src/runner.c
index c0939bb0bac4d44d51ab692a7119680b5a91316a..068360302580db33dc43f2b8c92162a598a94260 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -513,6 +513,7 @@ void runner_do_black_holes_ghost(struct runner *r, struct cell *c, int timer) {
 
   struct bpart *restrict bparts = c->black_holes.parts;
   const struct engine *e = r->e;
+  const int with_cosmology = e->policy & engine_policy_cosmology;
   const struct cosmology *cosmo = e->cosmology;
   const float black_holes_h_max = e->hydro_properties->h_max;
   const float black_holes_h_min = e->hydro_properties->h_min;
@@ -525,8 +526,6 @@ void runner_do_black_holes_ghost(struct runner *r, struct cell *c, int timer) {
   /* Running value of the maximal smoothing length */
   double h_max = c->black_holes.h_max;
 
-  double dt = 0.001;
-
   TIMER_TIC;
 
   /* Anything to do here? */
@@ -633,6 +632,19 @@ void runner_do_black_holes_ghost(struct runner *r, struct cell *c, int timer) {
           if (((bp->h >= black_holes_h_max) && (f < 0.f)) ||
               ((bp->h <= black_holes_h_min) && (f > 0.f))) {
 
+            /* Get particle time-step */
+            double dt;
+            if (with_cosmology) {
+              const integertime_t ti_step = get_integer_timestep(bp->time_bin);
+              const integertime_t ti_begin =
+                  get_integer_time_begin(e->ti_current - 1, bp->time_bin);
+
+              dt = cosmology_get_delta_time(e->cosmology, ti_begin,
+                                            ti_begin + ti_step);
+            } else {
+              dt = get_timestep(bp->time_bin, e->time_base);
+            }
+
             /* Compute variables required for the feedback loop */
             black_holes_prepare_feedback(bp, e->black_holes_properties,
                                          e->physical_constants, e->cosmology,
@@ -732,6 +744,19 @@ void runner_do_black_holes_ghost(struct runner *r, struct cell *c, int timer) {
         /* Check if h_max has increased */
         h_max = max(h_max, bp->h);
 
+        /* Get particle time-step */
+        double dt;
+        if (with_cosmology) {
+          const integertime_t ti_step = get_integer_timestep(bp->time_bin);
+          const integertime_t ti_begin =
+              get_integer_time_begin(e->ti_current - 1, bp->time_bin);
+
+          dt = cosmology_get_delta_time(e->cosmology, ti_begin,
+                                        ti_begin + ti_step);
+        } else {
+          dt = get_timestep(bp->time_bin, e->time_base);
+        }
+
         /* Compute variables required for the feedback loop */
         black_holes_prepare_feedback(bp, e->black_holes_properties,
                                      e->physical_constants, e->cosmology, dt);
diff --git a/src/units.c b/src/units.c
index 1194735aafd80d51897e4c4cb3e4b0976478145a..55a1c9e1bdf62fce97b22189b4741be0fea2abb4 100644
--- a/src/units.c
+++ b/src/units.c
@@ -377,6 +377,7 @@ void units_get_base_unit_exponants_array(float baseUnitsExp[5],
       break;
 
     case UNIT_CONV_SFR:
+    case UNIT_CONV_MASS_PER_UNIT_TIME:
       baseUnitsExp[UNIT_MASS] = 1.f;
       baseUnitsExp[UNIT_TIME] = -1.f;
       break;
diff --git a/src/units.h b/src/units.h
index 62669425e52c4e39800330a4150259856d8fc0bb..cc87696ccfdeb40b7cf4f23ccfcd853f8c8c1e56 100644
--- a/src/units.h
+++ b/src/units.h
@@ -96,7 +96,8 @@ enum unit_conversion_factor {
   UNIT_CONV_VOLUME,
   UNIT_CONV_INV_VOLUME,
   UNIT_CONV_SFR,
-  UNIT_CONV_SSFR
+  UNIT_CONV_SSFR,
+  UNIT_CONV_MASS_PER_UNIT_TIME
 };
 
 void units_init_cgs(struct unit_system*);