diff --git a/src/stars/Default/stars.h b/src/stars/Default/stars.h
index a8faf43c1aad4e9c9a300c99be75dede02e85565..67c665c76691321b9f06b84775278464752f7269 100644
--- a/src/stars/Default/stars.h
+++ b/src/stars/Default/stars.h
@@ -56,6 +56,7 @@ __attribute__((always_inline)) INLINE static void stars_init_spart(
     struct spart* sp) {
 
 #ifdef DEBUG_INTERACTIONS_STARS
+  for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i) sp->ids_ngbs_density[i] = -1;
   sp->num_ngb_density = 0;
 #endif
 
@@ -156,6 +157,7 @@ __attribute__((always_inline)) INLINE static void stars_evolve_spart(
 __attribute__((always_inline)) INLINE static void stars_reset_acceleration(
     struct spart* restrict p) {
 #ifdef DEBUG_INTERACTIONS_STARS
+  for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i) p->ids_ngbs_force[i] = -1;
   p->num_ngb_force = 0;
 #endif
 }
diff --git a/src/stars/Default/stars_iact.h b/src/stars/Default/stars_iact.h
index 7e36d4c4af6d121f804b9c43d942ed9c0abc2d0b..09c6f16c7972506a3a762e0a92bf3c85c4c66345 100644
--- a/src/stars/Default/stars_iact.h
+++ b/src/stars/Default/stars_iact.h
@@ -32,6 +32,10 @@ runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
   si->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
 
 #ifdef DEBUG_INTERACTIONS_STARS
+  /* Update ngb counters */
+  if (si->num_ngb_density < MAX_NUM_OF_NEIGHBOURS_STARS)
+    si->ids_ngbs_density[si->num_ngb_density] = pj->id;
+
   /* Update ngb counters */
   ++si->num_ngb_density;
 #endif
@@ -55,6 +59,10 @@ runner_iact_nonsym_stars_feedback(float r2, const float *dx, float hi, float hj,
                                   struct part *restrict pj, float a, float H) {
 
 #ifdef DEBUG_INTERACTIONS_STARS
+  /* Update ngb counters */
+  if (si->num_ngb_force < MAX_NUM_OF_NEIGHBOURS_STARS)
+    si->ids_ngbs_force[si->num_ngb_force] = pj->id;
+
   /* Update ngb counters */
   ++si->num_ngb_force;
 #endif
diff --git a/src/stars/Default/stars_io.h b/src/stars/Default/stars_io.h
index fbed04a00e2ab9b2aad0937b74d92612094a56c1..be849fa290ade621b9077e4189525ca0ead3bb42 100644
--- a/src/stars/Default/stars_io.h
+++ b/src/stars/Default/stars_io.h
@@ -78,12 +78,18 @@ INLINE static void stars_write_particles(const struct spart *sparts,
 #ifdef DEBUG_INTERACTIONS_STARS
 
   list += *num_fields;
-  *num_fields += 2;
+  *num_fields += 4;
 
   list[0] = io_make_output_field("Num_ngb_density", INT, 1, UNIT_CONV_NO_UNITS,
                                  sparts, num_ngb_density);
   list[1] = io_make_output_field("Num_ngb_force", INT, 1, UNIT_CONV_NO_UNITS,
                                  sparts, num_ngb_force);
+  list[2] =
+      io_make_output_field("Ids_ngb_density", LONGLONG, MAX_NUM_OF_NEIGHBOURS_STARS,
+                           UNIT_CONV_NO_UNITS, sparts, ids_ngbs_density);
+  list[3] =
+      io_make_output_field("Ids_ngb_force", LONGLONG, MAX_NUM_OF_NEIGHBOURS_STARS,
+                           UNIT_CONV_NO_UNITS, sparts, ids_ngbs_force);
 #endif
 }
 
diff --git a/src/stars/Default/stars_part.h b/src/stars/Default/stars_part.h
index 2ae3747b8756acddfe071d7f8d52721148b6a328..4f90fd75ae134ca15e4ef2122a0f368464dca6d7 100644
--- a/src/stars/Default/stars_part.h
+++ b/src/stars/Default/stars_part.h
@@ -89,8 +89,14 @@ struct spart {
   /*! Number of interactions in the density SELF and PAIR */
   int num_ngb_density;
 
-  /*! Number of interactions in the force SELF and PAIR */
+  /*! List of interacting particles in the density SELF and PAIR */
+  long long ids_ngbs_density[MAX_NUM_OF_NEIGHBOURS_STARS];
+
+/*! Number of interactions in the force SELF and PAIR */
   int num_ngb_force;
+
+  /*! List of interacting particles in the force SELF and PAIR */
+  long long ids_ngbs_force[MAX_NUM_OF_NEIGHBOURS_STARS];
 #endif
 
 } SWIFT_STRUCT_ALIGN;