diff --git a/src/hydro/SPHENIX/hydro.h b/src/hydro/SPHENIX/hydro.h
index 496a5b8e6b8f13e2ec8b13e4da7ec48b6c02ea62..b0da3ac888ec1a7421ee19471bb5736a01130682 100644
--- a/src/hydro/SPHENIX/hydro.h
+++ b/src/hydro/SPHENIX/hydro.h
@@ -774,6 +774,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
 
   /* Set our old div_v to the one for the next loop */
   p->viscosity.div_v_previous_step = p->viscosity.div_v;
+  p->viscosity.div_v_dt = div_v_dt;
 
   /* Now for the diffusive alpha */
 
@@ -1046,6 +1047,7 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
   p->viscosity.alpha = hydro_props->viscosity.alpha;
   /* Initialise this here to keep all the AV variables together */
   p->viscosity.div_v_previous_step = 0.f;
+  p->viscosity.div_v_dt = 0.f;
 
   /* Set the initial values for the thermal diffusion */
   p->diffusion.alpha = hydro_props->diffusion.alpha;
diff --git a/src/hydro/SPHENIX/hydro_io.h b/src/hydro/SPHENIX/hydro_io.h
index 1b2d31a0ff3d8b244372f538b785109d9f1191a8..ba4eb0f2c581ad42bfc951019d63d7009515194e 100644
--- a/src/hydro/SPHENIX/hydro_io.h
+++ b/src/hydro/SPHENIX/hydro_io.h
@@ -140,7 +140,7 @@ INLINE static void convert_part_potential(const struct engine* e,
 INLINE static void convert_viscosity(const struct engine* e,
                                      const struct part* p,
                                      const struct xpart* xp, float* ret) {
-  ret[0] = p->viscosity.alpha;
+  ret[0] = p->viscosity.alpha * p->force.balsara;
 }
 
 INLINE static void convert_diffusion(const struct engine* e,
@@ -161,7 +161,7 @@ INLINE static void hydro_write_particles(const struct part* parts,
                                          struct io_props* list,
                                          int* num_fields) {
 
-  *num_fields = 11;
+  *num_fields = 14;
   /* List what we want to write */
   list[0] = io_make_output_field_convert_part(
       "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, 1.f, parts, xparts,
@@ -203,11 +203,34 @@ INLINE static void hydro_write_particles(const struct part* parts,
 
   list[9] = io_make_output_field_convert_part(
       "ViscosityParameters", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, parts, xparts,
-      convert_viscosity, "Visosity coefficient (alpha_visc) of the particles");
+      convert_viscosity,
+      "Visosity coefficient (alpha_visc) of the particles, multiplied by the "
+      "balsara switch");
 
   list[10] = io_make_output_field_convert_part(
       "DiffusionParameters", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, parts, xparts,
       convert_diffusion, "Diffusion coefficient (alpha_diff) of the particles");
+
+  list[11] = io_make_output_field(
+      "LaplacianInternalEnergies", FLOAT, 1, UNIT_CONV_FREQUENCY_SQUARED,
+      1.f - 3.f * hydro_gamma, parts, diffusion.laplace_u,
+      "Laplacian (del squared) of the Internal Energy per "
+      "unit mass of the particles");
+
+  list[12] = io_make_output_field(
+      "VelocityDivergences", FLOAT, 1, UNIT_CONV_FREQUENCY, 0.f, parts,
+      viscosity.div_v,
+      "Local velocity divergence field around the particles. Provided without "
+      "cosmology, as this includes the Hubble flow. To return to a peculiar "
+      "velocity divergence, div . v_pec = a^2 (div . v - n_D H)");
+
+  list[13] = io_make_output_field(
+      "VelocityDivergenceTimeDifferentials", FLOAT, 1,
+      UNIT_CONV_FREQUENCY_SQUARED, 0.f, parts, viscosity.div_v_dt,
+      "Time differential (over the previous step) of the "
+      "velocity divergence field around the particles. Again, provided without "
+      "cosmology as this includes a Hubble flow term. To get back to a peculiar
+      velocity divergence time differential, x_pec = a^4 (x - a^{-2} n_D dH / dt)");
 }
 
 /**
diff --git a/src/hydro/SPHENIX/hydro_part.h b/src/hydro/SPHENIX/hydro_part.h
index e7f9a60545d7b0dedc3aa7163aaceb2e8047d234..8d7be455c2d8296b2cee8c0b89da7afa88ac9467 100644
--- a/src/hydro/SPHENIX/hydro_part.h
+++ b/src/hydro/SPHENIX/hydro_part.h
@@ -118,6 +118,9 @@ struct part {
     /*! Particle velocity divergence */
     float div_v;
 
+    /*! Time differential of velocity divergence */
+    float div_v_dt;
+
     /*! Particle velocity divergence from previous step */
     float div_v_previous_step;
 
diff --git a/src/units.c b/src/units.c
index ca9e61a29403fce4963a3907df8a5f2948f5bb6f..e6bba3ac6651d9dc6f8ef0eda5766a90bd3710e2 100644
--- a/src/units.c
+++ b/src/units.c
@@ -258,6 +258,10 @@ void units_get_base_unit_exponents_array(float baseUnitsExp[5],
       baseUnitsExp[UNIT_TIME] = -1.f;
       break;
 
+    case UNIT_CONV_FREQUENCY_SQUARED:
+      baseUnitsExp[UNIT_TIME] = -2.f;
+      break;
+
     case UNIT_CONV_DENSITY:
       baseUnitsExp[UNIT_MASS] = 1.f;
       baseUnitsExp[UNIT_LENGTH] = -3.f;
diff --git a/src/units.h b/src/units.h
index c1a47556df2cc78c872a24fb780ddc321605d099..898221d95b1a74bca1318ec6a2dbf71f8f3a4760 100644
--- a/src/units.h
+++ b/src/units.h
@@ -86,6 +86,7 @@ enum unit_conversion_factor {
   UNIT_CONV_POWER,
   UNIT_CONV_PRESSURE,
   UNIT_CONV_FREQUENCY,
+  UNIT_CONV_FREQUENCY_SQUARED,
   UNIT_CONV_ELECTRIC_CHARGE,
   UNIT_CONV_ELECTRIC_VOLTAGE,
   UNIT_CONV_ELECTRIC_CAPACITANCE,